pdMedian {pdSpecEst}

Weighted intrinsic median of HPD matrices

Description

pdMedian calculates a weighted intrinsic median of a sample of (d,d)-dimensional HPD matrices based on a Weiszfeld algorithm intrinsic to the chosen metric.In the case of the affine-invariant Riemannian metric as detailed in e.g., (Bhatia 2009)[Chapter 6] or (Pennec et al. 2006), the intrinsic Weiszfeld algorithm in (Fletcher et al. 2009) is used. By default, the unweighted intrinsic median is computed.

Usage

pdMedian(M, w, metric = "Riemannian", maxit = 1000, reltol)

Arguments

M

a (d,d,S)-dimensional array corresponding to a sample of (d,d)-dimensional HPD matrices of size S.

w

an S-dimensional nonnegative weight vector, such that sum(w) = 1.

metric

the distance measure, one of 'Riemannian', 'logEuclidean', 'Cholesky', 'Euclidean' or 'rootEuclidean'. Defaults to 'Riemannian'.

maxit

maximum number of iterations in gradient descent algorithm. Defaults to 1000

reltol

optional tolerance parameter in gradient descent algorithm. Defaults to 1E-10.

Note

The function does not check for positive definiteness of the input matrices, and (depending on the specified metric) may fail if matrices are close to being singular.

References

Bhatia R (2009). Positive Definite Matrices. Princeton University Press, New Jersey.

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

Pennec X, Fillard P, Ayache N (2006). “A Riemannian framework for tensor computing.” International Journal of Computer Vision, 66(1), 41–66.

See Also

pdMean

Examples

## Generate random sample of HPD matrices
m <- function(){
 X <- matrix(complex(real=rnorm(9), imaginary=rnorm(9)), nrow=3)
 t(Conj(X)) %*% X
}
M <- replicate(100, m())
## Generate random weight vector
z <- rnorm(100)
w <- abs(z)/sum(abs(z))
## Compute weighted intrinsic (Riemannian) median
pdMedian(M, w)