Co-rank 1 projections and the randomised Horn problem

TL;DR

This paper explains the similar eigenvalue distributions in rank 1 projections and perturbations of Hermitian matrices using secular equations, unifying different random matrix scenarios.

Contribution

It provides a common derivation for eigenvalue PDFs in rank 1 projection and perturbation cases, extending to multiplicative settings with unitary matrices.

Findings

Eigenvalue PDFs have similar structures in projection and perturbation cases.

Secular equations and Jacobians unify different random matrix models.

Results apply to complex and real Gaussian vectors, and multiplicative unitary models.

Abstract

Let $\hat{\boldsymbol x}$ be a normalised standard complex Gaussian vector, and project an Hermitian matrix $A$ onto the hyperplane orthogonal to $\hat{\boldsymbol x}$. In a recent paper Faraut [Tunisian J. Math. \textbf{1} (2019), 585--606] has observed that the corresponding eigenvalue PDF has an almost identical structure to the eigenvalue PDF for the rank 1 perturbation $A + b \hat{\boldsymbol x} \hat{\boldsymbol x}^\dagger$, and asks for an explanation. We provide this by way of a common derivation involving the secular equations and associated Jacobians. This applies too in related setting, for example when $\hat{\boldsymbol x}$ is a real Gaussian and $A$ Hermitian, and also in a multiplicative setting $A U B U^\dagger$ where $A, B$ are fixed unitary matrices with $B$ a multiplicative rank 1 deviation from unity, and $U$ is a Haar distributed unitary matrix. Specifically, in each case there is a dual eigenvalue problem giving rise to a PDF of almost identical structure.

Equations

$$B=\Pi A\Pi .$$

$$a_1>{\lambda }_1>a_2>{\lambda }_2>\cdots >{\lambda }_{n-1}>a_n$$

$$\Gamma (n)\frac{{\prod }_{1\le j<k\le n-1}({\lambda }_j-{\lambda }_k)}{{\prod }_{1\le j<k\le n}(a_j-a_k)}.$$

$$C=A+b\hat{\mathbf{x}}{\hat{\mathbf{x}}}^{†}$$

$${\mu }_1>a_1>{\mu }_2>\cdots >{\mu }_n>a_n$$

$$\sum _{i=1}^n{\mu }_i=\sum _{i=1}^n{a}_i+b,$$

$$\Gamma (n)\frac{1}{b^{n-1}}\frac{{\prod }_{1\le j<k\le n}({\mu }_j-{\mu }_k)}{{\prod }_{1\le j<k\le n}(a_j-a_k)}.$$

$$UA{U}^{†}+VB{V}^{†},$$

$$\det (\lambda {\mathbb{I}}_n-B)$$

$$=\det (\lambda {\mathbb{I}}_n-A)\det ({\mathbb{I}}_n+(\lambda {\mathbb{I}}_n-A{)}^{-1}A\hat{\mathbf{x}}{\hat{\mathbf{x}}}^{†})$$

$$\displaystyle=\det(\lambda\mathbb{I}_{n}-A)\Big{(}1+\hat{\boldsymbol{x}}^{\dagger}(\lambda\mathbb{I}_{n}-A)^{-1}A\hat{\boldsymbol{x}}\Big{)}.$$

$$\det ({\mathbb{I}}_p+C_{p\times q}{D}_{q\times p})=\det ({\mathbb{I}}_q+D_{q\times p}{C}_{p\times q})$$

$$0$$

$$=1-{\hat{\mathbf{x}}}^{†}\hat{\mathbf{x}}+\lambda {\hat{\mathbf{x}}}^{†}(\lambda \mathbb{I}-A{)}^{-1}\hat{\mathbf{x}}$$

$$=\lambda {\hat{\mathbf{x}}}^{†}(\lambda \mathbb{I}-A{)}^{-1}\hat{\mathbf{x}}.$$

$$\sum _{p=1}^n\frac{w_p}{\lambda -a_p},w_p:=|x_p{|}^2.$$

$$\frac{\Gamma (s_1+\cdots +s_n)}{\Gamma (s_1)\cdots \Gamma (s_n)}\prod _{j=1}^n{w}_j^{s_j-1}$$

$$\sum _{p=1}^n\frac{w_p}{\lambda -a_p}=\frac{{\prod }_{l=1}^{n-1}(\lambda -{\lambda }_l)}{{\prod }_{l=1}^n(\lambda -a_l)}.$$

$$w_j=\frac{{\prod }_{l=1}^{n-1}(a_j-{\lambda }_l)}{{\prod }_{l=1,l\ne j}^n(a_j-a_l)}.$$

$$\Gamma (n)d{w}_1\cdots d{w}_{n-1},$$

$$\Gamma(n)dw_{1}\cdots dw_{n-1}=\Gamma(n)\Big{(}\prod_{j=1}^{n-1}w_{j}\Big{)}\bigg{|}\det\Big{[}{1\over a_{j}-\lambda_{l}}\Big{]}_{j,l=1}^{n-1}\bigg{|}\,d\lambda_{1}\cdots d\lambda_{n-1}.$$

$$\det\Big{[}{1\over a_{j}-\lambda_{l}}\Big{]}_{j,l=1}^{n-1}={\prod_{1\leq j<k\leq n-1}(a_{j}-a_{k})(\lambda_{j}-\lambda_{k})\over\prod_{j,k=1}^{n-1}(a_{j}-\lambda_{k})}.$$

