# median of matrices

I have $n$ positive definite Hermitian matrices $M_n$ and I want to define and compute their median.

These matrices correspond to independent estimations of a covariance matrix in the presence of some noise I cannot quantify, hence the desire to use a median as opposed to a mean (non gaussian residuals / outliers).

I could:

(1) look for a positive definite Hermitian matrix that minimizes $d(M,(M_k)) = \sum \|M-M_k\|_1$.

(2) or I could look a the eigen decompositions $M_k = Q_k \Lambda_k Q_k^T$ (with $\Lambda_k$ sorted) and define $\Lambda = med\ \Lambda_k$ and $Q$ as the orthogonal matrix that minimizes $d(Q,(Q_k))$.

(3) or simply look for the closest (for norm $\|.\|_2$) matrix to the matrix of the element-wise medians.

Kjetil's idea of using the geometric median (=minimizer of the sum of distances) can also be combined with the classical Riemannian distance on positive definite matrices (see e.g. Chapter 6 of Bhatia, Positive Definite Matrices). There are several articles on computing the mean (=minimizer of the sum of squared distances) with respect to this distance; as far as I know the median is unexplored (although the same ideas could probably work, it's mainly optimization on manifolds through quasi-Newton methods, so it doesn't really matter what is the objective function).

With respect to the Euclidean metric, this choice has the advantage that you do not need to do any tricks to ensure positive-definiteness.

• I realize my question is ill-defined, but is there a strong reason to prefer the Riemannian distance on positive definite matrices (which is just the Euclidian distance on their logs if I understand well) to any other distance? If I were to look at the case of dimension 1, assuming that I have estimated the variance of a random variable using 2 different samples, producing estimators $\sigma_1^2$ and $\sigma_2^2$, I would combine them by taking the arithmetic mean, not the geometric mean. Am I missing something? – Bernard Apr 4 '12 at 15:43
• Good question. In some applications it makes more sense, especially when you are estimating matrices whose inverse has also a "meaning". It depends on how your errors are distributed I think. The main advantage that I envision, as I said, is that positive definiteness is nicely embedded in the model and you do not have to do anything artificial to preserve it. – Federico Poloni Apr 4 '12 at 20:32
• And if you want to calculate the geometric median using optimization, remember that unconstrained optimization is easier. So you want to find a good parametrization to use for posdef matrices. See the paper: "Unconstrained Parameterizations for Variance-Covariance Matrices " by Jose C. Pinheiro and Douglas M. Bates Google will find that, it is open-access. – kjetil b halvorsen Apr 5 '12 at 21:19

Here are some relevant Wikipedia articles: http://en.wikipedia.org/wiki/Frechet_mean

The generalization of this to a median is sometimes called the geometric median:

http://en.wikipedia.org/wiki/Geometric_median

Maybe this can point in a/the right direction:

http://www.pnas.org/content/97/4/1423.full.pdf

I think they do medians of vectors there but it's a start.

It's intuitively desirable for the answer not to depend on a unitary transform of the matrix. To estimate the distance of our estimate to the other matrices, a natural choice is the Kullback-Leibler divergence. The equivalent of a mean is then to pick:

$$\hat{\Sigma} = \text{argmin} \left( \sum _{k=1}^{n} \text{tr}\left(\Sigma^{-1}\Sigma_k\right)-\lg \left(\left|\Sigma^{-1}\Sigma_k\right|\right)-d\right)$$

Matrix calculus actually tells us that $$\hat{\Sigma} = \frac{1}{n}\sum _{k=1}^n \Sigma_k$$ too see why differentiate with respect to $\Sigma^{-1}$

Handwaving follows:

In a way, the KL-divergence plays the role of the squared distance here, since the average matrix minimizes the average KL-divergence. Note that this is similar to the Riemann metric, but instead of looking at $\sum_i \lg{(\lambda_i)}^2$ we're looking at $\sum_i \lambda_i-\lg{(\lambda_i)}-1$. If the matrix are contained in a small ball, the $\lambda_i$ are close to $1$ and the difference between the two functions - up to a scaling factor - is $O((\lambda_i-1)^3)$. The KL-divergence has a probabilistic interpretation which isn't clear with the Riemann metric.

We could get a median by using the square root of the KL-divergence.

$$\hat{\Sigma} = \text{argmin} \left( \sum _{k=1}^{n} \sqrt{ \text{tr}\left(\Sigma^{-1}\Sigma_k\right)-\lg \left(\left|\Sigma^{-1}\Sigma_k\right|\right)-d}\right)$$

It's easy to compute iteratively since

$$\frac{df}{d\Sigma} = \sum _{k=1}^{n} \frac{(\mathbf{I} - \Sigma^{-1}\Sigma_k)\Sigma^{-1}}{2\sqrt{ \text{tr}\left(\Sigma^{-1}\Sigma_k\right)-\lg \left(\left|\Sigma^{-1}\Sigma_k\right|\right)-d}}$$

Positive definite Hermitian matrices form a Cartan-Hadamard manifold the Riemannian distance and the explicite expressions of geodesics in this manifold are known (see R. Bhatia's book "Positive definite matrices"), so that your problem can be solved by the subgradient algorithm in the following article:

Riemannian median and its estimation, LMS J.Comput. Math. 13(2010), 461-479

Moreover, stochatic algorithms for finding p-means (p=1 for medians and p=2 for barycenters) in Riemannian manifolds can be found in the article:

Stochastic algorithms for computing means of probability measures, Stochastic Processes and their Applications, Volume 122, Issue 4, April 2012, Pages 1437–1455

Deterministic algorithms for computing p-means can be found in Chpter 4 of the thesis (French title but English contents):

http://tel.archives-ouvertes.fr/tel-00664188