RiemGrassmann (version 0.1.0)

gr.mean: Fr<U+00E9>chet Mean on Grassmann Manifold

Description

For manifold-valued data, Fr<U+00E9>chet mean is the solution of following cost function, $$\textrm{min}_x \sum_{i=1}^n \rho^2 (x, x_i),\quad x\in\mathcal{M}$$ for a given data \(\{x_i\}_{i=1}^n\) and \(\rho(x,y)\) is the geodesic distance between two points on manifold \(\mathcal{M}\). It uses a gradient descent method with a backtracking search rule for updating.

Usage

gr.mean(x, type = c("intrinsic", "extrinsic"), eps = 1e-06, parallel = FALSE)

Arguments

x

either an array of size \((n\times k\times N)\) or a list of length \(N\) whose elements are \((n\times k)\) orthonormal basis (ONB) on Grassmann manifold.

type

type of geometry, either "intrinsic" or "extrinsic".

eps

stopping criterion for the norm of gradient.

parallel

a flag for enabling parallel computation with OpenMP.

Value

a named list containing

mu

an estimated mean matrix for ONB of size \((n\times k)\).

variation

Fr<U+00E9>chet variation with the estimated mean.

Examples

# NOT RUN {
## generate a dataset with two types of Grassmann elements
#  first four columns of (8x8) identity matrix + noise
mydata = list()
sdval  = 0.1
diag8  = diag(8)
for (i in 1:10){
  mydata[[i]] = qr.Q(qr(diag8[,1:4] + matrix(rnorm(8*4,sd=sdval),ncol=4)))
}

## compute two types of means
mean.int = gr.mean(mydata, type="intrinsic")
mean.ext = gr.mean(mydata, type="extrinsic")

## visualize
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2))
image(mean.int$mu, main="intrinsic mean")
image(mean.ext$mu, main="extrinsic mean")
par(opar)

# }