$$\det (\lambda {\mathbb{I}}_n-C)$$

$$=\det (\lambda {\mathbb{I}}_n-A)\det ({\mathbb{I}}_n-b(\lambda {\mathbb{I}}_n-A{)}^{-1}\hat{\mathbf{x}}{\hat{\mathbf{x}}}^{†})$$

$$\displaystyle=\det(\lambda\mathbb{I}_{n}-A)\Big{(}1-b\hat{\boldsymbol{x}}^{\dagger}(\lambda\mathbb{I}_{n}-A)^{-1}\hat{\boldsymbol{x}}\Big{)},$$

$$1-b\sum _{l=1}^n\frac{w_l}{\lambda -a_l}=0,w_l:=|x_l{|}^2.$$

$$1-b\sum _{l=1}^n\frac{w_l}{\lambda -a_l}=\frac{{\prod }_{l=1}^n(\lambda -{\lambda }_l)}{{\prod }_{l=1}^n(\lambda -a_l)}.$$

$$-b{w}_j=\frac{{\prod }_{l=1}^n(a_j-{\lambda }_l)}{{\prod }_{l=1,l\ne j}(a_j-a_l)}.$$

$$\Gamma(n)dw_{1}\cdots dw_{n-1}=\Gamma(n)\Big{(}\prod_{j=1}^{n-1}w_{j}\Big{)}\bigg{|}\det\Big{[}{1\over a_{j}-\lambda_{l}}-{1\over a_{j}-\lambda_{n}}\Big{]}_{j,l=1}^{n-1}\bigg{|}\,d\lambda_{1}\cdots d\lambda_{n-1}.$$

$$\det\Big{[}{1\over a_{j}-\lambda_{l}}-{1\over a_{j}-\lambda_{n}}\Big{]}_{j,l=1}^{n-1}=\prod_{j=1}^{n-1}{(\lambda_{j}-\lambda_{n})\over(a_{j}-\lambda_{n})}\det\Big{[}{1\over a_{j}-\lambda_{l}}\Big{]}_{j,l=1}^{n-1}$$

$$\frac{1}{b}\frac{2x_1-(a_1+a_2+b)}{a_1-a_2}$$

$$\frac{\Gamma (m_1+\cdots +m_n)}{\Gamma (m_1)\cdots \Gamma (m_n)}\frac{{\prod }_{1\le j<k\le n-1}({\lambda }_j-{\lambda }_k)}{{\prod }_{1\le j<k\le n}(a_j-a_k{)}^{m_j+m_k-1}}\prod _{j=1}^{n-1}\prod _{p=1}^n|{\lambda }_j-a_p{|}^{m_p-1}.$$

$$w_l=\sum _{s=1}^{m_l}|x_l^{(s)}{|}^2,$$

Full text

1 Introduction

Let $A$ be an $n\times n$ complex Hermitian matrix with eigenvalues $a_{1}>a_{2}>\cdots>a_{n}$. Let $\hat{\boldsymbol{x}}$ denote a random $n\times 1$ vector of standard complex Gaussian entries, normalised to have unit length. The matrix $\Pi:=\mathbb{I}_{n}-\hat{\boldsymbol{x}}\hat{\boldsymbol{x}}^{\dagger}$ is then a co-rank 1 projection onto the hyperplane orthogonal to $\hat{\boldsymbol{x}}$. Define

[TABLE]

Interpreting a result of Baryshnikov [1], we know from [11] that the random matrix $B$ has one zero eigenvalue, and non-zero eigenvalues $\{\lambda_{j}\}_{j=1}^{n}$ supported on

[TABLE]

with probability density function (PDF)

[TABLE]

With $A$ and $\hat{\boldsymbol{x}}$ as above, consider next the random matrix

[TABLE]

with $b>0$ a parameter. Interpreting a result of Frumkin and Goldberger [19, Th. 6.1 and 6.7], restated in [10, Th. 5.2], $C$ has eigenvalues, $\{\mu_{i}\}_{i=1}^{n}$ say, supported on

[TABLE]

subject to the constraint

[TABLE]

and with PDF222The normalisation constant $\Gamma(n)$ is given as ${1\over n}$ in [19] and repeated in [10]. This is due to a different convention relating to the implementation of the delta function constraint on the Lebesgue measure in $\mathbb{R}^{n}$ which differs by a factor of $n!$; see the discussion in the second last paragraph of §1 of [10].

[TABLE]

Observing the similarity between (1.2), (1.3) and (1.5), (1.7), the recent paper of Faraut [10] in the concluding Section “Remarks”, asks for an explanation. Here we address this question, showing in Sections 2.1 and 2.2 how to give a unified derivation of both results.

The viewpoint of (1.4) taken in [12, 10] is that of a special case of the random matrix sum

[TABLE]

