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 thus 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}}$$
The median of positive reals commutes with inversion, it's not clear that this property translates to positive definite matrices with this definition, but it'd be way cool.
