rbase.median {RiemBase}

Geometric Median of Manifold-valued Data

Description

For manifold-valued data, geometric median is the solution of following cost function,

\textrm{min}_x \sum_{i=1}^n \rho (x, x_i) = \sum_{i=1}^n \| \log_x (x_i) \|,\quad x\in\mathcal{M}

for a given data \{x_i\}_{i=1}^n, \rho(x,y) the geodesic distance between two points on manifold \mathcal{M}, and \| \log_x (y) \| a logarithmic mapping onto the tangent space T_x \mathcal{M}. Weiszfeld's algorithms is employed.

Usage

rbase.median(input, maxiter = 496, eps = 1e-06, parallel = FALSE)

Arguments

input

a S3 object of riemdata class. See riemfactory for more details.

maxiter

maximum number of iterations for gradient descent algorithm.

eps

stopping criterion for the norm of gradient.

parallel

a flag for enabling parallel computation.

Value

a named list containing

x

an estimate geometric median.

iteration

number of iterations until convergence.

Author(s)

Kisung You

References

Fletcher PT, Venkatasubramanian S, Joshi S (2009). “The geometric median on Riemannian manifolds with application to robust atlas estimation.” NeuroImage, 45(1), S143–S152. ISSN 10538119.

Aftab K, Hartley R, Trumpf J (2015). “Generalized Weiszfeld Algorithms for Lq Optimization.” IEEE Transactions on Pattern Analysis and Machine Intelligence, 37(4), 728–745. ISSN 0162-8828, 2160-9292.

Examples

### Generate 100 data points on Sphere S^2 near (0,0,1).
ndata = 100
theta = seq(from=-0.99,to=0.99,length.out=ndata)*pi
tmpx  = cos(theta) + rnorm(ndata,sd=0.1)
tmpy  = sin(theta) + rnorm(ndata,sd=0.1)

### Wrap it as 'riemdata' class
data  = list()
for (i in 1:ndata){
  tgt = c(tmpx[i],tmpy[i],1)
  data[[i]] = tgt/sqrt(sum(tgt^2)) # project onto Sphere
}
data = riemfactory(data, name="sphere")

### Compute Geodesic Median
out1 = rbase.median(data)
out2 = rbase.median(data,parallel=TRUE) # test parallel implementation