for $U,V\in U(n)$, chosen with Haar measure. The matrices $A$ and $B$ are fixed Hermitian matrices, and the special case being considered is $B={\rm diag}\,(b,0,\dots,0)$. One remarks that $UAU^{\dagger}$ is the adjoint orbit of the matrix $A$, and similarly the meaning of $VBV^{\dagger}$. Also, since matrices in $U(n)$ diagonalise complex Hermitian matrices, the sum in (1.8) depends only on the eigenvalues of $A$ and $B$. Due to this, the question of the eigenvalue PDF of (1.8) is a randomised version of Horn’s problem [22]. This randomised version appears to have been first studied in [32], in the variant and specialisation of (1.8) for which $A,B$ are real symmetric matrices and $U,V\in{\rm SO}(3)$ (see also Section 3.2 below). Lie algebraic structures including and generalising (1.8) can be found in [7]. Recent years has seen a surge of interest in this problem; see e.g. [35, 10, 4, 33, 34, 3]. The currentness of this activity provides further motivation for extending our study beyond the question posed in [10].

The result (1.3) assumes the eigenvalues of $A$ are all distinct. In the case of the random matrix (1.2) it is known from [11] how to extend (1.3) to the case that $A$ has repeated eigenvalues. We will show in Section 3.1 how to use the methods [11] to calculate the eigenvalue PDF of (1.4) in this setting. In Section 3.2 the case of (1.8) with $A,B$ real symmetric and $U,V\in O(n)$ is considered in the case $B$ having rank 1. In the case $n=3$ this relates to the work [32]. The topic of Section 3.3 is a randomised multiplicative form of Horn’s problem, involving unitary matrices.

When both $A$ and $B$ in (1.8) have full rank, Zuber [35] has recently given a multiple integral formula for the corresponding eigenvalue PDF. This is based on a particular integral over the unitary group due to Harish-Chandra [21], and to Itzykson and Zuber [18], to be referred to as the HCIZ integral. In the case of (1.1) it is known how to use the latter do derive (1.3). In Section 4 we show how to use the result of [35] to reclaim (1.7). Following [33, 10], we also draw attention to the relevance of this integral to the computation of the diagonal entries of the random matrix $U_{p}AU_{p}^{\dagger}$, where $U_{p}$ is the $p\times n$ matrix formed from the first $p$ rows $(p\leq n)$ of $U\in U(n)$, chosen with Haar measure.

We conclude in Section 5 with some remarks relating to the analogue of (1.8) when $A,B$ are real anti-symmetric and $U$ real orthogonal.

2 A unified derivation of (1.3) and (1.7)

2.1 Derivation of (1.3)

We begin by recalling the derivation of the PDF (1.3) due essentially to Baryshnikov [1]; see also [11] and [13, §4.2]. A fundamental point is that the distribution of the random matrix $\hat{\boldsymbol{x}}\hat{\boldsymbol{x}}^{\dagger}$ in the definition of $\Pi$ in (1.1) is unchanged by multiplication on the left or on the right by a unitary matrix. This means that the eigenvalue distribution of $\Pi A\Pi$ is the same as that when $A$ is replaced by the diagonal matrix of its eigenvalues, which we henceforth assume.

Next, with $B$ as specified in (1.1), the fact that $\Pi$ is a projector can be used to check that $B$ and $A\Pi$ have the same eigenvalues. This can be seen from a manipulation of the characteristic polynomial, using (2.2) below with $p=q=n$, $C=-\Pi$, $D=A\Pi$. As a consequence

[TABLE]

To obtain the final equality, the well known formula (see e.g. [13, Exercises 5.2 q.2])

[TABLE]

has been used. The condition for an eigenvalue $\lambda$ of $\Pi A\Pi$ is thus

[TABLE]

We read off from (2.1) that one eigenvalue is always equal to 0, in keeping with $\Pi$ being of co-rank 1, and that the remaining eigenvalues are the zeros of the random rational function

[TABLE]

Here the $x_{p}$ are the components of $\hat{\boldsymbol{x}}$. The latter being a normalised standard complex Gaussian vector tells us that $(|x_{1}|^{2},|x_{2}|^{2},\dots,|x_{n}|^{2})$ is uniformly distributed on the simplex $\sum_{p=1}^{n}|x_{p}|^{2}=1$, or equivalently that this latter vector has Dirichlet distribution, as specified by the PDF

[TABLE]

with $w_{1},\dots,w_{n}>0$ and $\sum_{j=1}^{n}w_{j}=1$, in the case $s_{j}=1$ ($j=1,\dots,n$).

Denote the zeros of (2.4) and thus the non-zero eigenvalues of $\Pi A\Pi$ by $\{\lambda_{j}\}_{j=1}^{n-1}$. Consideration of the graph of (2.4) establishes the interlacing (1.3). Also, we can make use of the zeros of (2.4) to write

[TABLE]

Computing the residue at $\lambda=a_{j}$ gives

[TABLE]

The measure associated with the distribution of $\{w_{p}\}_{p=1}^{n}$ is read off from the special case $s_{p}=1$ of (2.5), and is thus equal to

[TABLE]

subject to the constraints $0<w_{j}<1$ ($j=1,\dots,n-1$) and $\sum_{j=1}^{n-1}w_{j}<1$. We want to change variables to $\{\lambda_{j}\}_{j=1}^{n-1}$. It follows from (2.7) by computing appropriate partial derivatives to form the Jacobian matrix that

[TABLE]

The determinant in (2.9) is referred to as the Cauchy double alternant, and has the well known evaluation (see e.g. [13, Eq. (4.14)])

[TABLE]

Substituting (2.10) in (2.9) gives (1.3).

2.2 Derivation of (1.7)

