Riemannian statistics meets random matrix theory: towards learning from high-dimensional covariance matrices

I Introduction

The present paper aims to bring together the tools of Riemannian statistics and those of random matrix theory, in order to unlock the computational issues which have stood in the way of applying Riemannian learning models to datasets of high-dimensional covariance matrices.

Over the past few years, Riemannian Gaussian distributions, and mixtures of these distributions, were introduced as a means of modeling the intrinsic structure of statistical populations of covariance matrices [Cheng2013, Sa16, Sa17]

. They have found successful applications in fields such as brain-computer interface analysis and artificial intelligence

[zanini2018][mathieu2019]. However, such applications could not be pursued for covariance matrices of relatively larger dimension (for example,

50 \times 50

). Indeed, there seemed to exist no practical method of computing the normalising factors associated with Riemannian Gaussian distributions on spaces of high-dimensional covariance matrices.

The theoretical contribution of this paper is to show that this missing method arises quite naturally, as soon as one realises that a strong connection exists between Riemannian Gaussian distributions and random matrix theory. Roughly speaking, Riemannian Gaussian distributions on spaces of real, complex, or quaternion covariance matrices correspond to log-normal orthogonal, unitary or symplectic matrix ensembles. These are similar to the classical, widely-known Gaussian orthogonal, unitary and symplectic ensembles, but with the normal weight function replaced with a log-normal weight function. Thanks to this new connection, the powerful tools of random matrix theory can be employed to uncover several original properties of Riemannian Gaussian distributions, especially for higher-dimensional covariance matrices.

The present paper also extends the equivalence between Riemannian Gaussian distributions and random matrix ensembles to the case of block-Toeplitz matrices. Instead of a log-normal weight function $exp (- {log}^{2} (x))$ , the corresponding weigh function is “ $a c o s h$ -normal" $exp (- {a c o s h}^{2} (x))$ , where $a c o s h$ is the inverse hyperbolic cosine. No attempt is made to develop these additional “ $a c o s h$ -normal" ensembles, at least not for now. In fact, there are already plenty of results to derive and apply with the log-normal ensembles, which cover the important cases of real and complex covariance matrices. This is in part because log-normal ensembles (also called Sitltjes-Wigert ensembles) are amenable to a direct analytic treatment, as will be seen below.

It is also because of a rather surprising connection also noted in [SH21]. The log-normal unitary ensemble had appeared in the theoretical physics literature about twenty years ago, as a random matrix model for the Chern-Simons quantum field theory [Ma05]. A seen in the present paper, this ensemble corresponds to Riemannian Gaussian distributions on the space of complex covariance matrices. Of course, such Riemannian Gaussian distributions have nothing to do with quantum field theory, but some of their valuable properties can be obtained by carefully re-adapting already existing results, found within this theory.

The present paper was developed independently from a theoretical physics paper, published only very recently [Ti20]. Both papers focus on the connection between Riemannian Gaussian distributions and log-normal matrix ensembles, but significantly differ in terms of their contribution and focus, as will be discussed below. The reader is also referred to [Forrester1989] for previous results on asymptotic computations of log-partition functions from the random matrix theory literature, as well as the recent paper [Forrester2021], which considers global and local scaling limits for the $β = 2$ Stieltjes-Wigert random matrix ensemble and associated physical interpretations.

The paper is organised as follows. Section II recalls basic background material on Riemannian Gaussian distributions. Section III presents the main original results, obtained by applying random matrix theory to the study of Riemannian Gaussian distributions. Section IV illustrates the importance of these original results to learning from datesets of high-dimensional covariance matrices. Finally, Section V considers Riemannian Gaussian distributions on the space of block-Toeplitz covariance matrices, and introduces the corresponding “ $a c o s h$ -normal" ensembles. Proofs of the propositions stated in Section III and V are provided in Appendix A.

Ii Notation and background

Let $P_{N}^{β}$ denote the space of $N \times N$ covariance matrices which are either real ( $β = 1$ ), complex ( $β = 2$ ) or quaternion ( $β = 4$ ). Precisely, the elements of $P_{N}^{1}$ are real, symmetric positive-definite matrices, while those of $P_{N}^{2}$ and $P_{N}^{4}$ are respectively complex and quaternion, hermitian positive definite matrices.

Of course, $P_{N}^{β}$

is an open convex cone, sitting inside the real vector space

S_{N}^{β}

of real, complex or quaternion self-adjoint matrices (the adjoint of a matrix being its conjugate transpose). Therefore,

P_{N}^{β}

is a real differentiable manifold, whose tangent space

T_{Y} P_{N}^{β}

at any point

Y \in P_{N}^{β}

is isomorphic to

S_{N}^{β}

. In particular, the dimension of

P_{N}^{β}

{dim}_{N}^{β} = N N_{β}

where

N_{β} = \frac{β}{2} (N - 1) + 1

The elements $Y$ of $P_{N}^{β}$ are in one-to-one correspondence with the centred (zero-mean) $N$

-variate normal distributions (real, complex circular, or quaternion circular, according to the value of

β

). Thus,

P_{N}^{β}

can be equipped with the Rao-Fisher information metric of the centred

N

-variate normal model [amaribook]. This is identical to the so-called affine-invariant metric introduced in [Pennec2006],

⟨ u, v ⟩_{Y} = R t r ((Y^{- 1} u) (Y^{- 1} v)) u, v \in T_{Y} P_{N}^{β} ≃ S_{N}^{β}

(1)

where $R$ denotes the real part and $t r$ denotes the trace. The Riemannian geometry of the metric (1) is quite well-known in the geometric information science community (see the recent book [Pennec2019]). Recall here the associated geodesic distance

(2)

where matrix logarithms and powers are understood as self-adjoint matrix functions, obtained by taking logarithms and powers of eigenvalues. The main advantage of the metric (

1) is its invariance under the action of the group

G_{N}^{β}

N \times N

invertible matrices with real, complex or quaternion entries (according to the value of

β

). For example, if

