3 Random Matrix Theory 3.1 Stochastic Properties of Random Matrix Factors 4 Errata for Publications

3.2 Ordered densities of squared singular values of a Gaussian matrix product

This rather simple result is a humble acknowledgment of the great work in finite-size random matrix theory (RMT) by Prof. Gernot Akemann and team at Uni Bielefeld. For finding the ordered densities, a straightforward recursive formulation, in terms of the $\mathrm{MeijerG}$ function, based on the work of Alberto Zanella at CNR in Italy, is utilized.

Note: Several integration formulas for the MeijerG function are known, e.g., see the $\mathrm{MeijerG}$ function reference. Although the expressions are complex, they can be numerically evaluated quite easily via Mathematica or MATLAB. It amazes me that these finite-size RMT densities are even analytically approachable, although, undoubtedly, the asymptotic RMT theory is “more elegant.”

The theorem is as follows. See below for Mathematica code and numerical simulations.

Theorem 7.

Let $\boldsymbol{X}$ and $\boldsymbol{Y}$ be independent $m\times l$ and $l\times q$ random matrices, respectively, $q\leq m\leq l,$ with identical and independently distributed elements $[\boldsymbol{X}]_{ij},\sim\mathcal{CN}(0,1),i=1,\dots,m,j=1,\dots,l,$ and $[\boldsymbol{Y}]_{ij},\sim\mathcal{CN}(0,1),i=1,\dots,l,j=1,\dots,q,$ respectively. Let $x_{1}\geq\dots\geq x_{q}$ denote the squared singular values of the product $\boldsymbol{X}\boldsymbol{Y}.$ The pdf of the $k$ -th squared singular value of $\boldsymbol{X}\boldsymbol{Y},$ denoted by $q_{k}(x_{k};m,l,q),$ is given by

\displaystyle q_{k}(x_{k};m,q)=K_{q_{k}}h^{[k,(),()]}_{k}(x_{k};m,l,q),

(3.1)

where $K_{q_{k}}$ is a constant ensuring that the integral over the pdf is equal to one, and function $h^{[d,\boldsymbol{n},\boldsymbol{m}]}_{k}(x_{k};m,q)$ is given by the recurrence relation

\displaystyle\sum_{i=1}^{|{\mathcal{I}}^{[d,\boldsymbol{n}]}|}\sum_{j=1}^{|{% \mathcal{I}}^{[d,\boldsymbol{m}]}|}h^{[d-1,\boldsymbol{n}^{\prime},\boldsymbol% {m}^{\prime}]}_{k}(x_{k};m,l,q),

(3.2)

$[l,(),()]$ denotes the initial value of $[d,\boldsymbol{n},\boldsymbol{m}],$ and “()” denotes the empty tuple. Tuples $\boldsymbol{n}$ and $\boldsymbol{m}$ are updated as $\boldsymbol{n}^{\prime}\coloneqq\boldsymbol{n}\cup\{(i,[{\mathcal{I}}^{[d,% \boldsymbol{n}]}]_{i})\}$ and $\boldsymbol{m}^{\prime}\coloneqq\boldsymbol{m}\cup\{(j,[{\mathcal{I}}^{[d,% \boldsymbol{m}]}]_{j})\},$ where $i$ and $j$ are the summation indices in (3.2), and $[{\mathcal{I}}^{[d,\boldsymbol{n}]}]_{i}$ and $[{\mathcal{I}}^{[d,\boldsymbol{m}]}]_{j}$ are the $i$ -th and $j$ -th elements of sets ${\mathcal{I}}^{[d,\boldsymbol{n}]}$ and ${\mathcal{I}}^{[d,\boldsymbol{m}]},$ respectively, defined as ${\mathcal{I}}^{[d,\boldsymbol{n}]}\coloneqq\{1,2,\dots,q\}\setminus\pi_{2}(% \boldsymbol{n})$ and ${\mathcal{I}}^{[d,\boldsymbol{m}]}\coloneqq\{1,2,\dots,q\}\setminus\pi_{2}(% \boldsymbol{m}).$ Next, the termination step is given below