As for (1.2), the invariance of $\hat{\boldsymbol{x}}\hat{\boldsymbol{x}}^{\dagger}$ under conjugation by unitary matrices implies that the eigenvalue problem for $C$ in (1.4) is the same as that when $A$ therein is replaced by its diagonal matrix of eigenvalues. We assume this form of $A$.

For the characteristic polynomial of $C$ we have

[TABLE]

where to obtain the final line use has been made of (2.2). It follows that the eigenvalues of $C$, $\{\lambda_{j}\}_{j=1}^{n}$ say, are the zeros of the random rational function

[TABLE]

As in (2.4) the variables $\{w_{j}\}_{j=1}^{n-1}$ have distribution (2.5) with parameters $s_{l}=1$ ($l=1,\dots,n$).

Consideration of the graph of the LHS of (2.12) implies, under the assumption $b>0$, that the interlacing condition (1.5) holds with $\{\mu_{j}\}_{j=1}^{n}$ relabelled $\{\lambda_{j}\}_{j=1}^{n}$. Next, analogous to (2.6), by regarding (2.12) as a partial fraction expansion involving $\{\lambda_{l}\}$ we have

[TABLE]

Note that equating the coefficient of $1/\lambda$ on both sides of this expression gives the constraint (1.6), telling us in particular that only $\{\lambda_{l}\}_{l=1}^{n-1}$ are independent.

Computing the residue at $\lambda=a_{j}$ gives

[TABLE]

Using this with $\lambda_{n}$ replaced by $b+a_{n}+\sum_{j=1}^{n-1}(a_{j}-\lambda_{j})$ in keeping with (1.6) allows us to compute the appropriate partial derivatives to form the Jacobian matrix for the change of variables from $\{w_{j}\}_{j=1}^{n-1}$ to $\{\lambda_{j}\}_{j=1}^{n-1}$ and so obtain

[TABLE]

Noting that

[TABLE]

and making use of the Cauchy double alternant determinant evaluation (2.10), then substituting the result in (2.15) gives (1.7).

We can readily integrate over $\{x_{j}\}_{j=1}^{n-1}$ for $n=2$ and $n=3$ and so check that the given normalisation is consistent with our conventions. In the case $n=2$, (1.7) with the substitution (2.15) reads

[TABLE]

while (1.5) and (1.6) together imply $a_{1}+b>x_{1}>{\rm max}\,(a_{1},a_{2}+b)$. There are thus two distinct cases: $0<b<a_{1}-a_{2}$ and $b>a_{1}-a_{2}$. Both integrate to give the value unity, in agreement with the normalisation $\Gamma(n)$. In the case $n=3$ we specialise to the choice $(a_{1},a_{2},a_{3})=(a,0,-a)$. The constraints (1.5) and (1.6) then imply $0<x_{2}<a$ and $a<x_{1}<a+b-x_{2}$. It is efficient to now use computer algebra to integrate (2.17) over these regions, with the value unity resulting, and again confirming the normalisation as stated in (1.7).

3 Generalisations

3.1 Degenerate eigenvalues

Suppose the matrix $A$ in (1.2) and (1.4) is of size $N=\sum_{l=1}^{n}m_{l}$ where $m_{l}$ is the multiplicity of the eigenvalue $\lambda_{l}$. In the previous sections it was assumed that $m_{l}=1$ $(l=1,\dots,n)$. With $\Pi$ defined as in (1.2) but now of size $N\times N$, we know from [11] that the matrix (1.2) has one zero eigenvalue, $m_{l}-1$ eigenvalues equal to $a_{l}$ $(l=1,\dots,n)$ and eigenvalues $\{\lambda_{l}\}_{l=1}^{n-1}$ supported on (1.3) with PDF

[TABLE]

It is of interest to compute the eigenvalue PDF of (1.4) in this setting, and so to extend (1.5)–(1.7). As in the case $m_{l}=1$, a minor modification of the working used in [11] to derive (3.1) suffices. This in turn implies that only a minor modification of the working of Section 2.2 is required.

With $n$ in (2.2) replaced by $N$ this equation again holds true in the setting of degenerate eigenvalues. This means that (2.12) is again valid, but now with

[TABLE]

where $x_{l}^{(s)}$ denotes the components of the vector $\hat{\boldsymbol{x}}$ in the same row as the eigenvalue $\lambda_{l}$ (multiplicity $m_{l}$) in the matrix of eigenvalues. Since $\hat{\boldsymbol{x}}$ is a vector of independent standard complex Gaussian entries normalised to have length unity, the variables (3.2) have Dirichlet distribution (2.5) with $s_{j}=m_{j}$ ($j=1,\dots,n)$.

The working of (2.13)–(2.16) is independent of the precise values of $\{m_{j}\}$ and so again applies. However, relative to the case $m_{j}=1$ ($j=1,\dots,n$) there is now a contribution to the PDF obtained by substituting (2.14) in (2.5), giving a final expression very similar to (3.1).

Theorem 1**.**

Let $A$ be a fixed diagonal matrix with diagonal entries $a_{l}$ $(a_{1}>a_{2}>\cdots>a_{n})$ each repeated $m_{l}$ times, and define $N=\sum_{l=1}^{n}m_{l}$. Let $\hat{\boldsymbol{x}}$ be a $N\times 1$ vector of independent standard complex Gaussians and consider the rank 1 perturbed matrix $C=A+b\hat{\boldsymbol{x}}\hat{\boldsymbol{x}}^{\dagger}$. This matrix has eigenvalues $a_{l}$ with multiplicity $m_{l}-1$, and remaining eigenvalues $\{\lambda_{l}\}_{l=1}^{n}$ say supported on (1.5) and subject to the constraint (1.6) with the eigenvalue PDF

