RiemStiefel (version 0.1.1)

st.mean: Fr<U+00E9>chet Mean on Stiefel 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. In the Stiefel manifold case, analytic formula is not known so we use numerical approximation scheme.

Usage

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

Arguments

x

either an array of size \((p\times r\times n)\) or a list of length \(n\) whose elements are \((p\times r)\) matrix on Stiefel manifold.

type

type of distance, 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 of size \((p\times r)\).

variation

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

References

bhattacharya_nonparametric_2012RiemStiefel

Examples

# NOT RUN {
#-------------------------------------------------------------------
#              Average Projection with 'iris' dataset
#-------------------------------------------------------------------
#  For PCA, take half of data from 'iris' data and repeat it 10 times.
#  We will compare naive PCA, intrinsic, and extrinsic mean.

data(iris)
label = iris$Species           # label information
idata = as.matrix(iris[,1:4])  # numeric data
ndata = nrow(idata)

# define a function for extracting pca projection
pcaproj <- function(X, p){
  return(eigen(stats::cov(X))$vectors[,1:p])
}

# extract embedding for random samples
proj10 = list()
for (i in 1:10){
   # index for random subsample
   rand.now    = base::sample(1:ndata, round(ndata/2)) 
   
   # PCA via personal tool
   proj10[[i]] = pcaproj(idata[rand.now,], p=2)
}

# compute intrinsic and extrinsic mean of projection matrices
mean.int = st.mean(proj10, type='intrinsic')$mu
mean.ext = st.mean(proj10, type='extrinsic')$mu

# compute 2-dimensional embeddings
fproj = pcaproj(idata, p=2)
f2 = idata%*%fproj     # PCA with full data
i2 = idata%*%mean.int  # projection with intrinsic mean
e2 = idata%*%mean.ext  # projection with extrinsic mean

# visualize
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3))
plot(f2, cex=0.5, pch=19, col=label, main='full PCA')
plot(i2, cex=0.5, pch=19, col=label, main='intrinsic projection')
plot(e2, cex=0.5, pch=19, col=label, main='extrinsic projection')
par(opar)

# }