$\displaystyle h^{[1,\boldsymbol{n},\boldsymbol{m}]}_{k}(x_{k};m,q)$	$\displaystyle=\sum_{i=1}^{\|{\mathcal{I}}^{[d,\boldsymbol{n}]}\|}\sum_{j=1}^{\|{% \mathcal{I}}^{[d,\boldsymbol{m}]}\|}s\left(\boldsymbol{n}^{\prime},\boldsymbol{% m}^{\prime}\right)\mathrm{G}^{\tiny\begin{bmatrix}2&0\\ 0&2\end{bmatrix}}\left(\begin{matrix}-\\ v_{2}+n-1,m+n+v_{1}-2\end{matrix}\;\middle\|\;x_{k}\right)$
	$\displaystyle\qquad\times\det{\boldsymbol{\Xi}\left(k,m,q,{\mathcal{I}}^{[d+1,% \boldsymbol{n}^{\prime}]},{\mathcal{I}}^{[d+1,\boldsymbol{m}^{\prime}]}\right)}$
	$\displaystyle\qquad\times\prod_{p=1}^{k-1}\mathrm{G}^{\tiny\begin{bmatrix}3&0% \\ 1&3\end{bmatrix}}\left(\begin{matrix}1\\ 0,v_{2}+[{\mathcal{I}}^{[d,\boldsymbol{n}]}]_{p},[{\mathcal{I}}^{[d,% \boldsymbol{m}]}]_{p}+[{\mathcal{I}}^{[d,\boldsymbol{n}]}]_{p}+v_{1}-1\end{% matrix}\;\middle\|\;x_{k}\right)$	(3.3)

where $v_{1}=l-q,v_{2}=m-q,$ $G$ is the Meijer G function [15], $\boldsymbol{\Xi}\left(k,m,q,{\mathcal{I}}^{[d+1,\boldsymbol{n}^{\prime}]},{% \mathcal{I}}^{[d+1,\boldsymbol{m}^{\prime}]}\right)$ is a $(q-l)\times(q-l)$ matrix with elements

\left[\mathrm{G}^{\tiny\begin{bmatrix}2&1\\ 1&3\end{bmatrix}}\left(\begin{matrix}1\\ v_{2}+n,m+n+v_{1}-1,0\end{matrix}\;\middle|\;x_{k}\right)\right]_{ij},

(3.4)

and $i,j=1,\dots,q-l,$

\displaystyle s\left(\boldsymbol{n}^{\prime},\boldsymbol{m}^{\prime}\right)

\displaystyle=(-1)^{\sum_{i=1}^{|\boldsymbol{n}^{\prime}|}[\pi_{1}(\boldsymbol% {n}^{\prime})]_{i}+[\pi_{1}(\boldsymbol{m}^{\prime})]_{i}},

(3.5)

and $\gamma(\cdot,\cdot)$ denotes the lower incomplete Gamma function.

Second projection notation: Let set $s:=\{(1,2),(3,4),(5,6)\}$ be a set containing three tuples. Then, $\pi_{2}(s)$ is the second projection of $s$ given by $\pi_{2}(s)=\{2,4,6\}.$ Below, a sketch of the proof for completeness.

Proof.

The joint pdf of the squares of the singular values of $\boldsymbol{X}\boldsymbol{Y}$ is given in [1, (18)],[9] as

\displaystyle K_{q}\prod_{1\leq i<j\leq q}(x_{j}-x_{i})\det{\boldsymbol{\Delta% }},

(3.6)

where $\boldsymbol{\Delta}$ is a $q\times q$ matrix with elements

\left[\mathrm{G}^{\tiny\begin{bmatrix}2&0\\ 0&2\end{bmatrix}}\left(\begin{matrix}-\\ v2,v1+j-1\end{matrix}\;\middle|\;x_{i}\right)\right]_{ij},

(3.7)