[TABLE]

(cf. (3.1) and its support (1.3)).

3.2 Adjoint orbits involving real orthogonal matrices

Consider the variant of (1.1) in which $A$ is an $n\times n$ real symmetric matrix with eigenvalues $a_{1}>a_{2}>\cdots>a_{n}$ and $\Pi=\mathbb{I}_{n}-\hat{\mathbf{x}}\hat{\mathbf{x}}^{T}$, with $\hat{\mathbf{x}}$ a random $n\times 1$ vector of standard real Gaussian entries, normalised to have unit length. We know from [11, Cor. 1 with $\beta=1$, $m_{i}=1$, $(i=1,\dots,n)$] that the eigenvalue PDF of the $n-1$ non-zero eigenvalues $\{\lambda_{j}\}$ is given by (3.1) with $m_{l}=1/2$, $l=1,\dots,n$, and is thus equal to

[TABLE]

with support given by (1.3).

The analogous variant of (1.8) is to choose the matrices $A,B$ as real symmetric, and $U,V\in O(n)$. Suppose furthermore that $B={\rm diag}\,(b,0,\dots,0)$. An essential point, already used in the derivation of (3.4) as given in [11], is that the joint distribution of the components an $n\times 1$ vector of standard real Gaussian entries is given by the Dirichlet distribution (2.5) with $s_{j}=1/2$ ($j=1,\dots,n)$. Consideration of the working needed to derive (3.3) then implies that the eigenvalue PDF is given by (3.3) specialised to $m_{l}=1/2$, $l=1,\dots,n$, and $N=n/2$. Explicitly, the eigenvalue PDF equals

[TABLE]

supported on (1.5) and subject to the constraint (1.6).

The first meaningful case of (3.5) is when $n=2$. Introducing the variable $s:=\lambda_{1}-\lambda_{2}$, and the constants $s_{\rm max}:=a_{1}-a_{2}+b$, $s_{\rm min}:=a_{2}-a_{1}+b$, a simple calculation gives that (3.4) can then be reduced to the density for $s$,

[TABLE]

This is a special case ($\beta_{2}=0$) of the density given in [35, Eq. (36)] for the setting under consideration but now with $B$ full rank, $B={\rm diag}\,(b,\beta_{2})$.

The case $n=3$ and $b=1$ was first considered in [32]. Noting the parametrisation [32, Eq. (3.1)], it appears that the computed density [32, Eq. (4.2)] agrees with our (3.5), except that the numerator is abscent. This would seem to be a misprint, as the specialisation $a_{2}=a_{3}=0$ given in [32, Eq. (6.3)] contains the denominator as is consistent with (3.5).

3.3 A multiplicative randomised Horn’s problem

Let $U,V\in U(N)$ be chosen with Haar measure, and let $A,B$ be fixed unitary matrices. Asking for the eigenvalues of the product matrix $UAU^{\dagger}VBV^{\dagger}$ is a randomised form of a multiplicative variant of Horn’s problem (for information and references relating to this multiplicative Horn’s problem without randomisation, see [2, Sec. 12]). The facts that unitary matrices are diagonalised by conjugation by other unitary matrices, and that the Haar measure is invariant under multiplication by fixed unitary matrices, tell us that the eigenvalue PDF of $UAU^{\dagger}VBV^{\dagger}$ depends only on the eigenvalues of $A$ and $B$. In the case that $B$ is of the form ${\rm diag}\,(t,1,\dots,1)$ with $|t|=1$, it is possible to adapt workings already in the literature [12] [13, Exercises 4.2 q.3] to deduce the eigenvalue PDF (cf. (3.3) and its support).

Proposition 2**.**

Let $A$ be a fixed $N\times N$ diagonal unitary matrix, with diagonal entries $e^{i\theta_{l}}$ ($0\leq\theta_{1}<\theta_{2}<\cdots<\theta_{n}<2\pi$), each repeated $m_{l}$ times so that $N=\sum_{l=1}^{n}m_{l}$. Let $t=e^{i\phi}$ and set $B={\rm diag}\,(t,1,\dots,1)$. The random unitary product matrix $UAU^{\dagger}VBV^{\dagger}$, where $U,V\in U(N)$ are chosen with Haar measure, has eigenvalues $e^{i\theta_{l}}$ of multiplicity $m_{l}-1$ $(l=1,\dots,n)$. The remaining $n$ eigenvalues, $\{e^{i\psi_{j}}\}_{j=1}^{n}$ say, are supported on

[TABLE]

and subject to the constraint

[TABLE]

They have eigenvalue PDF

[TABLE]

Proof.

The matrix $UAU^{\dagger}VBV^{\dagger}$ has the same eigenvalues as $AWBW^{\dagger}$ where $W=U^{\dagger}V\in U(N)$ chosen with Haar measure. Now $AWBW^{\dagger}=A(\mathbb{I}_{N}+(t-1)\hat{\mathbf{w}}\hat{\mathbf{w}}^{\dagger})$ where $\mathbf{w}$ denotes the first column of $W$. For the characteristic polynomial of the latter, manipulation analogous to that used in the final two equalities of (2.1) gives the factorised form