A \in G_{N}^{β}

and

X

is replaced by

A \cdot X = A X A^{†}

while

Y

is replaced by

A \cdot Y = A Y A^{†}

, then the distance

d (X, Y)

remains unchanged (note that

^{†}

denotes the adjoint, or conjugate transpose).

In terms of the distance (2), we define Riemannian Gaussian distributions on the space $P_{N}^{β}$

to be given by parameterised probability density function

p (Y | ¯ Y, σ)

, as in [Sa16][Sa17],

(3)

where $¯ Y \in P_{N}^{β}$ is the centre-of-mass parameter, and $σ > 0$ the dispersion parameter. The normalising factor $Z (σ)$ is given by the integral

(4)

with respect to the Riemannian volume $d v (Y) = det (Y)^{- N_{β}} {d Y}$ . Here,

\begin{matrix} {d Y} = & \prod_{i \leq j} d Y_{i j} & (β = 1) {d Y} = & \prod_{i \leq j} d Y_{i j}^{(a)} \prod_{i < j} d Y_{i j}^{(b)} & (β = 2) {d Y} = & \prod_{i \leq j} d Y_{i j}^{(a)} \prod_{i < j} d Y_{i j}^{(b)} \prod_{i < j} d Y_{i j}^{(c)} \prod_{i < j} d Y_{i j}^{(d)} & (β = 4) \end{matrix}

where $Y_{i j} = Y_{i j}^{(a)} + Y_{i j}^{(b)} i$ if $Y_{i j}$ is complex and if $Y_{i j}$ is quaternion (here, $i, j, k$ denote complex or quaternion imaginary units).

Both formulae (3) and (4) are greatly simplified by introducing “polar coordinates". Each $Y \in P_{N}^{β}$ can be diagonalised as $Y = U e^{r} U^{†}$ where $e^{r}$ is a diagonal matrix, with diagonal elements $e^{r_{i}}$ for $(r_{1}, \dots, r_{N}) \in R^{N}$ , and , which means $U U^{†} = I_{N}$ (the $N \times N$ identity matrix). Note that $K_{N}^{1}$ is the orthogonal group $O (N)$ , $K_{N}^{2}$ is the unitary group $U (N)$ , and $K_{N}^{4}$ is the symplectic (quaternion unitary) group $S p (N)$ (for deeper insight, see the recent review [edelman]).

Now, as in [Sa17], let $Y$ follow the Riemannian Gaussian density (3), and $Y^{'} = {¯ Y}^{- \frac{1}{2}} Y {¯ Y}^{- \frac{1}{2}}$ . If $Y^{'} = U e^{r} U^{†}$ , then $U$

is uniformly distributed on

