# R/gr_mean.R In RiemGrassmann: Inference, Learning, and Optimization on Grassmann Manifold

#### Documented in gr.mean

#' Fréchet Mean on Grassmann Manifold
#'
#' For manifold-valued data, Fréchet mean is the solution of following cost function,
#' \deqn{\textrm{min}_x \sum_{i=1}^n \rho^2 (x, x_i),\quad x\in\mathcal{M}}
#' for a given data \eqn{\{x_i\}_{i=1}^n} and \eqn{\rho(x,y)} is the geodesic distance
#' between two points on manifold \eqn{\mathcal{M}}. It uses a gradient descent method
#' with a backtracking search rule for updating.
#'
#' @param x either an array of size \eqn{(n\times k\times N)} or a list of length \eqn{N} whose elements are \eqn{(n\times k)} orthonormal basis (ONB) on Grassmann manifold.
#' @param type type of geometry, either \code{"intrinsic"} or \code{"extrinsic"}.
#' @param eps stopping criterion for the norm of gradient.
#' @param parallel a flag for enabling parallel computation with OpenMP.
#'
#' @return a named list containing
#' \describe{
#' \item{mu}{an estimated mean matrix for ONB of size \eqn{(n\times k)}.}
#' \item{variation}{Fréchet variation with the estimated mean.}
#' }
#'
#' @examples
#' ## 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
#' par(mfrow=c(1,2))
#' image(mean.int$mu, main="intrinsic mean") #' image(mean.ext$mu, main="extrinsic mean")
#' par(opar)
#'
#' @author Kisung You
#' @export
gr.mean <- function(x, type=c("intrinsic","extrinsic"), eps=1e-6, parallel=FALSE){
############################################################
# Preprocessing
x      = RiemBase::riemfactory(return_gr(x), name="grassmann")
mytype = match.arg(type)
myeps      = as.double(eps)
myparallel = as.logical(parallel)

############################################################
# Computation
output = RiemBaseExt::rstat.frechet(x, type=mytype, int.eps=myeps, parallel=myparallel)
return(output)
}


## Try the RiemGrassmann package in your browser

Any scripts or data that you put into this service are public.

RiemGrassmann documentation built on March 25, 2020, 5:07 p.m.