[TABLE]

where $q_{j}=\sum_{s=1}^{m_{l}}|w_{j}^{(s)}|^{2}$ with $w_{j}^{(s)}$ denoting the components of the vector $\hat{\mathbf{w}}$ in the same rows as the eigenvalue $e^{i\theta_{j}}$ (multiplicity $m_{j}$).

It follows immediately from (3.10) that the eigenvalues of $A$, $e^{i\theta_{l}}$, with multiplicity greater than 1 remain as eigenvalues of the product matrix, but now with multiplicity $m_{l}-1$. It follows too that the remaining eigenvalues are given by the zeros of the second factor. Writing $\lambda=e^{i\psi}$ and recalling $t=e^{i\phi}$ the implied equation can be written

[TABLE]

Consideration of the graph of the right hand side of this equation implies the interlacing (3.7). Also, with $S=AWBW^{\dagger}$ by taking the determinant we must have $\det S=\det A\det B$ which is the constraint (3.8).

Denoting the second factor in (3.10) by $C_{n}(\lambda)$, and setting $\tilde{\lambda}_{j}=e^{i\psi_{j}}$ we observe that it permits the rational function form

[TABLE]

Taking residues allows us to then deduce

[TABLE]

Using the above, we can read off from the working of [12, Lemma 2] that the Jacobian $J$ for the change of variables from $\{q_{j}\}_{j=1,\dots,n-1}\cup\{t\}$ to $\{\tilde{\lambda}_{j}\}_{j=1,\dots,n-1}\cup\{t\}$ is given by

[TABLE]

The probability density for $\{q_{j}\}$ is given by the Dirichlet distribution (2.5) with $w_{j}=q_{j}$, $s_{j}=m_{j}$ ($j=1,\dots,n)$. Substituting (3.11) for $q_{j}$ and using too (3.12), by changing variables in the corresponding probability measure, wedged with $dt$, we read off (3.9).

∎

4 Some applications of the HCIZ integral

4.1 Derivation of (1.7)

Zuber [35] has initiated a study of the eigenvalues of the random matrix sum (1.8) based on a matrix integral due to Harish-Chandra [21], and Itzykson and Zuber [18]. Let the eigenvalues of the Hermitian matrix $X$ be denoted $x:=(x_{1},\dots,x_{n})$ and those of the Hermitian matrix $Y$ be denoted $y=:(y_{1},\dots,y_{N})$. This matrix integral, referred to as the HCIZ integral for short, then reads (see e.g. [13, Proposition 11.6.1])

[TABLE]

where $[U^{\dagger}dU]$ denotes the normalised Haar measure for $U(n)$, and for an array $u=(u_{1},\dots,u_{n})$, $\Delta_{n}(u):=\prod_{1\leq j<k\leq n}(u_{k}-u_{j})$. In this section we will show how the PDF (1.7) can be derived making use of (4.1). Use of the later to derive (1.3) can be found in [30].

Let $X$ and $C$ be $n\times n$ Hermitian matrices $X=[x_{jk}]_{j,k=1}^{n}$, $C=[c_{jk}]_{j,k=1}^{n}$. Let $X$ be random with distribution having PDF $f(X)$, and define the Fourier-Laplace transform

[TABLE]

From this definition it is immediate that for $X$ and $Y$ independent

[TABLE]

Noting that

[TABLE]

and writing $(dX)=\prod_{1\leq j\leq k\leq n}dx_{jk}^{(r)}\prod_{1\leq j<k\leq n}dx_{jk}^{(i)}$ (here the superscripts $(r)$ and $(i)$ denote the real and imaginary parts) (4.2) can be rewritten

[TABLE]

It is assumed that $f$ decays fast enough that this integral converges. Making use of the usual multi-dimensional inverse Fourier transform shows that (4.5) can be inverted to give

[TABLE]

Suppose now that $\hat{f}_{X}(iC)=\hat{f}_{X}(iUCU^{\dagger})$ for all $U\in U(N)$, and thus is a function of the eigenvalues $c=(c_{1},\dots,c_{n})$ only, which is to be denoted by writing $\hat{f}_{X}(iC)=f_{X}(ic)$. Then (4.6) is a function of the eigenvalues of $X$ only and we write $f(X)=f(x)$, where $x=(x_{1},\dots,x_{n})$. In this setting it can be shown, by averaging over $U$ using the HCIZ integral (4.1), that (4.6) reduces to (see e.g. [27, Eq. (1.6)])

[TABLE]

Let $Z$ denote the random matrix sum (1.8). Making use of (4.3) and then the HCIZ integral to evaluate $\hat{f}_{UAU^{\dagger}}(C)$ and $\hat{f}_{VBV^{\dagger}}(C)$ gives that

[TABLE]

Replace $X$ by $Z$ in (4.7) and substituting (4.8) with $C$ replaced by $ic$ gives us the PDF of $Z$. However we seek not the PDF of $Z$ itself but rather the eigenvalue PDF. This can be read off from the former by recalling that associated with the diagonalisation formula $Z=W^{\dagger}{\rm diag}\,(z_{1},\dots,z_{n})W$, where $\{z_{j}\}$ are the eigenvalues and $W$ is the matrix of the corresponding eigenvectors, is the decomposition of measure (see e.g. [13, Eq. (1.27) with $\beta=2$])

