$ \DeclareMathOperator{\arccosh}{arccosh} \DeclareMathOperator*{\argmin}{arg\,min} \DeclareMathOperator{\Exp}{Exp} \newcommand{\geo}[2]{\gamma_{\overset{\frown}{#1,#2}}} \newcommand{\geoS}{\gamma} \newcommand{\geoD}[2]{\gamma_} \newcommand{\geoL}[2]{\gamma(#2; #1)} \newcommand{\gradM}{\nabla_{\M}} \newcommand{\gradMComp}[1]{\nabla_{\M,#1}} \newcommand{\Grid}{\mathcal G} \DeclareMathOperator{\Log}{Log} \newcommand{\M}{\mathcal M} \newcommand{\N}{\mathcal N} \newcommand{\mat}[1]{\mathbf{#1}} \DeclareMathOperator{\prox}{prox} \newcommand{\PT}[3]{\mathrm{PT}_{#1\to#2}#3} \newcommand{\R}{\mathbb R} \newcommand{\SPD}[1]{\mathcal{P}(#1)} \DeclareMathOperator{\Tr}{Tr} \newcommand{\tT}{\mathrm{T}} \newcommand{\vect}[1]{\mathbf{#1}} $

The Riemannian median

The Riemannian median can be defined similarly to the mean, the Riemannian center of mass. We compute for a matrix of m x n points m medians in parallel, each consisting of n points. The median of one such set is defined as a (not necessarily unique) minimizer of

For this implementation one can further compute a weighted mean, by introducing the summands . The stopping criterion is a functional based on the currend iterate, the last iterate and the number of iterations, where usually a maximum iterations number and a minimial change (of ) is used (see stopCritMaxIterEpsilonCreator).

The algorithm follows [1,2] using a subgradient descent algorithm.

Matlab Documentation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
% median(x) calculates the m medians of x ([.,m,n])
% of the input data with a gradient descent algorithm.
% This implementation is based on
%
% B. Afsari, Riemannian Lp center of mass: Existence,
%    uniqueness, and convexity,
%    Proc. AMS 139(2), pp.655-673, 2011.
% and adapted to the median defined in
% P. T. Fletcher, S. Venkatasubramanian, and S. Joshi:
%    The geometric median on Riemannian manifolds with
%    application to robust atlas estimation.
%    NeuroImage. 45 S143?S152
%
% INPUT
%    x :  m x n Data points ([this.Itemsize,m,n]) to compute
%         m means of n points each, pp.
% OUTPUT
%    y :  m data points of the medians calculated
%
% OPTIONAL
% 'Weights' : (1/n*ones([m,n]) 1xn or mxn weights for the mean
%            the first case uses the same weights for all means
% 'InitVal' : m Initial Data points for the gradient descent
% 'MaxIterations': (50) Maximal Number of Iterations
% 'Epsilon'      : (10^(-5)) Maximal change before stopping
% ---
% Manifold-valued Image Restoration Toolbox 1.0
% J. Persch, R. Bergmann 2015-07-24 | 2018-02-17

See also

References

  1. Afsari, B (2011). Riemannian \(L^p\)center of mass: Existence, uniqueness, and convexity. Proceedings of the American Mathematical Society. 139 655–73
  2. Fletcher, P T, Venkatasubramanian, S and Joshi, S (2009). The geometric meadian on Riemannian manifolds with application to robust atlas estimation. NeuroImage. 45 S143–S152