Geometric Matrix Midranges

TL;DR

This paper introduces geometric matrix midranges for positive definite Hermitian matrices, exploring their properties, computational aspects, and extensions from two matrices to multiple matrices, with insights from linear algebra, geometry, and optimization.

Contribution

It defines and analyzes geometric matrix midranges, extending the concept from two matrices to multiple matrices, and examines their properties from various mathematical perspectives.

Findings

Comparison of matrix midrange with scalar and vector cases

Significance of matrix midrange in computational contexts

Analysis using linear algebra, differential geometry, and convex optimization

Abstract

We define geometric matrix midranges for positive definite Hermitian matrices and study the midrange problem from a number of perspectives. Special attention is given to the midrange of two positive definite matrices before considering the extension of the problem to $N > 2$ matrices. We compare matrix midrange statistics with the scalar and vector midrange problem and note the special significance of the matrix problem from a computational standpoint. We also study various aspects of geometric matrix midrange statistics from the viewpoint of linear algebra, differential geometry and convex optimization.

Equations

$$x=\frac{1}{2}\left(\underset{i}{\min }{a}_i+\underset{i}{\max }{a}_i\right).$$

$$\underset{x\in \mathbb{R}}{\min }\underset{i}{\max }|x-a_i|.$$

$$x=\exp \left(\frac{1}{2}\left[\underset{i}{\min }\log y_i+\underset{i}{\max }\log y_i\right]\right)={\left(\underset{i}{\min }{y}_i\cdot \underset{i}{\max }{y}_i\right)}^{1/2}.$$

$$\min_{x>0}\;\max_{i}\;|\log x-\log y_{i}|=\min_{x>0}\;\max_{i}\;\bigg{|}\log\frac{x}{y_{i}}\bigg{|}.$$

$$\underset{X\succ 0}{\min }\underset{i}{\max }\Vert \log (Y_i^{-1/2}X{Y}_i^{-1/2}){\Vert }_{\infty },$$

$$\underset{X\succ 0}{\min }\underset{i}{\max }\Vert \log (Y_i^{-1/2}X{Y}_i^{-1/2}){\Vert }_2,$$

$$A\ast B=\frac{1}{\sqrt{{\lambda }_{\min }}+\sqrt{{\lambda }_{\max }}}\left(B+\sqrt{{\lambda }_{\min }{\lambda }_{\max }}A\right),$$

$$A\#B=A^{1/2}{\left(A^{-1/2}B{A}^{-1/2}\right)}^{1/2}{A}^{1/2}.$$

$$d_{\infty }(A,B)=\log \max \{{\lambda }_{\max }(A{B}^{-1}),{\lambda }_{\max }(B{A}^{-1})\}.$$

$$\underset{X\succ 0}{\min }\max \{d_{\infty }(A,X),d_{\infty }(B,X)\}.$$

$$\max \{\log {\lambda }_{\max }(X{A}^{-1}),\log {\lambda }_{\max }(X{B}^{-1}),-\log {\lambda }_{\min }(X{A}^{-1}),-\log {\lambda }_{\min }(X{B}^{-1})\}.$$

$$X{A}^{-1}$$

$$X{B}^{-1}$$

$$B\ast A$$

$$=\frac{\sqrt{{\lambda }_{\min }{\lambda }_{\max }}}{\sqrt{{\lambda }_{\min }}+\sqrt{{\lambda }_{\max }}}\left(\frac{1}{\sqrt{{\lambda }_{\min }{\lambda }_{\max }}}B+A\right)=A\ast B.$$

$${\lambda }_i\left(CB{C}^{\ast }(CA{C}^{\ast }{)}^{-1}\right)={\lambda }_i\left(CB{C}^{\ast }(C^{\ast }{)}^{-1}{A}^{-1}{C}^{-1})\right)={\lambda }_i(B{A}^{-1}).$$

$${\lambda }_i(I\ast \Sigma )-1$$

$$=\frac{{\lambda }_i(\Sigma )+\sqrt{{\lambda }_{\min }{\lambda }_{\max }}}{\sqrt{{\lambda }_{\min }}+\sqrt{{\lambda }_{\max }}}-1$$

$$=\frac{{\delta }_i+\left(\sqrt{{\lambda }_{\min }}-1\right)\left(\sqrt{{\lambda }_{\max }}-1\right)}{\sqrt{{\lambda }_{\min }}+\sqrt{{\lambda }_{\max }}}\ge 0,$$

$${\lambda }_i(\Sigma -I\ast \Sigma )$$

$$=\left(\frac{\sqrt{{\lambda }_{\min }}+\sqrt{{\lambda }_{\max }}-1}{\sqrt{{\lambda }_{\min }}+\sqrt{{\lambda }_{\max }}}\right){\lambda }_i(\Sigma )-\frac{\sqrt{{\lambda }_{\min }{\lambda }_{\max }}}{\sqrt{{\lambda }_{\min }}+\sqrt{{\lambda }_{\max }}}$$

$$=\sqrt{{\lambda }_{\min }}\left(\sqrt{{\lambda }_{\min }}-1\right)+\left(\frac{\sqrt{{\lambda }_{\min }}+\sqrt{{\lambda }_{\max }}-1}{\sqrt{{\lambda }_{\min }}+\sqrt{{\lambda }_{\max }}}\right){ϵ}_i\ge 0,$$

$$I⪯\Sigma ⟹I⪯I\ast \Sigma ⪯\Sigma ,$$

$${\Sigma }_1⪯{\Sigma }_2⟹I\ast {\Sigma }_1⪯I\ast {\Sigma }_2.$$

$${\lambda }_i(I\ast \Sigma )=\frac{{\lambda }_i(\Sigma )+\sqrt{{\lambda }_{\min }{\lambda }_{\max }}}{\sqrt{{\lambda }_{\min }}+\sqrt{{\lambda }_{\max }}},$$

$${\lambda }_{\min }(I\ast \Sigma )=\sqrt{{\lambda }_{\min }(\Sigma )}\mathrm{and}{\lambda }_{\max }(I\ast \Sigma )=\sqrt{{\lambda }_{\max }(\Sigma )},$$

$${\lambda }_2(I\ast \Sigma )=f(x):=\frac{b+\sqrt{ax}}{\sqrt{a}+\sqrt{x}}.$$

$$f^{\prime }(x)=\frac{a-b}{2\sqrt{x}(\sqrt{a}+\sqrt{x}{)}^2}<0,\forall x\ge b,$$

$$(aA)\ast (bB)=\sqrt{ab}(A\ast B).$$