[TABLE]

Since the PDF for $Z$ is dependent only on the eigenvalues $z$, we can integrate over $W$ using (see e.g.  [13, Eq. (1.28) with $\beta=2$])

[TABLE]

With $f(z)$ now denoting the eigenvalue PDF of (1.8), we obtain

[TABLE]

This is the result of Zuber [35, Proposition 1], obtained by following essentially the same steps.

Our specific interest is in the case $b_{1},b_{2},\dots,b_{n-1}\to 0$ and $b_{n}=b$. In this limit

[TABLE]

which follows by taking the limits successively; see also [9]. More explicitly, note that when it comes to taking $b_{l}\to 0$, the first $l-1$ rows of the determinant can be subtracted in appropriate multiples from row $l$ to reduce its leading term to the one proportional to $b_{l}^{l-1}$ in its Maclaurin expansion. The denominator at this stage consists of $1/((b_{l}b_{l+1}\cdots b_{n})^{l-1}\Delta_{N-l+1}(\{b_{j}\}_{j=l}^{N}))$, so the limit $b_{l}\to 0$ can now be taken immediately by operating on only row $l$ of the determinant.

Consider the product of the factor in the integrand of (4.9) $1/\Delta_{n}(c)$ times the determinant in (4.10). We see upon making of Laplace expansion of the latter, then evaluating the cofactors as Vandermonde products that this quantity, which is a symmetric function of $\{c_{j}\}_{j=1}^{n}$, that this simplifies to read

[TABLE]

For the product of the other factors in the integrand, we can write

[TABLE]

Multiplying together (4.11) and (4.12) we see, upon minor manipulation, that the integrand of (4.9) in the limiting case of interest reduces down to

[TABLE]

Consider term $p$ in this sum. The dependence on each $c_{l}$, $(l\neq p)$ occurs soley in column $l$, so for all these variables, the integrations can be done column by column. For these we require

[TABLE]

which follows by a residue computation. Hence, after simple manipulation of the determinant, and with the integration of each $c_{p}$ in the summand still remaining, we are left with

[TABLE]

Integrating over $c_{p}$ is now immediate, showing that the above expression reduces to

[TABLE]

Note that the delta function constraint is the requirement (1.6), with $\{\mu_{i}\}$ relabelled $\{z_{i}\}$.

Consider the determinant in (4.14). We can check that with the $z_{j}$’s ordered $z_{1}>z_{2}>\cdots>z_{n}$, if two of the $a_{j}$’s say $a_{q}$ and $a_{q^{\prime}}$ should fall between two consecutive $z_{j}$’s, or outside of $z_{1}$ or $z_{n}$, then rows $q$ and $q^{\prime}$ are identical, so the determinant vanishes. Considering too the requirement of the delta function, we must therefore have the ordering

[TABLE]

which with $\{\mu_{i}\}$ relabelled $\{z_{i}\}$ is (1.5). With this ordering we can check that only the $p=1$ term is non-zero, with the determinant therein equal to

[TABLE]

Hence (4.14), which is the multiple integral in (4.9) in the limiting case of interest, has the evaluation

[TABLE]

supported on (4.15). Substituting this in (4.9), together with the appropriate factors from (4.11), reclaims (1.7).

An outstanding question along these lines is to develop a method based on matrix transforms to similarly reclaim Proposition 2; see the works [25, 26] for recent results on transforms of random product matrices.

4.2 Distribution of the diagonal entries for $U_{p}AU_{p}^{\dagger}$

It is observed in [10, 33], and in fact much earlier in [16, Eqns. (3)–(5)] (see also [8, 15]) that the HCIZ integral (4.1) has the interpretation as the Fourier-Laplace transform of the distribution of the diagonal entries of the random matrix $UBU^{\dagger}$. Choosing $A={\rm diag}\,(a_{1},\dots,a_{p},0,\dots,0)$ corresponds to the Fourier-Laplace transform of the distribution of the diagonal entries of the random matrix $U_{p}BU_{p}^{\dagger}$ where $U_{p}$ is the $p\times n$ matrix formed by the first $p$ rows of $U$. Such distributions first appeared in a more general context in the work of Heckman [17], and are termed Heckman measures.

Let us consider first the case $p=1$. The matrix $U_{1}BU_{1}^{\dagger}$ is then a scalar quantity, corresponding to a particular random quadratic form.

Proposition 3**.**

Let $\mathbf{z}$ be a row vector chosen uniformly at random from the unit sphere in $\mathbb{C}^{n}$, and let $B$ be an Hermitian matrix with eigenvalues $\{b_{i}\}_{i=1}^{n}$, ordered $b_{1}<\cdots<b_{n}$. Let $h_{n}(x;b):=(b-x)^{n-2}{\rm sgn}\,(b-x)$. The PDF for the distribution of the random quadratic form $\mathbf{z}B\mathbf{z}^{\dagger}$ is supported on $(b_{1},b_{n})$ and is given by

[TABLE]

Proof.

