spd.wassbary {Riemann}

Wasserstein Barycenter of SPD Matrices

Description

Given N observations X_1, X_2, \ldots, X_N in SPD manifold, compute the L_2-Wasserstein barycenter that minimizes

\sum_{n=1}^N \lambda_i \mathcal{W}_2 (N(X), N(X_i))^2

where N(X) denotes the zero-mean Gaussian measure with covariance X.

Usage

spd.wassbary(spdobj, weight = NULL, method = c("RU02", "AE16"), ...)

Arguments

spdobj

a S3 "riemdata" class of SPD-valued data of (p\times p) matrices.

weight

weight of observations; if NULL it assumes equal weights, or a nonnegative length-N vector that sums to 1 should be given.

method

name of the algortihm to be used; one of the "RU02", "AE16".

...

extra parameters including

maxiter

maximum number of iterations to be run (default:20).

abstol

tolerance level for stopping criterion (default: 1e-8).

Value

a (p\times p) Wasserstein barycenter matrix.

Examples

#-------------------------------------------------------------------
#        Covariances from standard multivariate Gaussians.
#-------------------------------------------------------------------
## GENERATE DATA
ndata = 20
pdim  = 10
mydat = array(0,c(pdim,pdim,ndata))
for (i in 1:ndata){
  mydat[,,i] = stats::cov(matrix(rnorm(100*pdim), ncol=pdim))
}
myriem = wrap.spd(mydat)

## COMPUTE BY DIFFERENT ALGORITHMS
baryRU <- spd.wassbary(myriem, method="RU02")
baryAE <- spd.wassbary(myriem, method="AE16")

## VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
image(diag(pdim), axes=FALSE, main="True Covariance")
image(baryRU, axes=FALSE, main="by RU02")
image(baryAE, axes=FALSE, main="by AE16")
par(opar)