$$A\diamond B=\{\begin{array}{l}\frac{\sqrt{{\lambda }_{\max }}}{1+{\lambda }_{\max }}(A+B)\mathrm{if}{\lambda }_{\min }{\lambda }_{\max }\ge 1\\ \frac{\sqrt{{\lambda }_{\min }}}{1+{\lambda }_{\min }}(A+B)\mathrm{if}{\lambda }_{\min }{\lambda }_{\max }\le 1\end{array}$$

$$d_{\Phi }(A,B)=\Vert \log A^{-1/2}B{A}^{-1/2}{\Vert }_{\Phi },$$

$${\gamma }_{\mathcal{G}}(t)=A^{1/2}{\left(A^{-1/2}B{A}^{-1/2}\right)}^t{A}^{1/2}.$$

$$d_{\infty }(A,B)=\Vert \log A^{-1/2}B{A}^{-1/2}{\Vert }_{\infty }=\max \{\log {\lambda }_{\max }(B{A}^{-1}),\log {\lambda }_{\max }(A{B}^{-1})\}.$$

Full text

{AMS}

15B48, 53C22, 90C26

1 Introduction

The midrange of a collection of real numbers $a_{1},\dots,a_{N}$ is defined as the arithmetic average of the extremal values. That is,

[TABLE]

This is the unique solution to the optimization problem

[TABLE]

In this paper, we are interested in midrange statistics in convex cones and in particular the cone of positive definite Hermitian matrices of a fixed dimension.111This manuscript further develops ideas partially introduced by the authors in [39]. The midrange of scalar-valued data is sensitive to outliers and is therefore a non-robust statistic. Despite this, it can be a useful measure in some contexts. For instance, the midrange is the maximally efficient estimator for the center of a uniform distribution. Thus, it can be an appropriate tool for data that is devoid of extreme outliers. It can also be useful in clustering algorithms that require the isolation of outlying clusters [16, 47, 49]. It is an important notion in the statistics of extreme events [26].

Data representations based on symmetric positive definite matrices are common in a variety of applications from computer vision to machine learning. Often such matrices arise as covariance matrices that encode the correlations implicit in data and are thus highly structured [56]. Specific applications include brain-computer interface (BCI) systems [45, 55], radar data processing [3], and diffusion tensor imaging (DTI) [19]. It has been noted in many works that using nonlinear geometries related to generalized spectral properties of positive definite matrices yield significantly improved performance [44]. Indeed, Euclidean techniques for statistics and analysis on covariance matrices often result in poor accuracy and undesirable effects, such as swelling phenomena in DTI [5]. It is in this context that much attention has been paid to developing geometric statistical methods on the cone of positive definite matrices [2, 8, 9, 29, 38, 46]. A fundamental geometry that is associated to such spaces is the affine-invariant geometry [8, 40], whereby congruence transformations play the role of translations between matrices. The analogue of this geometry for scalars defined in the cone of positive real numbers $\mathbb{R}_{+}=\{x\in\mathbb{R}:x>0\}$ simply reduces to working with the logarithms of the data points and then mapping the result back to the positive cone $\mathbb{R}_{+}$ via the exponential map. Thus, we can define the affine-invariant midrange of $N$ positive numbers $y_{i}>0$ to be

[TABLE]

Note that Eq. 1 is the unique solution of the optimization problem

[TABLE]

In the matrix setting, we define the geometric midrange problem on the cone of positive definite matrices as

[TABLE]

where $\{Y_{1},\dots,Y_{N}\}$ are a collection of $N$ positive definite matrices of dimension $n$, $\|\cdot\|_{\infty}$ denotes the spectral operator norm on the space of Hermitian matrices of dimension $n$ defined by $\|A\|_{\infty}=\max\{|\lambda_{1}(A)|,\cdots,|\lambda_{n}(A)|\}$, and $X\succ 0$ denotes the positive definiteness of $X$. Note that Eq. 2 can be interpreted as the smallest enclosing ball problem for the collection of data $\{Y_{i}\}$ in Thompson geometry. In particular, Eq. 2 can be expressed as $\min_{X\succ 0}\;\max_{i}\;d_{\infty}(X,Y_{i})$, where $d_{\infty}$ denotes the Thompson metric (see Section 2). The smallest enclosing ball problem of a finite set of points in Euclidean space was first posed by Sylvester in [50] and is a fundamental problem in computational geometry. The problem is also known as the minimum enclosing ball, the 1-center, or the minimax optimization problem and has been studied by several authors [6, 54]. It is an important problem that finds many applications in computer graphics and machine learning, including in collision detection, support vector clustering and similarity search [41, 53]. On manifolds, the Riemannian smallest enclosing ball problem has been studied by Arnaudon and Nielsen in [4]. In particular, the authors consider

[TABLE]

which is the corresponding problem with respect to the standard affine-invariant Riemannian geometry of positive definite matrices (see Section 2.3). Note that $\|\cdot\|_{2}$ in Eq. 3 denotes the standard Frobenius norm. Although this problem is clearly closely related to Eq. 2, there are fundamental differences between them. For instance, as proved by Afsari in [1], there exists a unique point that minimizes the cost function in Eq. 3. In contrast, even in the case of two matrices $Y_{1}=A$ and $Y_{2}=B$, the solution to the optimization problem Eq. 2 is generally not unique. In Section 2.1, we provide an interpretation of the Riemannian distance $d_{2}$ and Thompson distance $d_{\infty}$ as members of a family of affine-invariant Finsler distances.

One particular analytic solution for the geometric midrange of two positive definite matrices $A$ and $B$ that will receive special consideration is $A*B$ defined by

[TABLE]

where $\lambda_{\max}$ and $\lambda_{\min}$ denote the maximum and minimum generalized eigenvalues of the pencil $(B,A)$, which are determined by the equation $\det(B-\lambda A)=0$. Note that $\lambda_{\max}$ and $\lambda_{\min}$ also coincide with the maximum and minimum eigenvalues of $BA^{-1}$, respectively. We will consider the properties of the expression in Eq. 4 in some detail in Section 2. For now, it is instructive to compare $A*B$ with the well-known geometric mean $A\#B$ of positive definite matrices $A$ and $B$ given by

[TABLE]

The matrix geometric mean has been studied in great detail by several authors and is used in a variety of applications. Much research has been devoted to extending the notion of a geometric mean from two matrices to an arbitrary number of matrices and finding efficient algorithms for computing such a mean [12, 27, 28]. These include optimization-based approaches [10] as well as inductive sequential constructions [2, 11, 35, 36].

It is noteworthy that the formula Eq. 4 for a matrix midrange of $A$ and $B$ is considerably less expensive to compute than the the geometric mean Eq. 5, particularly for high dimensional matrices. This is because $A*B$ mainly relies on the evaluation of extremal generalized eigenvalues that can be computed efficiently using a variety of techniques such as Krylov subspace methods [23, 24, 37, 48]. In contrast, the Cholesky-Schur algorithm for computing the geometric mean Eq. 5 of two matrices has a complexity of $O(n^{3})$ [27]. Thus, we already see an important difference between the scalar and matrix midrange problems: in the scalar case, the mean and midrange of two points are trivially the same, whereas a geometric midrange of two matrices may be much cheaper to compute than their geometric mean. The plots in Fig. 1 $(a)$ and $(b)$ illustrate how the computational cost of the midrange $A*B$ evolves with the matrix dimension as compared to the arithmetic and geometric means. The computations are based on a large number of randomly generated positive definite matrices of dimension $n=5$ to $n=10000$. Fig. 1 $(c)$ provides a similar comparison for the cost of computing the Thompson distance versus the Euclidean and affine-invariant Riemannian distances. The logarithms in Fig. 1 $(b)$ and $(c)$ refer to the natural logarithm. All computations were performed on a 2017 Apple MacBook Pro laptop in MATLAB. The extremal generalized eigenvalues of $(B,A)$ that appear in $A*B$ were computed using the eigs function in MATLAB with the default settings, which utilizes algorithms outlined in [32, 48].

1.1 Paper organization and contributions

The paper is organized as follows. In Section 2, the midrange of two positive definite matrices is studied in detail from a variety of perspectives. We begin by proving a number of key properties of Eq. 4 that are expected of a measure of central tendancy, including suitable order and monotonicity properties. In Section 2.1, we present an interpretation of the geometric midrange within a unified optimization framework alongside the geometric mean and median. In Section 2.2, we present a characterization of the midrange formula Eq. 4 based on an extremal ordering property defined using the Löwner order. In Section 2.3, we review the differential geometry of the manifold of $n\times n$ positive definite Hermitian matrices and consider midranges arising as midpoints of geodesics. In Section 3, we define the geometric midrange problem for $N$ positive definite matrices and study its properties in some detail. We offer a solution to the problem via convex optimization in Section 3.1 before proving a number of optimality conditions and related results in Section 3.2.

2 Midrange of two positive definite matrices

Let $\mathbb{P}_{n}$ denote the set of $n\times n$ positive definite Hermitian matrices, which is the interior of the pointed, closed and convex cone of positive semidefinite matrices of the same dimensions. A pointed, closed and convex cone $C$ in a vector space $V$ induces a partial order on $V$ given by $x\leq y$ if and only if $y-x\in C$. The Thompson metric [33, 51] on $C$ is defined to be $d_{\infty}(x,y)=\log\max\{M(x/y;C),M(y/x;C)\}$, where $M(y/x;C)=\inf\{\lambda\in\mathbb{R}:y\leq\lambda x\}$ for $x\in C\setminus\{0\}$ and $y\in V$. For $A,B\in\mathbb{P}_{n}$, we have $M(A/B)=\lambda_{\max}(AB^{-1})$, so that

[TABLE]

Noting that $\lambda_{i}(A^{-1/2}BA^{-1/2})=\lambda_{i}(BA^{-1})$ and $\lambda_{\max}(\Sigma^{-1})=1/\lambda_{\min}(\Sigma)$ for any $\Sigma\in\mathbb{P}_{n}$, we find that the 2-point midrange problem Eq. 2 for data $A$ and $B$ takes the form

[TABLE]

A point $X$ is said to be a Thompson midpoint of the pair $(A,B)$ if $d_{\infty}(A,X)=d_{\infty}(B,X)=\frac{1}{2}d_{\infty}(A,B)$. As $\left(\mathbb{P}_{n},d_{\infty}\right)$ forms a complete metric space [33], the minimizers of Eq. 7 coincide with the Thompson midpoints of $(A,B)$, which are generally non-unique. The geometry of the set of Thompson midpoints of a given pair of points $A,B\in\mathbb{P}_{n}$ is studied in detail in [34], where it is shown that the midpoint is unique if and only if the spectrum of $BA^{-1}$ lies in a set $\{\lambda,\lambda^{-1}\}$ for some $\lambda>0$. In this paper, we will pay special attention to the midrange $A*B$ given by Eq. 4 due to its scalable computational properties.

Note that Eq. 7 is equivalent to $\min_{X\succ 0}f(X)$, where $f(X)$ is given by

[TABLE]

Using this expression and the following elementary lemma, it is easy to verify that $A*B$ is indeed a Thompson metric midpoint of $(A,B)$.

Lemma 2.1**.**

*If $c_{1},c_{2}\in\mathbb{R}$ and $M$ is an $n\times n$ matrix with eigenvalues $\lambda_{i}(M)$, then $c_{1}M+c_{2}I$ has eigenvalues $c_{1}\lambda_{i}(M)+c_{2}$. *

Specifically, we find that for $X=A*B$, we have

[TABLE]

where $\lambda_{\max}$ and $\lambda_{\min}$ refer to the extremal eigenvalues of $BA^{-1}$. Using Lemma 2.1 and $\lambda_{\max}(BA^{-1})=1/\lambda_{\min}(AB^{-1})$, we find that $f(A*B)$ simplifies to $\frac{1}{2}d_{\infty}(A,B)$ as required.

We now consider the merits of the midrange $A*B$ as a measure of central tendency for $\{A,B\}$. The following are a number of properties that are desirable for such a mapping $\mu:\mathbb{P}_{n}\times\mathbb{P}_{n}\rightarrow\mathbb{P}_{n}$. We denote the conjugate transpose of $X$ by $X^{*}$ and the general linear group of $n\times n$ matrices by $GL(n)$.

Continuity: $\mu$ is a continuous map. 2. 2.

Symmetry: $\mu(A,B)=\mu(B,A)$ for all $A,B\in\mathbb{P}_{n}$. 3. 3.

Affine-invariance: $\mu(XAX^{*},XBX^{*})=X\mu(A,B)X^{*}$, for all $X\in GL(n)$. 4. 4.

Order property: $A\preceq B\implies A\preceq\mu(A,B)\preceq B$. 5. 5.

Monotonicity: $\mu(A,B)$ is monotone in its arguments.

We will now prove that $\mu(A,B):=A*B$ indeed satisfies properties 1-3 listed above before turning our attention to the order and monotonicity properties 4 and 5, which merit special consideration.

Proposition 2.2**.**

*The map $\mu(A,B)=A*B$ satisfies properties 1-3. *

Proof 2.3**.**

1. The continuity of $\mu$ follows directly from the expression for $A*B$ in Eq. 4, the invertibility of $A$, and the continuous dependence of eigenvalues on matrix entries, which itself follows from consideration of the roots of the characteristic polynomial of a matrix. 2. For symmetry, we note that $\lambda_{\min}(AB^{-1})=1/\lambda_{\max}(BA^{-1})$ and $\lambda_{\max}(AB^{-1})=1/\lambda_{\min}(BA^{-1})$, so that

[TABLE]

3. Affine-invariance follows immediately by noting that

[TABLE]

The order property is a generalization of the property of means of positive numbers whereby a mean of a pair of points is expected to lie between the two points on the number line. For Hermitian matrices, a standard partial order $\preceq$ exists according to which $A\preceq B$ if and only if $B-A$ is positive semidefinite. This partial order is known as the Löwner order and the monotonicity in condition 5 is also with reference to this order. Unlike the case of real positive numbers $a,b>0$, which always satisfy $a\leq b$ or $b\leq a$, two Hermitian matrices $A$ and $B$ may fail to satisfy both $A\preceq B$ and $B\preceq A$. It is well-known that the Löwner order is affine-invariant in the sense that for all $A,B\in\mathbb{P}_{n}$, $X\in GL(n)$, $A\preceq B$ implies that $XAX^{*}\preceq XBX^{*}$. In particular, $A\preceq B$ if and only if $I\preceq A^{-1/2}BA^{-1/2}$. Thus, by affine-invariance of $\mu$, it suffices to prove point 4 in the case where $A=I$ since $A\preceq\mu(A,B)\preceq B$ if and only if $I\preceq\mu(I,A^{-1/2}BA^{-1/2})\preceq A^{-1/2}BA^{-1/2}$. To establish the 4th property for $\mu(A,B)=A*B$, we make use of Lemma 2.1. Let $\Sigma\in\mathbb{P}_{n}$ be such that $I\preceq\Sigma$ and note that this is equivalent to $\lambda_{i}(\Sigma)\geq 1$ for $i=1,\ldots,n$. Writing $\lambda_{\min}=\lambda_{\min}(\Sigma)$, $\lambda_{\max}=\lambda_{\max}(\Sigma)$, and $\lambda_{i}(\Sigma)=1+\delta_{i}$ for $\delta_{i}\geq 0$, we have by Lemma 2.1 that

[TABLE]

since $\lambda_{i}(\Sigma)\geq 1$ implies that $\sqrt{\lambda_{i}(\Sigma)}\geq 1$. Thus, we have shown that $I\preceq\Sigma$ implies $I\preceq I*\Sigma$. To prove the other inequality, let $\lambda_{i}(\Sigma)=\lambda_{\min}(\Sigma)+\epsilon_{i}$ for $\epsilon_{i}\geq 0$, and note that

[TABLE]

as $I\preceq\Sigma$ ensures that $\sqrt{\lambda_{\min}}\geq 1$. Therefore, we have also shown that $\Sigma-I*\Sigma\succeq 0$. That is,

[TABLE]

for all $\Sigma\in\mathbb{P}_{n}$. In particular, upon substituting $\Sigma=A^{-1/2}BA^{-1/2}$ in Eq. 8 and using the affine-invariance properties of both the Löwner order and the mean $\mu(A,B)=A*B$, we establish the following important property.

Proposition 2.4**.**

*For $A,B\in\mathbb{P}_{n}$, $A\preceq B$ implies that $A\preceq A*B\preceq B$. *

We now consider the 5th and final desirable property of $\mu:\mathbb{P}_{n}\times\mathbb{P}_{n}\rightarrow\mathbb{P}_{n}$, which is monotonicity of $\mu$ in its arguments. First recall that a map $F:\mathbb{P}_{n}\rightarrow\mathbb{P}_{n}$ is said to be monotone if $\Sigma_{1}\preceq\Sigma_{2}$ implies that $F(\Sigma_{1})\preceq F(\Sigma_{2})$. By symmetry and affine-invariance, it is sufficient to consider monotonicity of $\mu(I,\Sigma)$ with respect to $\Sigma$. That is, monotonicity is established by showing that

[TABLE]

However, it turns out that $F(\Sigma):=I*\Sigma$ is not monotone with respect to $\Sigma$ as we demonstrate below. Nonetheless, $F$ is seen to enjoy certain weaker monotonicity properties. Considering the eigenvalues of $I*\Sigma$, we find that

[TABLE]

where $\lambda_{\min}$ and $\lambda_{\max}$ refer to the smallest and largest eigenvalues of $\Sigma$.

Proposition 2.5**.**

*The maximum and minimum eigenvalues of $F(\Sigma)=I*\Sigma$ are monotone with respect to $\Sigma$. *

Proof 2.6**.**

Considering the cases $i=1$ and $i=n$, we find that (9) yields

[TABLE]

*both of which are seen to be monotone functions of $\Sigma$. *

It is in the sense of the above that $\mu(A,B)=A*B$ inherits a weak monotonicity property. The monotonic dependence of the extremal eigenvalues of $I*\Sigma$ on $\Sigma$ ensures that if $\Sigma_{1}\preceq\Sigma_{2}$, then we can at least rule out the possibility that $I*\Sigma_{1}\succ I*\Sigma_{2}$, where $\succ 0$ here denotes positive definiteness. To prove that monotonicity is generally not satisfied in the full sense, consider a diagonal matrix $\Sigma=\operatorname{diag}(a,b,x)\in\mathbb{P}(3)$, where $\lambda_{\min}(\Sigma)=a<b\leq x=\lambda_{\max}(\Sigma)$ and $x$ is thought of as a variable. We have $I*\Sigma=\operatorname{diag}\left(\sqrt{a},f(x),\sqrt{x}\right)$, where

[TABLE]

Taking the derivative of $f$ with respect to $x$, we find that

[TABLE]

which shows that the second eigenvalue of $I*\Sigma$ decreases as $x$ increases. Thus, we see that $I*\Sigma$ cannot depend monotonically on $\Sigma$ in this example.

As a summary, we collect the main results so far in the following theorem.

Theorem 2.7**.**

*The midrange $\mu(A,B)=A*B$ defined in Eq. 4 yields a Thompson metric midpoint of $A,B\in\mathbb{P}_{n}$ that is continuous, symmetric and affine-invariant. Moreover, if $A\preceq B$, then $A\preceq\mu(A,B)\preceq B$, and the extremal eigenvalues of $\mu(I,\Sigma)$ depend monotonically on $\Sigma\in\mathbb{P}_{n}$. *

We also note that $A*B$ satisfies a key scaling property which suggests that it may be a plausible candidate for a computationally scalable substitute for the standard geometric mean $A\#B$ of two positive definite matrices.

Proposition 2.8**.**

For any real scalars $a,b>0$ and matrices $A,B\in\mathbb{P}_{n}$, we have

[TABLE]

Proof 2.9**.**

*The result follows upon substituting $\lambda_{i}\left((bB)(aA)^{-1}\right)=\frac{b}{a}\lambda_{i}(BA^{-1})$ into the formula Eq. 4. *

Remark 2.10**.**

The scaling in Eq. 10 of course does not generally hold for a mean of two matrices. Indeed, it does not generally hold for means arising as $d_{\infty}$-midpoints either. For instance, [34] identifies

[TABLE]

*as another $d_{\infty}$-midpoint of $A$ and $B$. Clearly $A\diamond B$ does not scale geometrically in the sense of Eq. 10. *

2.1 An optimization-based formulation

A norm $\|\cdot\|$ on the space of $n\times n$ complex matrices is said to be unitarily invariant if $\|UXV\|=\|X\|$ for all $n\times n$ matrices $X$ and unitary matrices $U,V$. A norm $\Phi$ on $\mathbb{R}^{n}$ is called a symmetric gauge norm if it is invariant under permutations and sign changes of coordinates. Consider the family of affine-invariant metric distances $d_{\Phi}$ on $\mathbb{P}_{n}$ defined as

[TABLE]

where $\|\cdot\|_{\Phi}$ is any unitarily invariant norm on the space of Hermitian matrices of dimension $n$ defined by $\|X\|_{\Phi}:=\Phi(\lambda_{1}(X),\ldots,\lambda_{n}(X)),$ with $\lambda_{\min}(X)=\lambda_{n}(X)\leq\ldots\leq\lambda_{1}(X)=\lambda_{\max}(X)$ denoting the $n$ real eigenvalues of $X$ and $\Phi$ a symmetric gauge norm on $\mathbb{R}^{n}$ [7]. The norms $\|\cdot\|_{\Phi}$ induced by the $l_{p}$-norms on $\mathbb{R}^{n}$ for $1\leq p\leq\infty$ are called the Schatten $p$-norms. For the choice of $\Phi(x_{1},\ldots,x_{n})=(\sum_{i}x_{i}^{2})^{1/2}$, $d_{2}:=d_{\Phi}$ corresponds to the metric distance generated by the standard affine-invariant Riemannian metric on $\mathbb{P}_{n}$ given by $\langle X,Y\rangle_{\Sigma}=\operatorname{tr}(\Sigma^{-1}X\Sigma^{-1}Y)$ for $\Sigma\in\mathbb{P}_{n}$ and Hermitian matrices $X,Y\in T_{\Sigma}\mathbb{P}_{n}$. The length element $ds$ of this geometry satisfies $ds^{2}=\operatorname{tr}\left(\Sigma^{-1}d\Sigma\right)^{2}$. The unique (up to parametrization) Riemannian geodesic from $A$ to $B$ is given by the curve $\gamma_{\mathcal{G}}:[0,1]\rightarrow\mathbb{P}_{n}$ defined by

[TABLE]

This curve is significant as a minimal geodesic for any of the affine-invariant metrics $d_{\Phi}$ [7]. The midpoint of $\gamma_{\mathcal{G}}$ is the matrix geometric mean $A\#B$ Eq. 5, which is a metric midpoint in the sense that $d_{\Phi}(A,A\#B)=d_{\Phi}(A\#B,B)=\frac{1}{2}d_{\Phi}(A,B)$ for any choice of symmetric gauge $\Phi$. With $\Phi(x_{1},\ldots,x_{n})=\max_{i}|x_{i}|$, $d_{\Phi}=d_{\infty}$ yields the distance function that coincides with the Thompson metric [51] on the cone $\mathbb{P}_{n}$

[TABLE]

Therefore, we see that the geometric mean $A\#B$ is also a geometric midrange of $A$ and $B$.

The invariant Finsler metrics Eq. 12 provide a route to geometrically generalize several measures of aggregation of data to the space of positive definite matrices $\mathbb{P}_{n}$. Specifically, the mean, median, and midrange of a collection of real numbers $a_{1},\dots,a_{N}$ can be defined as

[TABLE]

respectively, where $\boldsymbol{1}=(1,\dots,1)\in\mathbb{R}^{N}$ and $\boldsymbol{a}=(a_{1},\dots,a_{N})$. By analogy, one can extend these notions to geometric averages for a collection of data $Y_{1},\dots,Y_{N}\in\mathbb{P}_{n}$ arising as

[TABLE]

where $\Phi_{n}$ denotes the gauge norm on the space of $n\times n$ Hermitian matrices and $\Phi_{N}$ denotes the corresponding gauge function acting on the $N$ distances $d_{\Phi_{n}}(X,Y_{i})$. If $\Phi$ corresponds to the $l_{2}$ vector norm, Eq. 14 yields the geometric mean $\mathcal{G}_{2}$ of $Y_{1},\dots,Y_{N}$, also known as the Karcher mean [8, 10, 38]:

[TABLE]

If $N=2$, the unique solution $\mathcal{G}_{2}(A,B)$ of Eq. 15 coincides with the geometric mean $A\#B$. One can also define geometric medians of $Y_{1},\dots,Y_{N}$ to be solutions to Eq. 14 for the choice of $\Phi(x_{1},\dots,x_{n})=\sum_{i}|x_{i}|$:

[TABLE]

Note that the $d_{1}$ distance of $X\in\mathbb{P}_{n}$ to the identity $I$ takes the form

[TABLE]

It is interesting to compare Eq. 17 to the function $F(X)=\log\det(X)$, which plays an important role in convex optimization [14]. In particular, we have

[TABLE]

If $X\succeq I$, then $\log X\succeq 0$ and hence $d_{1}(X,I)=\operatorname{tr}(\log X)=\log\det X$.

If $\Phi$ corresponds to the $l_{\infty}$-norm, then Eq. 14 yields the geometric midrange problem Eq. 2. In the $N=2$ case, we have already seen that $A*B$ is a solution to the corresponding midrange optimization problem $\operatorname{min}_{X\succ 0}\max\{d_{\infty}(X,A),d_{\infty}(X,B)\}$. The $N$-point problem is studied in more detail in Section 3.

2.2 Extremal ordering property

Here we describe a characterization of the midrange $A*B$ of $A,B\in\mathbb{P}_{n}$ that does not rely on any additional structures on $\mathbb{P}_{n}$ except for the standard Löwner partial order $\succeq$. It is remarkable that such a characterization that is independent of any metric or differential geometric structure on $\mathbb{P}_{n}$ exists.

Theorem 2.11**.**

Let $A,B\in\mathbb{P}_{n}$. Then,

[TABLE]

Proof 2.12**.**

For any $X=X^{*}$, we have the congruence relation

[TABLE]

This matrix is clearly positive semidefinite if and only if $A\succeq XB^{-1}X$, which is equivalent to

[TABLE]

Using the monotonicity of the square root function $\Sigma\mapsto\Sigma^{1/2}$ on $\mathbb{P}_{n}$ and the affine-invariance of the Löwner order, Eq. 18 holds if and only if

[TABLE]

The expression on the right hand side is of course the geometric mean $A\#B$ and thus we have

[TABLE]

If we restrict $X$ to be in the real span of $A$ and $B$, Eq. 19 becomes

[TABLE]

which is equivalent to $a\Sigma+bI\preceq\Sigma^{1/2}$ for $\Sigma=B^{-1/2}AB^{-1/2}$. By diagonalizing $\Sigma^{1/2}$ we obtain a unitary matrix $V$ such that $\Sigma=V^{*}DV$ and $\Sigma^{1/2}=V^{*}D^{1/2}V$, where $D=\operatorname{diag}(\lambda_{1}(\Sigma),\dots,\lambda_{n}(\Sigma))$. Therefore, we have $aD+bI\preceq D^{1/2}$, which is equivalent to

[TABLE]

for $i=1,\dots,n$. If we require that equality hold in Eq. 20 for $i=1$ and $i=n$, so that

[TABLE]

we find that

[TABLE]

For this choice of $a$ and $b$, we have $aA+bB=A*B$. Moreover, Eq. 20 is satisfied for each $i=1,\dots,n$ since

[TABLE]

*Imposing equality in Eq. 20 for any pair of indices other than $i=1$ and $i=n$ would yield coefficients $a$ and $b$ that result in the violation of some of the other inequalities in Eq. 20. Therefore, our choice of $a$ and $b$ is indeed optimal. *

2.3 Differential geometric viewpoint

The set $\mathbb{P}_{n}$ is a smooth manifold whose tangent space $T_{\Sigma}\mathbb{P}_{n}$ at any point $\Sigma\in\mathbb{P}_{n}$ can be identified with the set of $n\times n$ Hermitian matrices $\mathbb{H}_{n}$. The matrix exponential map $X\mapsto e^{X}$ maps $\mathbb{H}_{n}$ bijectively onto $\mathbb{P}_{n}$. Its differential $de^{X}:\mathbb{H}_{n}\rightarrow\mathbb{H}_{n}$ at $X$ is the linear map given by $de^{X}(Z)=\frac{d}{dt}(e^{X+tZ})|_{t=0}$. The following exponential metric increasing property is established by Bhatia in [7]. See [30] for an earlier version of the theorem.

Theorem 2.13**.**

For any symmetric gauge norm $\Phi$ and Hermitian matrices $X$ and $Z$, we have

[TABLE]

*where $\|\cdot\|_{\Phi}$ denotes the unitarily invariant norm induced by $\Phi$. *

This theorem has several important consequences, which we will briefly review. First note that the distance functions $d_{\Phi}$ defined in Eq. 12 are induced by the affine-invariant Finsler structures on $\mathbb{P}_{n}$ given by

[TABLE]

for $\Sigma\in\mathbb{P}_{n}$ and $d\Sigma\in T_{\Sigma}\mathbb{P}_{n}$. For our purposes, we can think of a Finsler structure on $\mathbb{P}_{n}$ as a smoothly varying norm on the tangent bundle of $\mathbb{P}_{n}$. Such a structure can be used to calculate the length of any smooth curve $\gamma$ in $\mathbb{P}_{n}$. We can express any such curve as the image of a curve $\Gamma$ in $\mathbb{H}_{n}$ under the exponential map. In particular, any smooth curve $\gamma:[0,1]\rightarrow\mathbb{P}_{n}$ from $I$ to $\Sigma$ can be expressed as $\gamma(t)=e^{\Gamma(t)}$, where $\Gamma(0)=0$ and $\Gamma(1)=\log(\Sigma)\in\mathbb{H}_{n}$. The length of this curve with respect to the Finsler structure Eq. 21 is

[TABLE]

The last integral is simply the length of the curve $\Gamma$ in $\mathbb{H}_{n}$ and the least value it can take is $\|\log\Sigma\|_{\Phi}$, which is attained by the straight line segment from [math] to $\log\Sigma$ in $\mathbb{H}_{n}$. The distance between $I$ and $\Sigma$ is defined as $d_{\Phi}(I,\Sigma)=\inf_{\gamma}L_{\Phi}[\gamma]$, where the infimum is taken over all smooth curves $\gamma$ from $I$ to $\Sigma$. Therefore, we see that

[TABLE]

and note that this distance is attained by the curve $\gamma(t)=e^{t\log\Sigma}$, which is a geodesic from $I$ to $\Sigma$. See Fig. 2.

Since congruence transformations are isometries of $(\mathbb{P}_{n},d_{\Phi})$, it follows that the curve $\gamma_{\mathcal{G}}^{A,B}(t)=A^{1/2}\exp(t\log(A^{-1/2}BA^{-1/2}))A^{1/2}$, is a geodesic from $A$ to $B$. Note that this is precisely in agreement with Eq. 13. This geodesic is unique provided that the geodesics in $\mathbb{R}^{n}$ induced by $\Phi$ are unique. In particular, uniqueness of geodesics in $(\mathbb{P}_{n},d_{\Phi})$ is inherited from $\mathbb{R}^{n}$ when $\Phi$ corresponds to the $l_{p}$-norms for $1<p<\infty$, but not for $p=1,\infty$.

The exponential metric increasing property can also be used to show that the metric space $(\mathbb{P}_{n},d_{\Phi})$ is a space of non-positive curvature for any choice of $\Phi$. See [7, 15] for further details. A closely related result is the geodesic convexity [31] of $d_{\Phi}$, which follows from the inequality $d_{\Phi}(e^{tX},e^{tZ})\leq td_{\Phi}(e^{X},e^{Z})$, for all $X,Z\in\mathbb{H}_{n}$ and $0\leq t\leq 1$. More generally, the geodesic convexity theorem states that for all $A_{1},A_{2},B_{1},B_{2}\in\mathbb{P}_{n}$, the real function

[TABLE]

is convex for any symmetric gauge norm $\Phi$ [7].

In geometric midrange statistics we are interested in the distance $d_{\infty}$, which coincides with the Thompson metric on the cone of positive definite Hermitian matrices. It is known that the Thompson metric does not admit unique minimal geodesics. Indeed, a remarkable construction by Nussbaum in [42] describes a family of geodesics that generally consists of an infinite number of curves connecting a pair of points in a cone $C$. In particular, setting $\alpha:=1/M(x/y;C)$ and $\beta:=M(y/x;C)$, the curve $\phi:[0,1]\rightarrow C$ given by

[TABLE]

is always a minimal geodesic from $x$ to $y$ with respect to the Thompson metric. The curve $\phi$ defines a projective straight line in the cone. If we take $C$ to be the cone of positive semidefinite matrices with interior $\operatorname{int}C=\mathbb{P}_{n}$, then for a pair of points $A,B\in\mathbb{P}_{n}$, we have $\beta=M(B/A;C)=\lambda_{\max}(BA^{-1})$ and $\alpha=1/M(A/B;C)=\lambda_{\min}(BA^{-1})$. Thus, the minimal geodesic described by Eq. 23 takes the form

[TABLE]

where $\lambda_{\max}$ and $\lambda_{\min}$ denote the largest and smallest eigenvalues of $BA^{-1}$, respectively. Taking the midpoint $t=1/2$ of this geodesic, we recover the $d_{\infty}$-midpoint $A*B$ Eq. 4 of $A$ and $B$. Thus, we have arrived at another interpretation of $A*B$ as the midpoint of a suitable geodesic in $(\mathbb{P}_{n},d_{\infty})$. The result follows from elementary algebraic simplification upon setting $A*B=\phi(1/2;A,B)$ in the case $\lambda_{\min}\neq\lambda_{\max}$. If $\lambda_{\min}=\lambda_{\max}$, then $\phi(1/2;A,B)=\sqrt{\lambda_{\min}}A$ also agrees with the formula in Eq. 4. It is shown in [34] that $\phi=\phi(t;A,B)$ is the unique $d_{\infty}$ geodesic connecting $A$ to $B$ if and only if the spectrum of $BA^{-1}$ consists of at most two distinct eigenvalues, one of which is the reciprocal of the other. Moreover, it is shown that otherwise there are infinitely many $d_{\infty}$ minimal geodesics from $A$ to $B$, and that the set of $d_{\infty}$-midpoints of $A$ and $B$ is compact and convex in both Riemannian and Euclidean senses [34].

3 The $N$-point geometric midrange problem

Given a collection of $N$ points $Y_{1},\dots,Y_{N}$ in $\mathbb{P}_{n}$, the midrange problem can be formulated as the following optimization problem

[TABLE]

We call a solution $X^{\star}$ to the above problem a midrange of $\{Y_{i}\}$. Note that the cost function $f(X):=\max_{i}d_{\infty}(X,Y_{i})$ is continuous but not smooth. That is, Eq. 25 is a non-smooth continuous optimization problem on a smooth Finsler manifold.

Proposition 3.1**.**

The optimum cost $t^{\star}=\min_{X\succ 0}\max_{i}d_{\infty}(X,Y_{i})$ of Eq. 25 satisfies $l\leq t^{\star}\leq u$, where the lower and upper bounds are given by

[TABLE]

Proof 3.2**.**

Let $X^{\star}$ denote a midrange of $\{Y_{i}\}$ so that $t^{\star}=\max_{i}d_{\infty}(X^{\star},Y_{i})$. By the triangle inequality, we have for any $i,j=1,\ldots,N$,

[TABLE]

*Taking the maximum of the left-hand side over $i,j$, we arrive at $l=\frac{1}{2}\operatorname{diam}_{\infty}(\{Y_{i}\})\leq t^{\star}$. For the upper bound, note that taking $X=Y_{i}$ for each $i$, we obtain a cost $f(Y_{i})=\max_{j}d_{\infty}(Y_{i},Y_{j})$. The minimum value of these $N$ cost evaluations will clearly still yield an upper bound on the optimum cost $t^{\star}$. Thus we have $t^{\star}\leq u=\min_{i}\max_{j}d_{\infty}(Y_{i},Y_{j})$. *

Note that it is possible to have a collection of points $\{Y_{i}\}$ for which either $l$ or $u$ is attained. For instance, $l$ is clearly attained by the $d_{\infty}$-midpoint when $\{Y_{i}\}$ consists of a pair of points. Similarly, $u$ is attained if we have 3 points $Y_{1},Y_{2},Y_{3}$, where $Y_{3}$ happens to be a $d_{\infty}$-midpoint of $Y_{1}$ and $Y_{2}$. In general, the upper bound is attained if the midrange coincides with one of the data points.

It is instructive to consider the $N$-point affine-invariant midrange of vectors in the positive orthant. In the vector case, the midrange problem in $\mathbb{R}^{n}_{+}$ takes the form

[TABLE]

where $\boldsymbol{x}>0$ means that $\boldsymbol{x}=(x^{a})$ satisfies $x^{a}>0$ for $a=1,\ldots,n$ and $\boldsymbol{y_{i}}$ are a collection of $N$ given points in $\mathbb{R}^{n}_{+}$. As in the matrix case, the optimum cost $t^{\star}=\min_{\boldsymbol{x}>0}f(\boldsymbol{x})=\min_{\boldsymbol{x>0}}\max_{i}\|\log\boldsymbol{x}-\log\boldsymbol{y}_{i}\|_{\infty}$ has a lower bound

[TABLE]

Proposition 3.3**.**

*The lower bound Eq. 28 is attained by $\boldsymbol{x}^{\star}=(x^{a})\in\mathbb{R}^{n}_{+}$ defined by $x^{a}=\left(\min_{i}y_{i}^{a}\cdot\max_{i}y_{i}^{a}\right)^{1/2}$. *

Proof 3.4**.**

Note that

[TABLE]

With $\boldsymbol{x}^{\star}$ as defined in Proposition 3.3, we have

[TABLE]

Remark 3.5**.**

*The midrange problem does not generally have a unique solution in the vector case as can be readily seen through simple examples. For instance, the problem in $\mathbb{R}_{+}^{2}$ with $N=2$ and $\boldsymbol{y}_{1}=(a,1)$, $\boldsymbol{y}_{2}=(1/a,1)$ for some $a>1$ has the solution $\boldsymbol{x}=(1,s)$ for any $s$ satisfying $1/a<s<a$. *

3.1 Geometric midranges via convex optimization

The geometric matrix midrange problem Eq. 25 can be written as

[TABLE]

which has the equivalent epigraph formulation

[TABLE]

This can be rewritten as the quasiconvex problem

[TABLE]

While this problem is not convex due to the presence of the $\log$ function, the feasibility condition $e^{-t}Y_{i}\preceq X\preceq e^{t}Y_{i}$ is convex for fixed $t$ and can be solved using standard convex optimization packages such as CVX [25]. Given a $t$ that is greater than or equal to the optimum value $t^{\star}=\min_{X\succ 0}\max_{i}\;d_{\infty}(X,Y_{i})$, we can solve Eq. 29 using the bisection method [14] by successively solving the feasibility problem as we effectively decrease $t$. In the bisection method it is desirable to have a good estimate for the initial $t$ as the successive reductions in $t$ can be quite slow. In particular, if the lower bound $l=\frac{1}{2}\operatorname{diam}_{\infty}(\{Y_{i}\})$ is attained as in the vector case, then we can solve Eq. 29 in one step by taking $t=l$ and solving the feasibility condition once. However, rather remarkably, numerical examples show that unlike the scalar and vector case, the lower bound $l$ is not always attained in the geometric matrix midrange problem.

Proposition 3.6**.**

*The lower bound $l=\frac{1}{2}\operatorname{diam}_{\infty}(\{Y_{i}\})$ is not necessarily attained in Eq. 25. *

Proof 3.7**.**

Consider the $N=3$ geometric midrange problem in $\mathbb{P}_{2}$ for

[TABLE]

The lower bound $l$ is computed to be $\frac{1}{2}\operatorname{diam}_{\infty}(\{Y_{1},Y_{2},Y_{3}\})=0.7880$. On the other hand, solving the quasiconvex optimization problem Eq. 29 via the bisection method yields the midrange

[TABLE]

*with minimum cost $t^{\star}=0.7901>0.7880=l$. Indeed, we have $t^{\star}=d_{\infty}(X^{\star},Y_{1})=d_{\infty}(X^{\star},Y_{2})=d_{\infty}(X^{\star},Y_{3})$. *

Definition 3.8**.**

(Active matrices)* Let $N\geq 2$ and $(X^{\star},t^{\star})$ be a solution to Eq. 29. Then $Y_{j}$ is called active if $d_{\infty}(X^{\star},Y_{j})=t^{\star}$. In particular, at least one of the following must hold:*

[TABLE]

Proposition 3.6 suggests that the $N$-point matrix midrange problem is richer than the vector case in fundamental ways. While the bisection method applied to the quasiconvex problem Eq. 29 offers a solution, it can be quite slow due to the need for multiple bisection steps and the requirement to compute a reasonable upper bound estimate of the optimum cost for initialization. However, it is possible to recast Eq. 29 as a convex optimization problem and thereby obtain a dramatic improvement in efficiency by introducing new variables. Specifically, by introducing $\xi=e^{t}$ and $\tau=e^{-t}$, and adding the extra convex constraint that $1/\xi-\tau\leq 0$, we find that Eq. 29 can be reformulated as the convex optimization problem

[TABLE]

which can generally be solved much more efficiently than Eq. 29 using standard convex optimization techniques and software packages.

Example 3.9**.**

As an example, we use Eq. 30 to compute the geometric midrange of $N=1000$ real symmetric positive definite $2\times 2$ matrices. The data was generated as $Y_{j}=\Sigma+A_{j}^{T}A_{j}$, for $j=1,\cdots,N$, where the $A_{j}$ are $2\times 2$ matrices with normally distributed entries and $\Sigma\succ 0$ is a fixed matrix. The data matrices can be represented as points in a cone in $\mathbb{R}^{3}$ via the the bijection

[TABLE]

as described in [40]. Fig. 3 shows a visualization of the results of the computation in $\mathbb{R}^{3}$ from two perspectives. The surrounding open cone represents the boundary of the set of real symmetric positive definite matrices and the cloud of points in gray are the data points $Y_{i}$. At optimum, we find that there are 4 active points that are highlighted in red, two of which nearly coincide in the figure. The geometric midrange is highlighted in blue and is the center of the Thompson sphere of radius $t^{\star}$ that defines the smallest enclosing Thompson ball of the data. Note how the active points lie on this sphere. It is interesting that the Thompson ball is the intersection of two cones in this representation. For the sake of comparison, the Karcher mean of the data is also included as a solid black point.

Remark 3.10**.**

*The preceding analysis provides an interesting example of a nonconvex optimization problem that admits a reformulation as a convex optimization problem in the Euclidean sense through a nonlinear change of coordinates. *

While the convex formulation Eq. 30 offers a dramatic improvement to the bisection algorithm applied to the quasiconvex formulation of the problem, we expect that yet more efficient solutions to the problem can be found. In particular, conventional SDP-solvers are based on interior point methods with fast convergence, but high cost per iteration [43, 52], which makes them less suitable for matrices of larger size. Alternatively, one may consider so-called proximal splitting methods such as alternating projections, Douglas-Rachford, or the alternating direction method of multipliers (ADMM) [13, 17, 18, 20] applied to Eq. 29, which have cheap cost per iteration. Unfortunately, these methods tend to have poor convergence properties when the optimal solution is an intersection point of the boundaries of two convex sets with a small intersection angle [21, 22]. Indeed, our numerical experiments indicate that the rates of convergence of such methods degrade as $t$ gets close to the true minimum. This can be expected as each $Y_{i}$ in Eq. 29 defines a bounding box for $X$ through the inequality constraint and $X$ cannot be an interior point to all of them. Ideally, an efficient algorithm for solving this problem would principally rely on the computation of dominant generalized eigenpairs as in the $N=2$ case for which very efficient algorithms exist. In the next subsection, we will consider the optimality conditions for the geometric midrange problem in more detail. Before doing so, we note the following special case for which the $N$-point midrange problem reduces to the 2-point problem as in the scalar case.

Proposition 3.11**.**

*If $Y_{1},\dots,Y_{N}$ are such that $Y_{1}\preceq Y_{i}\preceq Y_{N}$ for all $i=1,\dots,N$, then the geometric midrange of $\{Y_{i}\}$ is given by the set of $d_{\infty}$-midpoints of $Y_{1}$ and $Y_{N}$. *

Proof 3.12**.**

The ordering $Y_{1}\preceq Y_{i}\preceq Y_{N}$ means that the intersection of the feasibility constraints $e^{-t}Y_{i}\preceq X\preceq e^{t}Y_{i}$ in the epigraph formulation Eq. 29 is simply

[TABLE]

*Thus, the optimization problem is unchanged following the elimination of all $Y_{i}$ for $i\neq 1,N$. Hence, the problem is equivalent to the midrange problem for $\{Y_{1},Y_{N}\}$ and is solved by any $d_{\infty}$-midpoint of this pair. Furthermore, the lower bound $l=\frac{1}{2}\operatorname{diam}_{\infty}(\{Y_{i}\})=\frac{1}{2}d_{\infty}(Y_{1},Y_{N})$ is trivially attained. *

Remark 3.13**.**

*Note that in the above we do not assume an order relation between $Y_{i}$ and $Y_{j}$ for $i,j\neq 1,N$. The value of this result lies in the insight that it provides in how and why the matrix $N$-point midrange problem diverges from the scalar and vector case. Fundamentally, no order relation need exist between a pair of matrices, whereas in the scalar case such an ordering is always possible, and similarly an unambiguous ordering is possible at the level of coordinates for vectors. *

3.2 Necessary optimality conditions

Finally, we prove a number of results on the optimality conditions of the geometric matrix midrange problem and the connection between the attainment of the lower bound $l=\frac{1}{2}\operatorname{diam}_{\infty}(\{Y_{i}\})$ and the number of active matrices at optimum.

Proposition 3.14**.**

Let $N\geq 2$ and $(X^{\star},t^{\star})$ be a solution to Eq. 29. Then there exist distinct $i^{\star},j^{\star}\in\{1,\dots,N\}$ such that

[TABLE]

Proof 3.15**.**

By the definition of $(X^{\star},t^{\star})$, there exists at least one index $i^{\star}$ or $j^{\star}$ such that $|\log(\lambda_{\max}(Y_{i^{\star}}^{-\frac{1}{2}}X^{\star}Y_{i^{\star}}^{-\frac{1}{2}}))|=t^{\star}$ or $|\log(\lambda_{\min}(Y_{j^{\star}}^{-\frac{1}{2}}X^{\star}Y_{j^{\star}}^{-\frac{1}{2}}))|=t^{\star}$. In particular, for such $i^{\star}$ and $j^{\star}$ it must hold that

[TABLE]

Next we will show that if there exists only one $i^{\star}$ or $j^{\star}$, then $(X^{\star},t^{\star})$ would not be a solution. To this end, assume that no index such as $j^{\star}$ exists, so that

[TABLE]

Then for sufficiently large $0<k<1$ and $\tilde{X}^{\star}:=kX^{\star}$, it holds that

[TABLE]

*which would mean that $\tilde{X}^{\star}$ is a feasible solution of smaller cost than $t^{\star}$. Analogously, it follows that there always exists an index $i^{\star}$ with the required property. *

Proposition 3.16**.**

Recall the $N=2$ geometric midrange problem for $Y_{1}$ and $Y_{2}$ in $\mathbb{P}_{n}$. Set $\alpha:=\lambda_{\max}(Y_{1}^{-\frac{1}{2}}Y_{2}Y_{1}^{-\frac{1}{2}})$ and $\beta:=\lambda_{\min}(Y_{1}^{-\frac{1}{2}}Y_{2}Y_{1}^{-\frac{1}{2}})$. If $\alpha\neq\beta$, then

[TABLE]

is the only midrange of $\{Y_{1},Y_{2}\}$ in $\{k_{1}Y_{1}+k_{2}Y_{2}:k_{1},k_{2}\geq 0\}$ for which the following is satisfied:

[TABLE]

*The optimal cost to Eq. 29 is given by $t^{\star}=\frac{1}{2}\max\{|\log(\alpha)|,|\log(\beta)|\}$. Furthermore, if $Y_{2}V=Y_{1}VD$ is a generalized eigenvalue decomposition such that $V^{*}Y_{1}V=I$ and $D$ is diagonal, then $V^{*}X^{\star}V$ is diagonal. *

Proof 3.17**.**

Using the linear ansatz $X^{\star}=k_{1}Y_{1}+k_{2}Y_{2}$, we obtain:

[TABLE]

Substituting these expressions into Eq. 32, we find that

[TABLE]

Hence

[TABLE]

*which is equivalent to $(k_{2}^{2}\alpha\beta-k_{1}^{2})(\alpha-\beta)=0$ and implies that $k_{1}^{2}=k_{2}^{2}\alpha\beta$ since $\alpha\neq\beta$. Substituting $k_{1}=\sqrt{\alpha\beta}k_{2}$ into the first equation of Eq. 33 gives $k_{2}^{2}\beta+2k_{2}^{2}\sqrt{\alpha\beta}+k_{2}^{2}\alpha=k_{2}^{2}(\sqrt{\alpha}+\sqrt{\beta})^{2}=1$ and thus Eq. 31. That the cost is given by $t^{\star}$ is trivial and optimality of $X^{\star}$ follows by the attainment of the lower bound. Finally, $V^{*}X^{\star}V$ is diagonal by Proposition 3.14. *

Remark 3.18**.**

*Note that if $\alpha=\beta$ in the statement of the previous proposition, then $Y_{1}^{-1/2}Y_{2}Y_{1}^{-1/2}=\alpha I$, which is equivalent to $Y_{2}=\alpha Y_{1}$. The midrange $Y_{1}\star Y_{2}=\sqrt{\alpha}Y_{1}$ can then be obtained as a conic combination of $Y_{1}$ and $Y_{2}=\alpha Y_{1}$ in a non-unique way. *

In the remainder of this section, we explore the significance of the number of active points at optimum for the attainment of the lower bound of Eq. 29 when $N\geq 2$.

Lemma 3.19**.**

Let $N\geq 2$ and $(X^{\star},t^{\star})$ be a solution to Eq. 29. Then the following are equivalent:

$Y_{j}$* is active with $\log\left(\lambda_{\max}\left(Y_{j}^{-\frac{1}{2}}{X}^{\star}Y_{j}^{-\frac{1}{2}}\right)\right)=t^{\star}$* 2. 2.

$\lambda_{\min}(e^{t^{\star}}Y_{j}-X^{\star})=0$** 3. 3.

$\nexists\varepsilon>0:X^{\star}\preceq e^{t^{\star}}Y_{j}-\varepsilon I$**

Analogously, we have the equivalences:

$Y_{j}$* is active with $-\log\left(\lambda_{\min}\left(Y_{j}^{-\frac{1}{2}}{X}^{\star}Y_{j}^{-\frac{1}{2}}\right)\right)=t^{\star}$* 2. 2.

$\lambda_{\max}(X^{\star}-e^{-t^{\star}}Y_{j})=0$** 3. 3.

$\nexists\varepsilon>0:X^{\star}\succeq e^{t^{\star}}Y_{j}+\varepsilon I$**

Lemma 3.20**.**

Let $D=\operatorname{diag}(d_{1},\dots,d_{n})$ with $d_{1}\leq\dots\leq d_{n}$ and $D\preceq X\preceq d_{n}I$. Then,

[TABLE]

*and thus $d_{n}=\lambda_{\max}(X)\geq\lambda_{\max}(X_{11})$. *

Proof 3.21**.**

*From the inequality it follows that $X_{nn}=d_{n}$ and $\lambda_{\max}(X)\leq d_{n}$. Thus, by Courant-Fischer, $e_{n}=(0,0,\dots,1)$ is an eigenvector of $X$ with eigenvalue $d_{n}$ and thus $X$ has the required form. *

Proposition 3.22**.**

Let $(X^{\star},t^{\star})$ be a solution to Eq. 29 and

[TABLE]

Then, $Y_{1}$ and $Y_{2}$ are active with

[TABLE]

Further, if $Y_{2}V=Y_{1}VD$ is a generalized eigenvalue decomposition such that $V^{*}Y_{1}V=I$ and $D=\operatorname{diag}(d_{1},\dots,d_{n})$ with $d_{1}\leq\dots\leq d_{n}=\lambda_{\max}\left(Y_{1}^{-\frac{1}{2}}Y_{2}Y_{1}^{-\frac{1}{2}}\right)$, then

[TABLE]

Proof 3.23**.**

Let $X_{V}:=V^{*}X^{\star}V$. We first show that $X_{V}$ has the claimed structure. To this end, note that by Eq. 29

[TABLE]

*which implies that $e^{-t^{\star}}D\preceq X_{V}\preceq e^{t^{\star}}I$ with $e^{t^{\star}}=\sqrt{d_{n}}$. Therefore, Lemma 3.20 implies the required structure for $X_{V}$. Then by Lemma 3.19 it follows that $e^{t^{\star}}I$ and $e^{-t^{\star}}D$ are active matrices for $(X_{V},t^{\star})$ and thus $Y_{1}$ and $Y_{2}$ are active matrices for $(X^{\star},t^{\star})$. Then by Lemma 3.20 we can conclude the remaining claim. *

Proposition 3.24**.**

*If there are only two active matrices at an optimum $(X^{\star},t^{\star})$ of Eq. 29, then the lower bound $l=\frac{1}{2}\operatorname{diam}_{\infty}(\{Y_{i}\})$ is attained. *

Proof 3.25**.**

Suppose that $Y_{1}$ and $Y_{2}$ are the only two active matrices at $(X^{\star},t^{\star})$ and assume that the lower bound $l$ is not attained so that $l<t^{\star}$. Denote the geodesic Eq. 13 from $X^{\star}$ to the $d_{\infty}$-midpoint $Y_{1}*Y_{2}$ of $Y_{1}$ and $Y_{2}$ by $\gamma_{\mathcal{G}}(s)=\gamma_{\mathcal{G}}(s,X^{\star},Y_{1}*Y_{2})$. By the geodesic convexity of $(\mathbb{P}_{n},d_{\infty})$, the function $s\mapsto d_{\infty}(Y_{j},\gamma_{\mathcal{G}}(s,X^{\star},Y_{1}*Y_{2}))$ is convex for $j=1,2$. Thus, we have

[TABLE]

*for any $s>0$ and $j=1,2$. As all matrices other than $Y_{1}$ and $Y_{2}$ are inactive at $s=0$, we can achieve a local reduction in the cost function by moving a sufficiently small $s>0$ along the geodesic from $X^{\star}$ to $Y_{1}*Y_{2}$, which would contradict the optimality of $(X^{\star},t^{\star})$. Thus, we have $t^{\star}=l$. *

Remark 3.26**.**

*Note that although Proposition 3.24 implies that the optimum $X^{\star}$ will lie in the $d_{\infty}$-midpoint set of the active pair of matrices, it does not imply that any $d_{\infty}$-midpoint of $Y_{1},Y_{2}$ will be a solution. In particular, $Y_{1}*Y_{2}$ may not be a solution even if the only active matrices at optimum are $Y_{1}$ and $Y_{2}$ since $Y_{1}*Y_{2}$ may fail to satisfy one or more of the constraints in Eq. 29. This is in contrast to the scenario in Proposition 3.11, where any midrange of $Y_{1}$ and $Y_{N}$ will be a solution. *

4 Conclusion

We have introduced a theory of geometric midrange statistics for positive definite Hermitian matrices within an optimization framework. We have also established a number of key results including bounds on the optimization problem as well as necessary conditions for optimality. Furthermore, a solution to the $N$-point problem is offered via convex optimization. Special consideration has been given to the $2$-point midrange problem, which was studied in detail from a number of complementary perspectives. The existence of solutions to the $2$-point problem that can be computed using only extremal generalized eigenvalues has significant implications for computational scalability of matrix midrange statistics. We expect this work to offer a solid foundation for future research in statistics based on Thompson geometry and related topics such as $K$-midranges [16, 47] for matrix-valued data. The development of a fast algorithm for the computation of a midrange of $N$ matrices would be an important step in this direction, with weighted inductive schemes and stochastic algorithms offering a promising angle of attack.

Acknowledgments

We are most grateful to Yurii Nesterov for suggesting the change of variables that facilitated the conversion of the quasiconvex formulation of the $N$-point geometric matrix midrange problem to a convex optimization problem.