Any one row or column of a Haar distributed member of $U(n)$ is uniformly distributed on the complex unit sphere in $\mathbb{C}^{n}$; see e.g. [6] and references therein. Hence with $U_{1}$ defined as in the text above the statement of the proposition, $U_{1}BU_{1}^{\dagger}\mathop{=}\limits^{\rm d}\mathbf{z}B\mathbf{z}^{\dagger}$. Furthermore, with $A={\rm diag}\,(ia,0,\dots,0)$ we see that ${\rm Tr}\,AUBU^{\dagger}=iaU_{1}BU_{1}^{\dagger}$, so in the limit $a_{1},a_{2},\dots,a_{n-1}\to 0$ with $a_{n}=ia$ the LHS of the HCIZ integral (4.1) can be written

[TABLE]

Taking the limit on the RHS gives

[TABLE]

The PDF is obtained by multiplying (4.18) by ${1\over 2\pi}e^{-iax}$ and integrating over $a$. For the latter task, we observe that the only dependence on $a$ in the determinant is in the final row, so we can effectively integrate this row. However, the integrals must then be considered as generalised functions due to the singularity at the origin (alternatively the factor $1/a^{n-1}$ can be replaced by $1/(a+i\delta)^{n-1}$, and the limit $\delta\to 0^{+}$ be taken at the end). Adapting the former viewpoint (this was done is a similar context in the recent work [33]), and thus using the generalised integral

[TABLE]

gives (4.16).

For $x$ outside the interval $(b_{1},b_{n})$ the $h_{n}(x;b)$ in the last row simplify to $h_{n}(x;b)=(b-x)^{n-2}$ (after possibly removing an overall sign from the row). The determinant can then be seen to vanish, so the support is restricted to $(b_{1},b_{n})$ in keeping with the definition of the quadratic form. ∎

The PDF (4.16) is a piecewise polynomial of degree $n-2$ in $x$. Such a simple structure is to be contrasted with the PDF of the random quadratic form $\mathbf{x}B\mathbf{x}^{\dagger}$, where $\mathbf{x}$ is a real random vector sampled uniformly at random from the sphere in $\mathbb{R}^{n}$ [31, 28], which is a far more complicated function of $x$.

From the original work [17] the PDF for the distribution of the diagonal entries of $UBU^{\dagger}$ is known as a particular $\binom{n}{2}$-fold convolution. For small $n$ more explicit calculations are also possible. For example, with $n=3$, taking the inverse transform of the HCIZ integral we find the PDF

[TABLE]

where we have ordered $b_{3}<b_{2}<b_{1}$ and similarly $x_{3}<x_{2}<x_{1}$. Here $\delta(u)$ denotes the Dirac delta function as in (4.14), while $\chi_{A}=1$ if $A$ is true, and zero otherwise.

There is a well studied Gaussian version of the above diagonal entries problem, which in the case of complex entries has attracted attention for its application to wireless communications [20, 29]. Thus let $\Sigma$ be a positive definite $p\times p$ matrix and $G_{p\times n}$ be a standard Gaussian matrix. The $p\times p$ matrix $X=\Sigma^{1/2}GG^{T}\Sigma^{1/2}$ is termed a correlated Wishart matrix. It is straightforward to show that the Fourier-Laplace transform of the distribution of the diagonal entries of $X$ is equal to

[TABLE]

where $\beta=1$ (2) in the case of real (complex) entries, and $A={\rm diag}\,(a_{1},a_{2},\dots,a_{p})$. For general $\Sigma$, $p,n$ there is no known structured formulae for the inverse transform, except in the case $p=2$ when the distribution can be expressed in terms of a Bessel function. The reference [20] gives this formula in the context of a study of the complexities faced in analysing the case $p=3$.

5 Concluding remarks

Unitary matrices diagonalise complex Hermitian matrices, while real orthogonal matrices diagonalise real symmetric matrices. According to the more general Lie algebraic view of [7] it is natural to consider the random matrix sum (1.8) in other circumstances which share an analogous relation between $A$ and $U$, $B$ and $V$. For example, suppose $A$ (and also $B$) is a real anti-symmetric matrix. It is well known (see e.g. [23]) that for $n$ even ($n=2N$ say) there exists an element of $O(2N)$ such that conjugation by this matrix puts $A$ into the block diagonal form

[TABLE]

The same holds true for $n$ odd ($n=2N+1$ say) with the block diagonal form now reading

[TABLE]

An integral formula for the eigenvalue PDF of (1.8), which makes use of Harish-Chandra’s [21] extension of (4.1) to such settings, has been given in [35], and the $N=2$ case has been made explicit.

It is known from the works [5, 14, 24] how to compute an explicit eigenvalue PDF for the randomised sum $A+G^{T}BG$, where the pair $(A,B)$ are of the form (5.1), or (5.2), with the further requirement that $N$ is rank 2 (the smallest rank compatible with the structures), and $G$ is a standard real Gaussian matrix. Replacing $G$ by a Haar distributed real orthogonal matrix, we have found that specialisation of the integral formulas in [35] to this low rank perturbation setting does not lead to a simple structured formula analogous to (1.7) or (3.5). Rather the fact that the perturbation is rank 2 leads to much more complicated structures involving a vast number of terms, conveying little information as to the salient analytic features (the support, singularities near the boundary etc.).

Acknowledgements

This work is part of a research program supported by the Australian Research Council (ARC) through the ARC Centre of Excellence for Mathematical and Statistical frontiers (ACEMS). PJF also acknowledges partial support from ARC grant DP170102028, and JZ acknowledges the support of a Melbourne postgraduate award, and an ACEMS top up scholarship. We thank J. Faraut for helpful remarks.