(for a rigorous definition of this “uniform distribution", see [meckes]), while

(r_{1}, \dots, r_{N})

have joint probability density function

(5)

where $\propto$ indicates proportionality and $sinh$ the hyperbolic sine. In addition, the normalising factor $Z (σ)$ in (4) reduces to a certain multiple integral $z_{β} (σ)$ . Specifically,

Z (σ) = Ω_{N} z_{β} (σ)

(6)

where $Ω_{N}$ is a numerical constant, which appears after the uniformly distributed matrix $U$ is integrated out of (4), and where

(7)

One of the main issues addressed in the present paper is the efficient approximation of the multiple integral $z_{β} (σ)$ . Until now (in [Sa16][Sa17]), this was done using a Monte Carlo technique which involved a smoothing method (containing certain arbitrarily fixed parameters) and which failed to produce coherent results when the dimension $N$ increased beyond $N = 20$ .

Before proceeding to present our main results, we should briefly note the existence of several alternative proposals for the extension of Gaussian distributions to Riemannian manifolds, including constructions based on heat flows and diffusion processes [Sommer2015]. Further discussion of these alternative formulations is beyond the scope of this paper and the interested reader is referred to the literature on the topic for further information. See [Pennec2019] (Section 3.4.3 and Chapter 10) and the references therein for a recent and comprehensive account.

Iii Main results

The main results of the present paper stem from the equivalence between Riemannian Gaussian distributions on the spaces of real, complex, or quaternion covariance matrices and log-normal orthogonal, unitary, and symplectic matrix ensembles.

Let $Y$ be a random matrix in $P_{N}^{β}$ which follows the Riemannian Gaussian density $p (Y | ¯ Y, σ)$ of (3). It is always possible to assume that $¯ Y = I_{N}$ , since this holds after replacing $Y$ by $Y^{'} = {¯ Y}^{- \frac{1}{2}} Y {¯ Y}^{- \frac{1}{2}}$

. Then, the probability distribution of

Y

is described by the following proposition.

Proposition 1.

Let $Y$ follow the Riemannian Gaussian density (3) with $¯ Y = I_{N}$ . If $X = e^{N_{β} σ^{2}} Y$ , then the probability distribution of $X$ is given by

P (X \in B) \propto \int_{B} e t r [- \frac{{log}^{2} (X)}{2 σ^{2}}] {d X}

(8)

for any measurable subset $B$ of $P_{N}^{β}$ . Here, $e t r (\cdot) = exp (t r (\cdot))$ and the notation ${d X}$ was introduced after (4).

In plain words, this proposition states that $X$ follows a log-normal matrix ensemble. If $X$ is diagonalised as $X = U x U^{†}$ , then $U$ is uniformly distributed on , and the eigenvalues $(x_{1}, \dots, x_{N})$ , which are all positive, have joint probability density function

p (x | σ) \propto | V (x) |^{β} \prod i ρ (x_{i}, 2 σ^{2})

(9)

where $V (x) = \prod_{i < j} (x_{j} - x_{i})$ is the Vandermonde determinant, and $ρ (x, k) = exp (- {log}^{2} (x) / k)$ is the log-normal weight function.

In essence, Proposition 1 is already contained in [Ti20]. As a consequence of this proposition, the multiple integral $z_{β} (σ)$ of (7) may be expressed as follows

z_{β} (σ) = {(2 π σ^{2})}^{N / 2} \times exp [- N N_{β}^{2} (σ^{2} / 2)] \times \frac{1}{N!} \int_{R_{+}^{N}} | V (x) |^{β} ω (d x)

(10)

where $ω (d x) = \prod ¯ ρ (d x_{i})$ (product over $i = 1, \dots, N$ ) with $¯ ρ (d x_{i}) = {(2 π σ^{2})}^{- 1 / 2} ρ (x_{i}, 2 σ^{2}) d x_{i}$ . In [Ti20], integrals as the one in (10) are expressed using the Andreev or De Bruijn identities, often employed in random matrix theory. These yield somewhat cumbersome formulae, involving determinants of size $N \times N$ . On the other hand, the present paper focuses on an alternative approach, which turns out to be more suitable from a practical point of view. Instead of expressing $z_{β} (σ)$ exactly, by means of complicated formulae, the aim is to use a highly efficient approximation, which involves a single analytic expression. This is indicated by the following proposition.

Proposition 2.

In the limit where $N \to \infty$ and $σ \to 0$ , while the product $t = N σ^{2}$ remains constant,

\frac{1}{N^{2}} log z_{β} (σ) ⟶ \frac{β}{2} Φ (\frac{β}{2} t) Φ (ξ) = \frac{ξ}{6} - \frac{{L i}_{3} (e^{- ξ}) - {L i}_{3} (1)}{ξ^{2}}

(11)

where ${L i}_{3}$ is the trilogarithm function, ${L i}_{3} (η) = \sum_{k = 1}^{\infty} η^{k} / k^{3}$ .

In [Ti20], the limit in (11) is only mentioned in passing, and not stated under the same form. Here, this limit will be given centre stage. Proposition 2 states that (11) is valid in the “double-scaling regime" ( $N \to \infty$ and $σ \to 0$ ), but numerical experiments have shown that $log z_{β} (σ)$ can be replaced by the expression afforded by (11) without notable loss of accuracy, whenever $σ$ is small in comparison with $N$ (this is further illustrated below).

At least informally, this can be justified by appealing to arguments originating in theoretical physics [Ma05]. Considered as a function of $t = N σ^{2}$ , $F (t) = log z_{β} (σ)$ is called the Free energy (log of partition function). This free energy can be expanded in an asymptotic series (see Section 1.3 of [Ma05]),

F (t) \sim \infty \sum g = 0 F_{g} (t) {(\frac{t}{N})}^{2 g - 2}

(12)

which is obtained by summing over Feynman diagrams, or so-called fatgraphs. Each coefficient $F_{g} (t)$ is itself a series $F_{g} (t) \sim \sum_{h} F_{g, h} t^{h}$ , where $F_{g, h}$ counts fatgraphs which are said to have $h$ holes and genus $g$ (this means that one thinks of a fatgraph as a graph with $h$ loops, drawn on a surface of genus $g$ , such as a sphere or torus, etc). Now, accepting (12), it follows that for each fixed value of $t$ ,

\frac{1}{N^{2}} F (t) = F_{0} (t) + O (\frac{1}{N^{2}})

(13)

so that $F_{0} (t)$ approximates the left-hand side up to an error of the order of $1 / N^{2}$ . Finally, recalling that $F (t) = log z_{β} (σ)$ , it is clear that $F_{0} (t)$ is the right-hand side of (11) — this follows by uniqueness of asymptotic expansions.

In addition to the asymptotic form of $log z_{β} (σ)$ , another quantity of interest is the asymptotic empirical distribution of eigenvalues, of a random matrix $Y$ which follows the Riemannian Gaussian density $p (Y | ¯ Y, σ)$ with $¯ Y = I_{N}$ . Let $(y_{1}, \dots, y_{N})$ denote the eigenvalues of $Y$ , and consider their empirical distribution

{^ν}_{β} (I) = E [\frac{| y_{i} \in I |}{N}]

(14)

for any open interval $I \subset R_{+}$ , where $E$ denotes expectation, and $| y_{i} \in I |$ the number of $y_{i}$ which belong to $I$ . In [Ti20], the probability density function of ${^ν}_{2}$ was expressed as a weighted sum of Gaussian distributions, by a direct application of the Christoffel-Darboux formula, as in [deift]. It is possible to do so for any value of $N$ , only because the $β = 2$ case can be studied using a well-known family of orthogonal polynomials, called Stieltjes-Wigert polynomials. No such analytic tool is available when $β = 1$ or $4$ .

To make up for this issue, the following proposition provides an asymptotic expression of the distribution ${^ν}_{β}$ , valid for all values $β = 1, 2, 4$ .

Proposition 3.

In the limit where $N \to \infty$ and $σ \to 0$ , while the product $t = N σ^{2}$ remains constant, the empirical distribution ${^ν}_{β}$ converges weakly to a distribution with probability density function $n (y | β t / 2)$ , where

(15)

on the interval $a (ξ) \leq y \leq b (ξ)$ , where $a (ξ) = c (1 + \sqrt{1 - c})^{- 2}$ and $b (ξ) = c (1 - \sqrt{1 - c})^{- 2}$ , with $c = e^{- ξ}$ .

One hopes that, similar to the situation discussed after Proposition 2, the asymptotic density (15) approximates the finite- $N$ empirical distribution ${^ν}_{β}$ to such a good accuracy that one can replace ${^ν}_{β}$ by this asymptotic density, for many practical purposes. This possibility will not be further investigated in the present paper.

Iv Towards learning applications

Riemannian Gaussian distributions were initially proposed as basic building blocks for learning models which aim to capture the intrinsic structure of statistical populations of covariance matrices. These include the mixture models and hidden Markov models, introduced in

[Sa16][mtnshmm] and further developed in [QT21]. In order to make use of these models in real-world applications, it is indispensable to know how to effectively compute the logarithm of the multiple integral

z_{β} (σ)

of (7), for a dimension

N

which may be in the tens or hundreds.

Knowledge of $log z_{β} (σ)$ is already crucial in the simplest situation, where one tries to fit a single Riemannian Gaussian density $p (Y | ¯ Y, σ)$ (rather than a whole mixture) to data . Indeed, this will require setting $¯ Y = {^Y}_{M}$ and $σ = {^σ}_{M}$ , where ${^Y}_{M}$ and ${^σ}_{M}$

are the maximum-likelihood estimates, given in

[Sa16][Sa17]. Specifically,

(16)

is the Fréchet mean of the data $Y_{m}$ , with respect to the distance (2), and ${^σ}_{M}$ is the solution of the nonlinear equation

(17)

Therefore, it is already impossible to solve a toy problem, with a single Riemannian Gaussian density, without having some kind of hold on $log z_{β} (σ)$ .

Until now (in [Sa16][Sa17]), $z_{β} (σ)$ was approximated using an ad hoc Monte Carlo technique, which failed to produce coherent results for a dimension $N$ just above $N = 20$ . In the present section, the aim will be to show that significantly improved results can be obtained by solving Equation (17) after approximating $log z_{β} (σ)$ with the expression afforded by (11), according to Proposition 2.

First, numerical experiments were conducted to verify the validity of this new approximation. According to (13), the right-hand side of (11) should approximate $log z_{β} (σ) / N^{2}$ up to an error of the order of $1 / N^{2}$ . To see that this is correct, the right-hand side of (11) was compared to certain exact expressions of $log z_{β} (σ)$ . Namely, for the $β = 2$ case, one has the following expression, obtained in [habilitation],

log z_{2} (σ) = \frac{N}{2} log (2 π σ^{2}) + N (N^{2} - 1) (σ^{2} / 6) + N - 1 \sum n = 1 (N - n) log (1 - e^{- n σ^{2}})

(18)

and for the $β = 1$ case, when the dimension $N$ is even, the following expression, based on [Ti20],

log z_{1} (σ) = \frac{N}{2} log (2 π σ^{2}) - N (N + 1)^{2} (σ^{2} / 8) + log P f [M (σ)] % (for N even)

(19)

where $P f$ denotes the Pfaffian, equal to the square root of the determinant, and the matrix $M (σ)$ has entries

M_{i j} (σ) = exp [(i^{2} + j^{2}) (σ^{2} / 2)] e r f ((j - i) (σ / 2)) i, j = 1, \dots, N

(20)

with $e r f$ the error function.

Numerical evaluation of (18) for large values of $N$ or $σ$ (up to $σ = 10$ ) is quite straightforward. Moreover, it immediately shows that (11) and (18) agree very closely when $σ$ is smaller than $N$ , and then gradually diverge away from one another as $σ$ increases. Figure 1 provides graphical illustration for $N = 10$ and $20$ . Still larger values of $N$ yield an even stronger match between (11) and (18).

(a) $β = 2$ and $N = 10$ : plot of $log z_{β} (σ) / N^{2}$

It is not equally straightforward to numerically evaluate (19). Even at $N = 10$ , we begin encountering overflow problems for moderate values of $σ$ when performing the computations in MATLAB. Still, as long as these overflow problems do not appear, it is possible to observe a close agreement between (11) and (19), as in the $β = 2$ case. This is shown in Figure 2 for $N = 6$ and $12$ .

(a) $β = 1$ and $N = 6$ : plot of $log z_{β} (σ) / N^{2}$

Based on the numerical results summarised in Figures 1 and 2, it seems possible to use the right-hand side of (11), multiplied by $N^{2}$ , as a substitute for $log z_{β} (σ)$ . While this is only an approximation, it is a highly efficient one, and has the advantage of being given by a rather simple analytic expression. In the $β = 1$ case, direct numerical evaluation of $log z_{β} (σ)$ is unstable for larger values of $σ$ , and (11) offers a practical way out of this problem.

This can be verified by using symbolic computation (Mathematica), which allows us to extend the evaluation of (19) to larger values of $σ$ for larger $N$ . It is nonetheless associated with the drawback that the resulting curves tend to be artificially non-smooth as in Figure 2(a), due to numerical artifacts. These non-smooth features also exist in the curve depicted in Figure 2(b), but do not appear visible at the given resolution. This behaviour is particularly problematic in the context of Equation (17), due to the presence of the derivative of $log z_{β} (σ)$ . Furthermore, symbolic computation of (19) becomes exceedingly slow for sufficiently large $N$ , and this can only be overcome by relying on approximations such as (11).

(a) $β = 1$ and $N = 50$ : plot of $log z_{β} (σ) / N^{2}$ up to $σ = 2$

The second set of experiments directly addressed Equation (17), in the $β = 1$ case, for a dimension $N$ ranging between $10$ and $25$ . This equation was solved using the Newton method, after its right-hand side was approximated according to (11). Here, the solution obtained in this way will be denoted ${~ σ}_{M}$ . This ${~ σ}_{M}$ is an approximation of the maximum-likelihood estimate ${^σ}_{M}$ (the exact solution of (17)). If this approximation is any good, ${~ σ}_{M}$ should approach the true value of $σ$ for sufficiently large $M$ (the number of data points $Y_{m}$ ). This can already be observed at $M = 10^{3}$ when $N = 10$ , as reported in the following Table I

. Each entry in this table gives the average value and standard deviation of

{~ σ}_{M}

, calculated over

20

independent trials (average

\pm

standard deviation). It is clear that increasing

M

from

10^{3}

10^{4}

only reduces the standard deviation of

{~ σ}_{M}

, without really affecting its average.

true $σ$	1	2	3	4	5	6	7
$M = 10^{3}$	$1.09 \pm 0.00$	$2.01 \pm 0.01$	$3.11 \pm 0.30$	$3.56 \pm 0.60$	$4.58 \pm 0.60$	$5.13 \pm 0.50$	$5.60 \pm 0.75$
$M = 10^{4}$	$1.09 \pm 0.00$	$2.01 \pm 0.00$	$3.11 \pm 0.15$	$3.80 \pm 0.16$	$4.50 \pm 0.23$	$5.37 \pm 0.31$	$5.60 \pm 0.37$

TABLE I: The solution

{~ σ}_{M}

of (17), for

β = 1

and

N = 10

(r.h.s. approximated using (11))

When $N = 10$ , the right-hand side of (17) can still be approximated using the Monte Carlo technique mentioned in [Sa16][Sa17]. The solution obtained with this approximation will be denoted . The following Table II shows that the estimation error from ${~ σ}_{M}$ is quite improved, in comparison with the estimation error from . Moreover, ${[{~ σ}_{M}]}_{M C}$ seems to systematically overestimate the true value of $σ$ .

true $σ$	1	2	3	4	5	6	7
$M = 10^{3}$	$1.16 \pm 0.01$	$2.90 \pm 0.01$	$5.63 \pm 0.93$	$7.38 \pm 1.35$	$8.35 \pm 1.50$	$9.40 \pm 1.11$	$9.95 \pm 1.12$
$M = 10^{4}$	$1.16 \pm 0.01$	$2.90 \pm 0.00$	$5.40 \pm 0.37$	$7.00 \pm 0.60$	$9.00 \pm 0.50$	$9.50 \pm 0.31$	$9.77 \pm 0.35$

TABLE II: The solution

{[{~ σ}_{M}]}_{M C}

of (17), for

β = 1

and

N = 10

(MC approximation of r.h.s.)

When $N$ is above $20$ , Monte Carlo approximation of the right-hand side of (17) is not feasible anymore, and one is left only with the approximation using (11). The solution ${~ σ}_{M}$ obtained from this approximation is shown in Table III, for $N = 20$ and $M = 10^{4}$ . Here (contrary to Table II), ${~ σ}_{M}$ systematically underestimates the true value of $σ$ . In fact, identical behavior was observed for $N$ between $20$ and $25$ , along with similar values of ${~ σ}_{M}$ .

A practical means of overcoming this issue would be to include a penalty term into Equation (17), in order to enforce greater values of its solution. Then, the approximation using (11) can be successfully implemented for larger dimension $N$ , where Monte Carlo approximation becomes useless (the present discussion stopped at $N = 25$ , because the sampling algorithms used to generate the data points $Y_{m}$ could not be taken any further).

true $σ$	1	2	3	4	5	6	7
$M = 10^{4}$	$1.02 \pm 0.00$	$1.51 \pm 0.00$	$1.87 \pm 0.01$	$2.35 \pm 0.02$	$2.82 \pm 0.02$	$3.24 \pm 0.03$	$3.63 \pm 0.02$

TABLE III: The solution

{~ σ}_{M}

of (17), for

β = 1

and

N = 20

(r.h.s. approximated using (11))

In conclusion, the present section has demonstrated the new approximation of $log z_{β} (σ)$ , based on (11), significantly improves the existing Monte Carlo approximation, and also extends it to higher dimension. In upcoming work, the experiments conducted in this section will be generalised to more realistic learning models which may be applied to real-world data (for example, the mixture models considered in [zanini2018]).

V Block-Toeplitz covariance matrices

Similar to the spaces $P_{N}^{β}$ introduced in Section II, consider the spaces $D_{N}^{β}$ , defined as follows.

$D_{N}^{2}$ is the space of $N \times N$ matrices with complex entries, whose operator norm is $< 1$ .
$D_{N}^{1}$ is the space of $N \times N$ symmetric matrices with complex entries, whose operator norm is $< 1$ .

Here, operator norm means the largest singular value. The space

D_{N}^{2}

will be called the Hermitian Siegel domain, and

D_{N}^{1}

the symmetric Siegel domain. The focus will be restricted to these

β = 1

2

cases, since they are closely related to time-series analysis and signal processing (the

β = 4

case is similar, but involves quaternion matrices).

Specifically, assume a wide-sense stationary $N$ -variate time series of length $T$ is described by its autocovariance matrix $Γ$ . Wide-sense stationarity implies $Γ$ has a block-Toeplitz structure, with $T \times T$ blocks of size $N \times N$ . When solving an optimal prediction problem, one may apply a multidimensional Szegö-Levinson algorithm to the autocovariance $Γ$ [jeuris1][yannthesis], and obtain a family of matrices $(Γ_{0}, Ω_{1}, \dots, Ω_{T - 1})$ , where $Γ_{0}$ is the zero-lag autocovariance of the original time series (this is a complex covariance matrix, $Γ_{0} \in P_{N}^{2}$ ), and each $Ω_{t}$ is a so-called matrix reflection coefficient, which belongs to $D_{N}^{2}$ . If the autocovariance $Γ$ has Toeplitz blocks, then each $Ω_{t}$ moreover belongs to $D_{N}^{1}$ (which is a subspace of $D_{N}^{2}$ ).

Riemannian geometry enters this picture in the following way [jeuris1][yannthesis]. Consider an information metric on the space of block-Toeplitz autocovariance matrices $Γ$ , equal to the Hessian of the entropy function $S (Γ) = log det (Γ)$ . In terms of the new coordinates $(Γ_{0}, Ω_{1}, \dots, Ω_{T - 1})$ , this information metric is a direct product of the affine-invariant metric (1) on the first coordinate $Γ_{0}$ , and of a scaled copy of the Siegel domain metric on each of the remaining coordinates $Ω_{t}$ (see formula (3.11) in [jeuris1]). Dropping the subscript $t$ , this Siegel domain metric is given by,

(21)

where the tangent space $T_{Ω} D_{N}^{β}$ is isomorphic to the space of $N \times N$ complex matrices if $β = 2$ , and to the space of $N \times N$ symmetric complex matrices if $β = 1$ (this isomorphism shows that the dimension of $D_{N}^{β}$ equals $2 N N_{β}$ ). The geodesic distance associated with the Siegel metric (21) has the following expression [jeuris1][yannthesis],

(22)

where $R (Ξ, Ω) = (Ξ - Ω) (I_{N} - Ω^{†} Ξ)^{- 1} (Ξ^{†} - Ω^{†}) (I_{N} - Ω Ξ^{†})^{- 1}$ is called the matrix cross-ratio.

Riemannian Gaussian distributions on the space $D_{N}^{β}$ are given by their probability density function $p (Ω | ¯ Ω, σ)$ which is of the same form as in (3), but with the distance $d (Ω, ¯ Ω)$ determined by (22), and the normalising factor

(23)

with respect to the Riemannian volume $d v (Ω) = det (I_{N} - Ω Ω^{†})^{- 2 N_{β}} {d Ω}$ . Here, ${d Ω} = \prod_{i j} R d Ω_{i j} I d Ω_{i j}$ where $R$ and $I$ denote the real and imaginary parts (the product is over $i \leq j$ if $β = 1$ and over all $i, j$ if $β = 2$ ).

As shown in [Sa17], a Riemannian Gaussian distribution on the space of block-Toeplitz covariance matrices $Γ$ is just a product of independent Riemannian Gaussian distributions, one for each coordinate $(Γ_{0}, Ω_{1}, \dots, Ω_{T - 1})$ . For $Γ_{0}$ , this is a Riemannian Gaussian distribution on $P_{N}^{2}$ , already considered in Section II. Thus, to understand Riemannian Gaussian distributions of block-Toeplitz covariance matrices, it only remains to study Riemannian Gaussian distributions on the Siegel domain $D_{N}^{β}$ .

In Section III, it was seen that Riemannian Gaussian distributions on $P_{N}^{β}$ are equivalent to log-normal matrix ensembles. On the other hand, Riemannian Gaussian distributions on $D_{N}^{β}$ are equivalent to “acosh-normal" ensembles, which will be described in Proposition 4 below.

Note first that each matrix $Ω \in D_{N}^{β}$ can be factorised in the following way

Ω = U λ V^{†} (if β = 2) Ω = U λ U^{T} (if β = 1)

(24)

where $U$ and $V$ are unitary, $^{T}$ denotes the transpose, and $λ$ is diagonal, with diagonal elements $λ_{i} = t a n h (r_{i})$ for $(r_{1}, \dots, r_{N}) \in R^{N}$ . The $β = 2$

case follows from the singular value decomposition of

Ω

, and the

β = 1

case follows from the Takagi decomposition of

Ω

[edelman]. In either case,

λ_{i}

is of the form

t a n h (r_{i})

because the singular values of

Ω

are all

< 1

Proposition 4.

Let $Ω$ follow a Riemannian Gaussian distribution on the Siegel domain $D_{N}^{β}$ with $¯ Ω = 0_{N}$ ( $N \times N$ zero matrix). If $Ω$ is factorised as in (24), then $U$ and $V$ are uniformly distributed on the unitary group $U (N)$ . Moreover, if $x_{i} = cosh (2 r_{i})$ then $(x_{1}, \dots, x_{N})$ have joint probability density function

p (x | σ) \propto | V (x) |^{β} \prod i ρ (x_{i}, 8 σ^{2}) (x_{1}, \dots, x_{N}) \in (1, \infty)^{N}

(25)

where $V (x) = \prod_{i < j} (x_{j} - x_{i})$ , and where $ρ (x, k) = exp (- {a c o s h}^{2} (x) / k)$ is the “acosh-normal" weight function.

Note that the assumption that the centre-of-mass parameter $¯ Ω$ is equal to $0_{N}$ does not entail any loss of generality. The transformation [yannthesis]

Ω ⟼ Ψ (Ω) = (I_{N} - ¯ Ω {¯ Ω}^{†})^{- \frac{1}{2}} (Ω - ¯ Ω) (I_{N} - {¯ Ω}^{†} Ω)^{- 1} (I_{N} - {¯ Ω}^{†} ¯ Ω)^{\frac{1}{2}}

(26)

maps $¯ Ω$ to $0_{N}$ , while preserving the Siegel domain metric (21) and the associated distance and Riemannian volume. Thus, if $Ω$ follows a Riemannian Gaussian density $p (Ω | ¯ Ω, σ)$ , it is enough to replace $Ω$ by $Ψ (Ω)$ , which will have a Riemannian Gaussian density with $¯ Ω = 0_{N}$ and with the same $σ$ .

Proposition 4 implies that the transformed singular values $x_{i} = cosh (2 r_{i})$ follow the classical eigenvalue distribution of an orthogonal ( $β = 1$ ) or unitary ( $β = 2$ ) matrix ensemble with “acosh-normal" weight function. Therefore, in particular, the normalising factor $Z (σ)$ of (23) reduces to a multiple integral (compare to (6) and (10))

z_{β} (σ) = \frac{1}{N!} \int_{(1, \infty)^{N}} | V (x) |^{β} ω (d x)

(27)

where $ω (d x) = \prod ρ (d x_{i})$ , with $ρ (d x_{i}) = ρ (x_{i}, 8 σ^{2}) d x_{i}$ for the weight function $ρ (x, k) = exp (- {a c o s h}^{2} (x) / k)$ .

Thus, Riemannian Gaussian distributions on the Siegel domain $D_{N}^{β}$ can be studied by applying the tools of random matrix theory to new “acosh-normal" matrix ensembles. Hopefully, in the near future, this will lead to similar results to the ones obtained with log-normal ensembles in Section III, paving the way to learning from datesets of high-dimensional block-Toeplitz covariance matrices.

Appendix A Proofs of Propositions 1 to 4

A-a Proof of Propositon 1

If $Y$ follows the Riemannian Gaussian density (3) with $¯ Y = I_{N}$ , then

P (Y \in B) = {(Z (σ))}^{- 1} \int_{B} exp [- \frac{d^{2} (Y, I_{N})}{2 σ^{2}}] d v (Y)

(28)

Let $(y_{1}, \dots, y_{N})$ denote the eigenvalues of $Y$ . Using (2) and the fact that $d v (Y) = det (Y)^{- N_{β}} {d Y}$ , (28) becomes

(29)

Recall that $X = e^{N_{β} σ^{2}} Y$ . Accordingly, if $(x_{1}, \dots, x_{N})$ are the eigenvalues of $X$ , an elementary calculation yields

exp [- \frac{{log}^{2} (y_{i})}{2 σ^{2}}] y_{i}^{- N_{β}} = exp [N_{β}^{2} (σ^{2} / 2)] \times exp [- \frac{{log}^{2} (x_{i})}{2 σ^{2}}]

Therefore, (29) can be written

P (Y \in B) = exp [N N_{β}^{2} (σ^{2} / 2)] \times {(Z (σ))}^{- 1} \int_{B} e t r [- \frac{{log}^{2} (X)}{2 σ^{2}}] {d Y}

(30)

To conclude, it is enough to use once more the definition $X = e^{N_{β} σ^{2}} Y$ . This implies

{d Y} = exp [- 2 N N_{β}^{2} (σ^{2} / 2)] {d X}

Thus, using the fact that

P (X \in B) = P (e^{N_{β} σ^{2}} Y \in B)

and changing the variable of integration from $Y$ to $X$ in (30), it follows that

P (X \in B) = exp [- N N_{β}^{2} (σ^{2} / 2)] \times {(Z (σ))}^{- 1} \int_{B} e t r [- \frac{{log}^{2} (X)}{2 σ^{2}}] {d X}

as required in (8).

A-B Proof of Proposition 2

The $β = 2$ case follows from (18), by an elementary calculation, after noting that

\frac{1}{N^{2}} N - 1 \sum n = 1 (N - n) log (1 - e^{- n σ^{2}})

is a Riemann sum for the improper integral

For other values of $β$ , the result will be obtained from the reasoning presented in the proof of Proposition 3, based on the scaling equation (33).

A-C Proof of Proposition 3

The $β = 2$ case has already been proved in [habilitation]. To deal with the general case, write (10) under the form $z_{β} (σ) = c_{β} (σ) I_{β} (σ)$ , where

c_{β} (σ) = \frac{1}{N!} exp [- N N_{β}^{2} (σ^{2} / 2)]

and $I_{β} (σ)$ is the multiple integral

I_{β} (σ) = \int_{R_{+}^{N}} exp [- \frac{1}{2 σ^{2}} N \sum i = 1 {log}^{2} (x_{i}) + β \sum i < j log | x_{j} - x_{i} |] d x

Now, if $μ_{N} = (1 / N) \sum_{i = 1}^{N} δ_{x_{i}}$ is the empirical distribution of $x_{i}$ (where $δ_{x_{i}}$ denotes the Dirac measure at $x_{i}$ ), then

I_{β} (σ) = \int_{R_{+}^{N}} exp [- N^{2} E_{β} (μ_{N}, t)] d x

(31)

where $E_{β} (μ, t)$ is the so-called energy functional

(32)

defined for any probability distribution $μ$ on $R_{+}$ . This energy functional satisfies the following scaling equation

(33)

which is easily obtained after dividing (32) by $β / 2$ . The proof can now be completed by applying to (31) and (33) the arguments in [deift]. First [deift] (Corollary 6.90, Page 155),

\frac{1}{N^{2}} log I_{β} (σ) ⟶ - E_{β} (μ^{*}, t)

where $μ^{*} = μ_{β}^{*} (t)$ is the so-called equilibrium distribution, the unique minimiser of the energy (32) among probability distributions on $R_{+}$ . From the scaling equation (33), it is now clear

\frac{1}{N^{2}} log I_{β} (σ) ⟶ - (β / 2) E_{2} (μ^{*}, (β / 2) t)

(34)

Therefore, by adding the limit of $log c_{β} (σ) / N^{2}$ ,

(35)

However, when $β = 2$ , it is already known that this limit is $Φ (t)$ , where $Φ$ is defined in (11). This provides $E_{2} (μ^{*}, t) = - t / 2 - Φ (t)$ , which can be replaced back into (35), yielding the general case of Proposition 2.

To prove Proposition 3, it is possible to use the scaling equation (33), once more. From [deift] (Section 6.4), if ${^μ}_{N}$ is defined as in (14), but with the $x_{i}$ instead of the $y_{i}$ , then ${^μ}_{N}$ converges weakly to the equilibrium distribution $μ^{*} = μ_{β}^{*} (t)$ . On the other hand, the scaling equation (33) implies $μ^{*} = μ_{2}^{*} (β t / 2)$ , and $μ_{2} (t)$ is already known to be the image of the distribution with density $n (y | t)$ , defined as in (15), under the change of variables $x = e^{\frac{β}{2} t} y$ . This provides $μ^{*}$ for any value of $β$ , and the Proposition then follows by changing the variables back from $x$ to $y$ .

A-D Proof of Proposition 4

The proof relies on the general theory of Gaussian distributions on Riemannian symmetric spaces, as outlined in [Sa17][habilitation]. Here, to keep the proof self-contained, it will be helpful to briefly recall certain aspects of this theory (for a more detailed, in-depth discussion, the reader is referred to [habilitation], Sections 1.9 and 3.3).

A Gaussian distribution on a Riemannian symmetric space $M$ is defined by its probability density function

(36)

with respect to the Riemannian volume $d v (x)$ on $M$ , with $¯ x \in M$ and $σ > 0$ .

It is always assumed $M$ is a Riemannian symmetric space of non-positive curvature, associated to a symmetric pair $(G, K)$ . Precisely, $G$ is a connected Lie group which acts transitively and isometrically on $M$ (this action is denoted $x \mapsto g \cdot x$ , for $g \in G$ and $x \in M$ ), and $K$ is a compact subgroup of $G$ , made up of those elements $k \in G$ which fix a certain point $o \in M$ (that is, $k \in K$ if and only if $k \cdot o = o$ ).

For a concrete understanding of Gaussian distributions on $M$ , it is necessary to take a closer look at the Lie algebras of $G$ and $K$ , denoted $g$ and $k$ , respectively. These are related together by the so-called Cartan decomposition, $g = k + p$ (direct sum), where the subspace $p$ of $g$ can be identified with the tangent space $T_{o} M$ (tangent space to $M$ at at $o$ ). In terms of this decomposition, the main construction needed for the proof can be described as follows.

Without any loss of generality, $G$ and $K$ are taken to be matrix Lie groups, and $g$ and $k$ their matrix Lie algebras. Let $a$ be a maximal abelian subspace of $p$ (that is, all the matrices $a \in a$ commute with each other). Then, each matrix $u \in p$ can be written under the form $u = k a k^{- 1}$ , where $k \in K$ and $a \in a$ (this factorisation is the general template for the fifty three matrix factorisations, outlined in [edelman]). Moreover, any $x \in M$ admits a representation

x = exp (k a k^{- 1}) \cdot o k \in K and a \in a

(37)

where $exp$ denotes the matrix exponential. Incidentally, this representation is not unique, but for almost all $x \in M$ (for all $x \in M$ , except a subset of zero volume) $x$ has exactly $| W |$ couples $(k, a)$ which satisfy (37), where $| W |$ is the number of elements of the Weyl group of the symmetric couple $(G, K)$ .

It is now possible to state the following Lemma 1 [Sa17][habilitation], which will yield the entire proof, by direct application.

Lemma 1.

Let $x$ follow the Gaussian density (36) on $M$ , with $¯ x = o$ . If $x$ is represented as in (37), then $k$ is uniformly distributed on the compact group $K$ . Moreover, $a$ has the following probability density function on $a$ (recall that $a$ is a real vector space, so the density is with respect to the usual Lebesgue measure on $a$ )

p (a | σ) \propto exp [- \frac{∥ a ∥_{o}^{2}}{2 σ^{2}}] \prod ρ \in Δ_{+} {(sinh | ρ (a) |)}^{m_{ρ}}

(38)

where $∥ a ∥_{o}$ is the Riemannian norm of $a$ (since $a \in p$ , it can be identified with a vector in $T_{o} M$ ), and $Δ_{+}$ is a set of positive roots $ρ$ on $a$ , with respective multiplicities $m_{ρ}$ (each $ρ \in Δ_{+}$ is a certain linear function $ρ : a \to R$ ).

Recall that positive roots $ρ : a \to R$ are any set of linear functions on $a$ such that, for any $a \in a$ , the eigenvalues of the linear operator ${a d}_{a} : g \to g$ , given by ${a d}_{a} (ξ) = [a, ξ]$ (this is $a ξ - ξ a$ ), are equal to $\pm ρ (a)$ with respective multiplicities $m_{ρ}$ .

The proof may now begin in earnest. It merely consists in identifying $G$ , $K$ , $g$ , $k$ , $p$ , and $a$ , for each Siegel domain $D_{N}^{β}$ , and then writing down the corresponding version of Lemma 1, which directly yields Proposition 4. Fortunately, all the necessary information can be found in [piatetski][helgason].
The symmetric pair : $D_{N}^{β}$ is associated to the symmetric pair $(G, K)$ , where $G$ and $K$ are groups of $2 N \times 2 N$ complex matrices $g$ , defined in the following way (here, $^{*}$ denotes the complex conjugate).

K_{β = 2} = {k = (\begin{matrix} U V \end{matrix}) : U, V \in U (N)}; K_{β = 1} = {k = (\begin{matrix} U U^{*} \end{matrix}) : U \in U (N)}

where $P$ and $S$ denote the following matrices

These groups $G$ and $K$ act on $D_{N}^{β}$ by matrix fractional transformations

g \cdot Ω = (A Ω + B) (C

Riemannian statistics meets random matrix theory: towards learning from high-dimensional covariance matrices

Riemannian Gaussian distributions, random matrix ensembles and diffusion kernels

Gaussian distributions on Riemannian symmetric spaces in the large N limit

A Canonical Representation of Block Matrices with Applications to Covariance and Correlation Matrices

Intrinsic Analysis of the Sample Fréchet Mean and Sample Mean of Complex Wishart Matrices

Group kernels for Gaussian process metamodels with categorical inputs

Manifold Optimisation Assisted Gaussian Variational Approximation

Graph connection Laplacian and random matrices with random blocks

I Introduction

Ii Notation and background

Iii Main results

Proposition 1.

Proposition 2.

Proposition 3.

Iv Towards learning applications

V Block-Toeplitz covariance matrices

Proposition 4.

Appendix A Proofs of Propositions 1 to 4

A-a Proof of Propositon 1

A-B Proof of Proposition 2

A-C Proof of Proposition 3

A-D Proof of Proposition 4

Lemma 1.

Riemannian statistics meets random matrix theory: towards learning from high-dimensional covariance matrices

Related Research

Riemannian Gaussian distributions, random matrix ensembles and diffusion kernels

Gaussian distributions on Riemannian symmetric spaces in the large N limit

A Canonical Representation of Block Matrices with Applications to Covariance and Correlation Matrices

Intrinsic Analysis of the Sample Fréchet Mean and Sample Mean of Complex Wishart Matrices

Group kernels for Gaussian process metamodels with categorical inputs

Manifold Optimisation Assisted Gaussian Variational Approximation

Graph connection Laplacian and random matrices with random blocks

I Introduction

Ii Notation and background

Iii Main results

Proposition 1.

Proposition 2.

Proposition 3.

Iv Towards learning applications

V Block-Toeplitz covariance matrices

Proposition 4.

Appendix A Proofs of Propositions 1 to 4

A-a Proof of Propositon 1

A-B Proof of Proposition 2

A-C Proof of Proposition 3

A-D Proof of Proposition 4

Lemma 1.