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.

Not sure if anything smart can be said about this... thanks for any help.