with $i,j=1,\dots,q.$ The pdf of the $k$ -th eigenvalue can be obtained via the procedure given in [19, Sec. IV-B] utilizing the function

\varphi(n,m,x)=\mathrm{G}^{\tiny\begin{bmatrix}2&0\\ 0&2\end{bmatrix}}\left(\begin{matrix}-\\ v_{2}+n-1,m+n+v_{1}-2\end{matrix}\;\middle|\;x\right),

(3.8)

and the integrals in [1, (A7)] and [15] to integrate over the Meijer G function. Lastly, the multiple summations in [19, Sec. IV-B] can be equivalently reformulated to obtain the recursive definition given in (3.2) and (3.3). ∎

3.2.1 Mathematica Reference Code

Now, a Mathematica package for evaluating the ordered densities is as follows. Here, in xyz[symb_, l, v_2, v_1, q_], $s y m b$ is the desired symbol for the variable, e.g., $x, l$ is the index (corresponding to $k$ in the theorem), $v_{1},v_{2},$ and $q$ are as in the theorem above.

⬇

(* ::Package:: *)

BeginPackage["xyl`"]

(*** Exported Symbol ***)

xyl

Begin["`Private`"]

\[CurlyPhi][v2_, v1_, q_, n_, m_] := MeijerG[{{}, {}}, {{v2+n-1, m + n+v1 - 2}, {}}, s];

Iup[v2_, v1_, q_, n_, m_] := MeijerG[{{}, {1}}, {{0, v2+n, m + n+v1 - 1}, {}}, s];

Idown[v2_, v1_, q_, n_, m_] := MeijerG[{{1}, {}}, {{v2+n, m + n+v1 -1}, {0}}, s];

r[i_, n_, q_] := Complement[Range[q], n][[i]]

sg[n_, m_, q_] := Product[(-1)^(FirstPosition[Complement[Range[q], m[[1 ;; i - 1]]], m[[i]]] + FirstPosition[Complement[Range[q], n[[1 ;; i - 1]]], n[[i]]]), {i, 1, Length[n]}]

Dmat[l_, v2_, v1_, q_, ns_, ms_] := Piecewise[{{1, l == q}}, Det[Table[Idown[v2, v1, q, r[n, ns, q], r[m, ms, q]], {n, 1, q - l}, {m, 1, q - l}]]]

xylpdf[l_, v2_, v1_, q_, d_, ns_, ms_] := Sum[xylpdf[l, v2, v1, q, d - 1, Append[ns, n], Append[ms, m]], {n, Complement[Range[q], ns]}, {m, Complement[Range[q], ms]}]

xylpdf[l_, v2_, v1_, q_] := xylpdf[l, v2, v1, q, l, {}, {}]

xylpdf[l_, v2_, v1_, q_, 1, ns_, ms_] := Product[Iup[v2, v1, q, ns[[i]], ms[[i]]], {i, 1, l - 1}]*

Sum[\[CurlyPhi][v2, v1, q, n, m]*sg[Append[ns, n], Append[ms, m], q]*Dmat[l, v2, v1, q, Append[ns, n], Append[ms, m]], {n, Complement[Range[q], ns]}, {m, Complement[Range[q], ms]}]

xyl[symb_,l_?NumericQ,v2_?NumericQ,v1_?NumericQ,q_?NumericQ] := xylpdf[l, v2, v1, q]/(Product[Gamma[n] Gamma[n+v1] Gamma[n+v2],{n,1,q}] Factorial[l-1]) /. s->symb ;

End []

EndPackage[]

3.2.2 Numerical Results

Lastly, the numerical simulations for $v_{1}=6,v_{2}=2,$ and $q=2.$ We see that the results based on the above expressions and those based on Monte-Carlo simulations agree well.

- ARK

3.2.3 Version History

1.

First published: 9th Jul. 2022 on aravindhk-math.blogspot.com
2.

Edited: 16th Dec. 2023 – converted theorem and proof images to LaTeX