Asymptotics of finite system Lyapunov exponents for some random matrix ensembles

While for a general ensemble of $d\times d$ random matrices the Lyapunov spectrum is difficult to compute, there is a class of probability densities on the space of matrices—referred to as isotropic—for which the computation becomes relatively simple [5]. The defining feature of an isotropic matrix distribution on $\{{{A}_{i}}\}$ is that for a particular vector $\vec{x}$, ${{A}_{i}}\vec{x}/|\vec{x}|$ has distribution independent of $\vec{x}$. Examples are when ${{A}_{i}}={{\Sigma }^{1/2}}G_{i}^{(\beta ,d)}$, where Σ is a fixed $d\times d$ positive definite matrix and $G_{i}^{(\beta ,d)}$ is a $d\times d$ standard Gaussian matrix with real $(\beta =1)$, complex $(\beta =2)$ or real quaternion $(\beta =4)$ entries. In the latter case, these entries are themselves $2\times 2$ matrices with complex entries, and the eigenvalues are calculated from the corresponding complex matrix which has a doubly degenerate spectrum; see e.g. [7, section 1.3.2].

The utility of an isotropic matrix distribution on $\{{{A}_{i}}\}$ can be seen by revising the formalism devised in [22, 24] for the computation of the Lyapunov spectrum. The main formula specifies the partial sum of Lyapunov exponents according to

$$\begin{eqnarray}&&{{\mu }_{1}}+\cdots +{{\mu }_{k}}={\rm sup} \mathop{{\rm lim} }\limits_{N\to \infty }\frac{1}{N}{\rm logVo}{{{\rm l}}_{k}}\{{{y}_{1}}(N),...,{{y}_{k}}(N)\}\qquad (k=1,...,d),\end{eqnarray} \tag{ 2.1 }$$

where ${{y}_{j}}(N){{P}_{N}}{{y}_{j}}(0)$, the sup operation is over all sets of linearly independent vectors $\{{{y}_{1}}(0),...,{{y}_{k}}(0)\}$ and ${\rm Vo}{{{\rm l}}_{k}}$ refers to the volume of the parallelogram generated by the given set of k vectors. To quantify the latter, with

$$\begin{eqnarray*}&&{{B}_{N}}[{{y}_{1}}(N){{y}_{2}}(N)\cdots {{y}_{k}}(N)],\end{eqnarray*}$$

so that B_N is the $d\times k$ matrix with its columns given by the k vectors $\{{{y}_{j}}(N)\}$, we have

$$\begin{eqnarray}&&{\rm Vo}{{{\rm l}}_{k}}\{{{y}_{1}}(N),...,{{y}_{k}}(N)\}={\rm det} {{(B_{N}^{\dagger }{{B}_{N}})}^{1/2}}={\rm det} {{(B_{0}^{\dagger }P_{N}^{\dagger }{{P}_{N}}{{B}_{0}})}^{1/2}}.\end{eqnarray} \tag{ 2.2 }$$

Suppose now that $\{{{y}_{1}}(0),...,{{y}_{k}}(0)\}$ is a set of k linearly independent unit vectors, and let ${{E}_{d\times k}}$ denote the $d\times k$ matrix with 1's in the diagonal positions of the k rows, and 0's elsewhere. The fact that the distribution of $G_{i}^{(\beta ,d)}\vec{x}$ is independent of $\vec{x}$ tells us that

$$\begin{eqnarray}&&B_{0}^{\dagger }P_{N}^{\dagger }{{P}_{N}}{{B}_{0}}\mathop{=}\limits^{{\rm d}}\mathop{\prod }\limits_{j=1}^{N}E_{d\times k}^{ \sf T}G_{j}^{(\beta ,d)\,\dagger }\Sigma G_{j}^{(\beta ,d)}{{E}_{d\times k}}=\mathop{\prod }\limits_{j=1}^{N}G_{j,k}^{(\beta ,d)\dagger }\Sigma G_{j,k}^{(\beta ,d)},\end{eqnarray} \tag{ 2.3 }$$

where $G_{j,k}^{(\beta ,d)}$ denotes the $d\times k$ matrix formed by the first k columns of $G_{j}^{(\beta ,d)}$, and it tells us furthermore that the sup operation is redundant. Consequently

$$\begin{eqnarray}&&{{\mu }_{1}}+\cdots +{{\mu }_{k}}=\mathop{{\rm lim} }\limits_{N\to \infty }\mathop{\sum }\limits_{j=1}^{N}\frac{1}{N}{\rm log} \;{\rm det} {{(G_{j,k}^{(\beta ,d)\,\dagger }\Sigma G_{j,k}^{(\beta ,d)})}^{1/2}},\end{eqnarray} \tag{ 2.4 }$$

while the law of large numbers gives that this can be evaluated as the average over $\{G_{k}^{(\beta ,d)}\}$

$$\begin{eqnarray}&&{{\mu }_{1}}+\cdots +{{\mu }_{k}}=\langle {\rm log} \;{\rm det} {{(G_{\cdot ,k}^{(\beta ,d)\,\dagger }\Sigma G_{\cdot ,k}^{(\beta ,d)})}^{1/2}}\rangle .\end{eqnarray} \tag{ 2.5 }$$

Computation of this average in the case $\Sigma ={{\mathbb{I}}_{d}}$ implies [8, 18, 21]

$$\begin{eqnarray}&&{{\mu }_{i}}=\frac{1}{2}{\rm log} \frac{2}{\beta }+\frac{1}{2}\Psi (\beta (d-i+1)/2).\end{eqnarray} \tag{ 2.6 }$$

Suppose (1.1) is generalized so that each A_i is of size $(d+{{\nu }_{i}})\times (d+{{\nu }_{i-1}})$, ${{\nu }_{0}}=0$. It has recently been noted in [17] that the calculations leading to (2.6) generalize and we have

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\mu }_{i}} & = & \frac{1}{2}{\rm log} \frac{2}{\beta }+\frac{1}{2}\langle \langle \Psi (\beta (d-i+1)/2)\rangle \rangle , \\ \langle \langle \Psi (\beta j/2)\rangle \rangle & & \mathop{{\rm lim} }\limits_{N\to \infty }\frac{1}{N}\mathop{\sum }\limits_{i=1}^{N}\Psi (\beta ({{\nu }_{i}}+j)/2). \\ \end{array}\end{eqnarray} \tag{ 2.7 }$$

Note that (2.7) is a symmetric function of $\{{{\nu }_{i}}\}$. This is in fact a general feature of the joint distribution of the eigenvalues of $P_{N}^{\dagger }{{P}_{N}}$ [15], and so without loss of generality we are free to suppose $0={{\nu }_{0}}\leqslant {{\nu }_{1}}\leqslant {{\nu }_{2}}\leqslant \cdots$. The existence of the limit $\mathop{{\rm lim} }\limits_{N\to \infty }{{\nu }_{N}}$ implies that the limit in (2.7) is well defined. Hence only a finite number of the matrices A_i can be rectangular, with all the rest being square. In [17] it was furthermore argued, using exact calculations based on knowledge of the joint distribution of the spectrum of P_N [1, 2, 9, 15, 16] , that the asymptotic form of the eigenvalues ${{{\rm e}}^{2N{{y}_{i}}}}$ in this setting is given by the form (1.4) with ${{\mu }_{i}}$ given by (2.7) and

$$\begin{eqnarray}&&\sigma _{i}^{2}=\frac{1}{4N}\langle \langle {{\Psi }^{{\prime} }}(\beta (d-i+1)/2)\rangle \rangle \end{eqnarray} \tag{ 2.8 }$$

(note that for $\beta =2$—the case of complex entries—this reproduces the formula for $\sigma _{i}^{2}$ in (1.5)). Our aim in the next section is to show how (2.8) can be deduced from analogous calculations as those leading to (2.6).

Our starting point is the fact that (2.1) is valid for finite N, provided that the limiting operation is removed on the rhs, and with $\{{{\mu }_{i}}\}$ such that the eigenvalues of $P_{N}^{\dagger }{{P}_{N}}$ are equal to $\{{{{\rm e}}^{2N{{\mu }_{i}}}}\}$ (see e.g. the discussion in [3, section 6.1]). In the Gaussian case this simplifies to (2.4) without the limiting operation on the rhs. The argument leading to (2.5) then gives that to leading order in N

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \langle {{({{\mu }_{1}}+\cdots +{{\mu }_{k}})}^{2}}\rangle -{{\langle ({{\mu }_{1}}+\cdots +{{\mu }_{k}})\rangle }^{2}} \\ \quad =\frac{1}{N}(\langle {{({\rm log} \;{\rm det} {{(G_{\cdot ,k}^{(\beta ,d)\,\dagger }\Sigma G_{\cdot ,k}^{(\beta ,d)})}^{1/2}})}^{2}}\rangle -{{\langle {\rm log} \;{\rm det} {{(G_{\cdot ,k}^{(\beta ,d)\,\dagger }\Sigma G_{\cdot ,k}^{(\beta ,d)})}^{1/2}}\rangle }^{2}}). \\ \end{array}\end{eqnarray} \tag{ 2.9 }$$

The lhs can be rewriten

$$\begin{eqnarray*}&&\mathop{\sum }\limits_{j=1}^{k}\sigma _{j}^{2}+2{{\sum }_{1\leqslant j\lt l\leqslant k}}(\langle {{\mu }_{j}}{{\mu }_{l}}\rangle -\langle {{\mu }_{j}}\rangle \langle {{\mu }_{l}}\rangle ).\end{eqnarray*}$$

But for large N, and a non-degenerate Lyapunov spectrum, we expect the covariance $\langle {{\mu }_{j}}{{\mu }_{l}}\rangle -\langle {{\mu }_{j}}\rangle \langle {{\mu }_{l}}\rangle$ to decay exponentially fast. This is in keeping with the joint probability density being of the form (1.4). Ignoring such terms, for large N we can substitute

$$\begin{eqnarray}&&\mathop{\sum }\limits_{j=1}^{k}\sigma _{j}^{2}\end{eqnarray} \tag{ 2.10 }$$

for the lhs of (2.9). It thus remains to compute the rhs. In the case $\Sigma ={{\mathbb{I}}_{k}}$ this is a straightforward task.

Lemma 2.1. We have

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle {{({\rm log} \;{\rm det} {{(G_{\cdot ,k}^{(\beta ,d)\,\dagger }G_{\cdot ,k}^{(\beta ,d)})}^{1/2}})}^{2}}\rangle & = & \frac{1}{4}\left( \mathop{\sum }\limits_{j=1}^{k}{{\Psi }^{{\prime} }}((\beta /2)(d-k)+j\beta /2) \right. \\ {} & {} & +\left. {{(\mathop{\sum }\limits_{j=1}^{k}({\rm log} \frac{2}{\beta }+\Psi ((\beta /2)(d-k)+j\beta /2)))}^{2}} \right) \\ \end{array}\end{eqnarray} \tag{ 2.11 }$$

and thus the rhs of (2.9) is equal to

$$\begin{eqnarray}&&\frac{1}{4N}\mathop{\sum }\limits_{j=1}^{k}{{\Psi }^{{\prime} }}((\beta /2)(d-k)+j\beta /2).\end{eqnarray} \tag{ 2.12 }$$

Proof. We follow the strategy of the proof given in [8, section 2.2] of the evaluation of the rhs of (2.5) in the case $\Sigma ={{\mathbb{I}}_{k}}$ and $\beta =2$.

The average in (2.11) is over $d\times k$ standard Gaussian matrices with real $(\beta =1)$, complex $(\beta =2)$ or real quaternion $(\beta =4)$ entries. Such matrices have probability density function proportional to ${{{\rm e}}^{-(\beta /2){\rm Tr}\,(G_{\cdot ,k}^{(\beta ,d)\,\dagger }G_{\cdot ,k}^{(\beta ,d)})}}$, where in the case $\beta =4$ the trace operation takes only one of the independent eigenvalues (recall the eigenvalues in this case are doubly degenerate). In particular, the average is a function of $G_{\cdot ,k}^{(\beta ,d)\,\dagger }G_{\cdot ,k}^{(\beta ,d)}$, suggesting that we change variables by introducing the $k\times k$ Wishart matrix ${{W}^{(\beta ,k)}}=G_{\cdot ,k}^{(\beta ,d)\,\dagger }G_{\cdot ,k}^{(\beta ,d)}$. The corresponding Jacobian is proportional to

$$\begin{eqnarray}&&{{({\rm det} {{W}^{(\beta ,k)}})}^{(\beta /2)(d-k+1-2/\beta )}}\end{eqnarray} \tag{ 2.13 }$$

(see e.g. [7, equation (3.22)]). Making use too of the simple identity

$$\begin{eqnarray*}&&\frac{{{{\rm d}}^{2}}}{{\rm d}{{\mu }^{2}}}{{({\rm det} {{W}^{(\beta ,k)}})}^{\mu }}{{|}_{\mu =0}}={{({\rm log} \;{\rm det} {{W}^{(\beta ,k)}})}^{2}}\end{eqnarray*}$$

we see that the lhs of (2.11) is equal to

$$\begin{eqnarray}&&\frac{1}{4}\frac{{{{\rm d}}^{2}}}{{\rm d}{{\mu }^{2}}}{{\langle {{({\rm det} {{W}^{(\beta ,k)}})}^{\mu }}\rangle }_{{{W}^{(\beta ,k)}}\gt 0}}{{|}_{\mu =0}}.\end{eqnarray} \tag{ 2.14 }$$

This average is over positive definite Hermitian matrices with real ($\beta =1$), complex ($\beta =2$) or real quaternion ($\beta =4$) entries, and distribution having a density function proportional to

$$\begin{eqnarray}&&{{({\rm det} {{W}^{(\beta ,k)}})}^{(\beta /2)(d-k+1-2/\beta )}}{{{\rm e}}^{-(\beta /2){\rm Tr}\,{{W}^{(\beta ,k)}}}}.\end{eqnarray} \tag{ 2.15 }$$

We see that (2.14) is a function only of the eigenvalues of ${{W}^{(\beta ,k)}}$, suggesting that we change variables to the eigenvalues and the eigenvectors. The Jacobian—which has the feature that the eigenvalue and eigenvector terms factorize—is given in e.g. [7, proposition 1.3.4], and we deduce that (2.14) is equal to

$$\begin{eqnarray}&&\frac{1}{4}\frac{{{{\rm d}}^{2}}}{{\rm d}{{\mu }^{2}}}\frac{Z_{k,c}^{(\beta )}}{Z_{k,{{c}_{0}}}^{(\beta )}}{{|}_{\mu =0}},\end{eqnarray} \tag{ 2.16 }$$

where $c=(\beta /2)(d-k+\mu +1-2/\beta )+\mu$, ${{c}_{0}}=(\beta /2)(d-k+1-2/\beta )$ and

$$\begin{eqnarray}&&Z_{k,a}^{(\beta )}\int _{0}^{\infty }{\rm d}{{x}_{1}}\cdots \int _{0}^{\infty }{\rm d}{{x}_{k}}\;\mathop{\prod }\limits_{j=1}^{k}{{{\rm e}}^{-(\beta /2){{x}_{j}}}}x_{j}^{a}\mathop{\prod }\limits_{1\leqslant j\lt l\leqslant k}|{{x}_{l}}-{{x}_{j}}{{|}^{\beta }}.\end{eqnarray} \tag{ 2.17 }$$

The multidimensional integral (2.17) is familiar in random matrix theory as a limiting case of the Selberg integral, and as such is evaluated as a product of gamma functions (see e.g. [7, proposition 4.7.3]). This shows that (2.16) is equal to

$$\begin{eqnarray*}&&\frac{1}{4}\frac{{{{\rm d}}^{2}}}{{\rm d}{{\mu }^{2}}}{{\mathop{\prod }\limits_{j=1}^{k}}^{{}}}{{(\frac{\beta }{2})}^{-\mu }}\frac{\Gamma ((\beta /2)(d-k)+\mu +j\beta /2)}{\Gamma ((\beta /2)(d-k)+j\beta /2)}{{|}_{\mu =0}}.\end{eqnarray*}$$

Evaluating the derivative in terms of the digamma function gives (2.11).

If we now subtract the second term in (2.9), recalling (2.6), we see that the term in the final line of (2.11) cancels, leaving (2.12).□

Equating (2.10) as the evaluation of the lhs of (2.9) for large N with the evaluation (2.12) of the rhs we obtain

$$\begin{eqnarray*}&&\mathop{\sum }\limits_{j=1}^{k}\sigma _{j}^{2}=\frac{1}{4N}\ \mathop{\sum }\limits_{j=1}^{k}{{\Psi }^{{\prime} }}(\beta (d-j+1)/2),\end{eqnarray*}$$

which is equivalent to (2.8) in the case of square matrices.

A minor modification of the above reasoning allows for the derivation of (2.8) in the case of rectangular matrices. An important feature is that the non-negative integer parameters ${{\nu }_{i}}$ ($i=1,...,N$) can only take on a finite number of different values. Let these values be ${{\gamma }_{1}},...,{{\gamma }_{s}}$ and suppose they are in proportion ${{\alpha }_{1}},...,{{\alpha }_{s}}$ where $\sum _{i=1}^{s}{{\alpha }_{i}}=1$. The rhs of the formula (2.9) with $\Sigma ={{\mathbb{I}}_{d}}$ then generalizes to read

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \frac{1}{N}\left( \mathop{\sum }\limits_{i=1}^{s}{{\alpha }_{i}}\langle {{({\rm log} \;{\rm det} {{(G_{\cdot ,k}^{(\beta ,d+{{\gamma }_{i}})\,\dagger }G_{\cdot ,k}^{(\beta ,d+{{\gamma }_{i}})})}^{1/2}})}^{2}}\rangle \right. \\ \quad -\left. \mathop{\sum }\limits_{i=1}^{s}{{\alpha }_{i}}{{(\langle {\rm log} \;{\rm det} {{(G_{\cdot ,k}^{(\beta ,d+{{\gamma }_{i}})\,\dagger }G_{\cdot ,k}^{(\beta ,d+{{\gamma }_{i}})})}^{1/2}}\rangle )}^{2}} \right). \\ \end{array}\end{eqnarray} \tag{ 2.18 }$$

Changing variables ${{W}^{(\beta ,k)}}=G_{\cdot ,k}^{(\beta ,d+{{\gamma }_{i}})\,\dagger }G_{\cdot ,k}^{(\beta ,d+{{\gamma }_{i}})}$ for $i=1,...,s$ we see that the exponent in the Wishart ensemble obtained from the corresponding Jacobian is now $(\beta /2)(d+{{\gamma }_{i}}-k+1-2/\beta )$. The result of lemma 2.1 shows that (2.11) holds with d replaced by $d+\nu$, and all quantities averaged over ν in the sense of (2.7).

A similar averaging of the mean and standard deviation holds in the case that P_N consists of both Gaussian and inverse Gaussian matrices (see [10] for a study of the singular values in this setting in the complex case with N finite). From the definitions, if the Lyapunov exponents for products of $d\times d$ Gaussian random matrices are $\mu _{1}^{+}\gt \mu _{2}^{+}\gt \cdots \gt \mu _{d}^{+}$, then the Lyapunov exponents for products of $d\times d$ inverse Gaussian random matrices are $-\mu _{d}^{+}\gt -\mu _{d-1}^{+}\gt \cdots \gt -\mu _{1}^{+}$. Equivalently, denoting the Lyapunov exponents for products of inverse Gaussian matrices by $\mu _{1}^{-}\gt \mu _{2}^{-}\gt \cdots \gt \mu _{d}^{-}$, we thus have

$$\begin{eqnarray}&&\mu _{i}^{-}=-\mu _{d+1-i}^{+},\end{eqnarray} \tag{ 2.19 }$$

while for the associated variances we have

$$\begin{eqnarray}&&{{(\sigma _{i}^{-})}^{2}}={{(\sigma _{d+1-i}^{+})}^{2}}.\end{eqnarray} \tag{ 2.20 }$$

Consider now a matrix product P_N with proportion ${{\alpha }_{+}}$ Gaussian matrices and ${{\alpha }_{-}}$ inverse Gaussian matrices, and denote the Lyapunov exponents and corresponding variances by $\{\mu _{i}^{(+,-)}\}$, $\{{{(\sigma _{i}^{(+,-)})}^{2}}\}$. Consideration of the derivation of (2.18), and making use of (2.19) and (2.20), we see that

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \mu _{i}^{(+,-)} & = & {{\alpha }_{+}}\mu _{i}^{+}+{{\alpha }_{-}}\mu _{i}^{-}={{\alpha }_{+}}\mu _{i}^{+}-{{\alpha }_{-}}\mu _{d+1-i}^{+}, \\ {{(\sigma _{i}^{(+,-)})}^{2}} & = & {{\alpha }_{+}}{{(\sigma _{i}^{+})}^{2}}+{{\alpha }_{-}}{{(\sigma _{i}^{-})}^{2}}={{\alpha }_{+}}{{(\sigma _{i}^{+})}^{2}}+{{\alpha }_{-}}{{(\sigma _{d+1-i}^{+})}^{2}}. \\ \end{array}\end{eqnarray} \tag{ 2.21 }$$

It has been shown in [8] that for ${{A}_{i}}={{\Sigma }^{1/2}}G_{i}^{(2,d)}$ with $G_{i}^{(2,d)}$ a $d\times d$ standard complex Gaussian random matrix and Σ a fixed positive definite matrix with eigenvalues of ${{\Sigma }^{-1}}$ equal to $\{{{y}_{1}},...,{{y}_{d}}\}$, the Lyapunov spectrum corresponding to the product (1.1) is given by

$$\begin{eqnarray}&&{{\mu }_{k}}=-\frac{1}{2\prod _{1\leqslant i\lt j\leqslant d}^{{}}({{y}_{j}}-{{y}_{i}})}{\rm det} \left[ \begin{array}{ccccccccccccccc} {{[y_{j}^{i-1}]}_{\frac{i=1,...,k-1}{j=1,...,d}}} \\ {{[({\rm log} {{y}_{j}})y_{j}^{k-1}]}_{j=1,...,d}} \\ {{[y_{j}^{i-1}]}_{\frac{i=k+1,...,d}{j=1,...,d}}} \\ \end{array} \right]+\frac{1}{2}\Psi (k).\end{eqnarray} \tag{ 2.22 }$$

This was possible because of the availability of an explicit formula for the average

$$\begin{eqnarray}&&\langle {{({\rm det} {{(G_{\cdot ,k}^{(2,d)\,\dagger }\Sigma G_{\cdot ,k}^{(2,d)})}^{1/2}})}^{\mu }}\rangle .\end{eqnarray} \tag{ 2.23 }$$

Thus, from the works [12, 25] we can read off that this average is equal to (see e.g. [8, equation (2.21)])

$$\begin{eqnarray}&&\frac{{{(-1)}^{k((k+1)/2-d)}}}{\prod {{_{i=1}^{k-1}}^{{}}}i!}\frac{\prod _{j=1}^{d}y_{j}^{k}}{\mathop{\prod }\limits_{1\leqslant j\lt l\leqslant d}}({{y}_{j}}-{{y}_{l}})\ {\rm det} \left[ \begin{array}{ccccccccccccccc} {{[\int _{0}^{\infty }{{t}^{\nu /2+k-i}}{{{\rm e}}^{-{{y}_{j}}t}}\;{\rm d}t]}_{\frac{i=1,...,k}{j=1,...,d}}} \\ {{[y_{j}^{d-i}]}_{\frac{i=k+1,...,d}{j=1,...,d}}} \\ \end{array} \right].\end{eqnarray} \tag{ 2.24 }$$

The Harish-Chandra/Itzykson–Zuber matrix integral (see e.g. [7, proposition 11.6.1]), giving a determinant formula for a particular integral over the unitary group is an essential component of the derivation of (2.24). The analogous matrix integral over the orthogonal or symplectic group does not permit such a structured evaluation, so impeding further progress for calculating the Lyapunov spectrum of the Gaussian ensembles with general Σ and real or real quaternion entries. Differentiating (2.24) with respect to μ and setting $\mu =0$ gives a formula for the right-hand side of (2.9) in the case $\beta =2$. This formula is equivalent to (2.22).

For purposes of calculating the variances associated with the Lyapunov spectrum, following the strategy of section 2.2 we require the analogous evaluation of

$$\begin{eqnarray}&&\langle {{({\rm log} \;{\rm det} {{(G_{\cdot ,k}^{(2,d)\,\dagger }\Sigma G_{\cdot ,k}^{(2,d)})}^{1/2}})}^{2}}\rangle .\end{eqnarray} \tag{ 2.25 }$$

Moreover, according to (2.5) the second term in (2.9) is equal to ${{(\sum _{j=1}^{k}\langle {{\mu }_{j}}\rangle )}^{2}}$, and we know from (2.22) the value of $\{\langle {{\mu }_{j}}\rangle \}$. Although the computation of (2.25) is elementary, simply requiring a double differentiation of (2.24) with respect to μ, unfortunately for general k the resulting formula is not very revealing, consisting of a number of terms involving either summations or double summations up to k, with some of the summands being determinants. We will therefore restrict attention to the cases k = 1, when there is much simplification.

Proposition 2.2. Let $\{{{y}_{1}},...,{{y}_{d}}\}$ be the eigenvalues of ${{\Sigma }^{-1}}$, assumed distinct. The average (2.25) in the case k = 1 is equal to

$$\begin{eqnarray}&&\mu _{1}^{2}+\frac{{{\Psi }^{{\prime} }}(1)}{4}+\frac{1}{4}\mathop{\sum }\limits_{j=1}^{d}\frac{{{({\rm log} {{y}_{j}})}^{2}}}{\prod _{l=1,l\ne j}^{d}{{(1-{{y}_{j}}/{{y}_{l}})}^{{}}}}-\frac{1}{4}{{(\mathop{\sum }\limits_{j=1}^{d}\frac{{\rm log} {{y}_{j}}}{\prod _{l=1,l\ne j}^{d}{{(1-{{y}_{j}}/{{y}_{l}})}^{{}}}})}^{2}}.\end{eqnarray} \tag{ 2.26 }$$

It follows that for large N

$$\begin{eqnarray}&&\sigma _{1}^{2}=\frac{1}{4N}({{\Psi }^{{\prime} }}(1)+\mathop{\sum }\limits_{j=1}^{d}\frac{{{({\rm log} {{y}_{j}})}^{2}}}{\prod _{l=1,l\ne j}^{d}{{(1-{{y}_{j}}/{{y}_{l}})}^{{}}}}-{{(\mathop{\sum }\limits_{j=1}^{d}\frac{{\rm log} {{y}_{j}}}{\prod _{l=1,l\ne j}^{d}{{(1-{{y}_{j}}/{{y}_{l}})}^{{}}}})}^{2}}).\end{eqnarray} \tag{ 2.27 }$$

Proof. Our task is to differentiate (2.24) with k = 1 twice with respect to μ and to set $\mu =0$. Since the only dependence on ν for k = 1 is in the first row, this reduces to computing the second derivative of each term in the first row. Noting

$$\begin{eqnarray*}&&\frac{{{{\rm d}}^{2}}}{{\rm d}{{\mu }^{2}}}\int _{0}^{\infty }{{t}^{\mu /2}}{{{\rm e}}^{-{{y}_{j}}t}}\;{\rm d}t{{|}_{\mu =0}}=\frac{1}{4{{y}_{j}}}({{({\rm log} {{y}_{j}})}^{2}}-2({\rm log} {{y}_{j}})\Psi (1)+{{\Psi }^{{\prime} }}(1)+{{(\Psi (1))}^{2}})\end{eqnarray*}$$

and substituting for the first row of (2.24) we see after some simple manipulation that the sought average is equal to

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \frac{1}{4\prod _{1\leqslant j\lt l\leqslant d}^{{}}({{y}_{l}}-{{y}_{j}})}({\rm det} \left[ \begin{array}{ccccccccccccccc} {{[{{({\rm log} {{y}_{j}})}^{2}}]}_{j=1,...,d}} \\ {{[y_{j}^{i-1}]}_{\frac{i=2,...,d}{j=1,...,d}}} \\ \end{array} \right]-2\Psi (1){\rm det} \left[ \begin{array}{ccccccccccccccc} {{[{\rm log} {{y}_{j}}]}_{j=1,...,d}} \\ {{[y_{j}^{i-1}]}_{\frac{i=2,...,d}{j=1,...,d}}} \\ \end{array} \right]) \\ \quad +\frac{{{\Psi }^{{\prime} }}(1)+{{(\Psi (1))}^{2}}}{4}, \\ \end{array}\end{eqnarray} \tag{ 2.28 }$$

where in eliminating the determinant which would otherwise appear on the second line use has been made of the Vandermonde determinant evaluation

$$\begin{eqnarray}&&{\rm det} {{[y_{j}^{i-1}]}_{i=1,...,d}}=\mathop{\prod }\limits_{1\leqslant j\lt l\leqslant d}({{y}_{l}}-{{y}_{j}}).\end{eqnarray} \tag{ 2.29 }$$

Recalling (2.22) in the case k = 1 allows (2.30) to be equivalently written

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \mu _{1}^{2}+\frac{{{\Psi }^{{\prime} }}(1)}{4}+\frac{1}{4\prod _{1\leqslant j\lt l\leqslant d}^{{}}({{y}_{l}}-{{y}_{j}})}{\rm det} \left[ \begin{array}{ccccccccccccccc} {{[{{({\rm log} {{y}_{j}})}^{2}}]}_{j=1,...,d}} \\ {{[y_{j}^{i-1}]}_{\frac{i=2,...,d}{j=1,...,d}}} \\ \end{array} \right] \\ -\frac{1}{4\prod _{1\leqslant j\lt l\leqslant d}^{{}}{{({{y}_{l}}-{{y}_{j}})}^{2}}}{{({\rm det} \left[ \begin{array}{ccccccccccccccc} {{[{\rm log} {{y}_{j}}]}_{j=1,...,d}} \\ \quad {{[y_{j}^{i-1}]}_{\frac{i=2,...,d}{j=1,...,d}}} \\ \end{array} \right])}^{2}}. \\ \end{array}\end{eqnarray} \tag{ 2.30 }$$

Expanding the determinants by the first row and making use of (2.29) to evaluate the cofactors we obtain (2.26).

To deduce (2.27) from (2.26) we subtract $\mu _{1}^{2}$ in keeping with (2.9) and equate with (2.10) in the case k = 1.□

Our derivation of the formula (2.27) for $\sigma _{1}^{2}$ in the case of general variance complex Gaussian matrices has made use of the determinant formula (2.24), which as already mentioned has no analogue in the real or real quaternion cases. This was similarly the case for the derivation given in [8] of the Lyapunov exponents (2.22). However, with regards to the latter with k = 1 (largest Lyapunov exponent) it has recently been shown by Kargin [18] that an alternative calculation is possible, which furthermore does generalize to the real $(\beta =1)$ and real quaternion $(\beta =4)$ cases, giving the formula [18, theorem 1.1]

$$\begin{eqnarray}&&2{{\mu }_{1}}=-\gamma +{\rm log} (\frac{2}{\beta })+\int _{0}^{\infty }({{\chi }_{x\in (0,1)}}-\mathop{\prod }\limits_{i=1}^{d}{{(1+{{\frac{x}{{{y}_{i}}}}^{{}}})}^{-\beta /2}})\;\frac{{\rm d}x}{x}.\end{eqnarray} \tag{ 2.31 }$$

In (2.31) γ denotes Euler's constant and ${{\chi }_{J}}=1$ for J true and ${{\chi }_{J}}=0$ otherwise. Here we adapt the method of [18] to the computation of the rhs of (2.9) with k = 1.

Proposition 2.3. Let $\mathcal{C}$ denote the simple closed contour starting at large $R\gt {\rm max} \;\{{{y}_{i}}\}$, goes along the upper edge of the real axis to $\epsilon \gt 0$, where $\epsilon \lt {\rm min} \;\{{{y}_{i}}\}$, crosses to the lower edge of the real axis and returns to R. Define

$$\begin{eqnarray}&&{{J}_{p}}=\frac{1}{2\pi i}{{\int }_{\mathcal{C}}}{{({\rm log} z)}^{p}}\mathop{\prod }\limits_{i=1}^{d}{{(1-{{\frac{z}{{{y}_{i}}}}^{{}}})}^{-\beta /2}}\;\frac{{\rm d}z}{z}.\end{eqnarray} \tag{ 2.32 }$$

We have

$$\begin{eqnarray}&&\sigma _{1}^{2}=\frac{1}{4N}(\frac{{{\pi }^{2}}}{6}-{{J}_{2}}-J_{1}^{2}).\end{eqnarray} \tag{ 2.33 }$$

Proof. For k = 1 the quantity $(G_{\cdot ,k}^{(\beta ,d)\,\dagger }\Sigma G_{\cdot ,k}^{(\beta ,d)})$ is a positive scalar random variable. Let its density function be denoted ${{p}_{\beta }}(\lambda )$. We then have that

$$\begin{eqnarray}&&\langle {{({\rm log} \;{\rm det} {{(G_{\cdot ,k}^{(\beta ,d)\,\dagger }\Sigma G_{\cdot ,k}^{(\beta ,d)})}^{1/2}})}^{2}}\rangle =\frac{1}{4}\int _{0}^{\infty }{{({\rm log} \lambda )}^{2}}{{p}_{\beta }}(\lambda )\;{\rm d}\lambda .\end{eqnarray} \tag{ 2.34 }$$

It is shown in [18] that

$$\begin{eqnarray}&&{{p}_{\beta }}(\lambda )=\frac{1}{{{\tilde{c}}_{\beta }}}\frac{1}{2\pi i}{{\int }_{\mathcal{C}}}{{{\rm e}}^{-\beta \lambda z/2}}\mathop{\prod }\limits_{i=1}^{d}{{(z-{{y}_{i}})}^{-\beta /2}}\;{\rm d}z,\qquad {{\tilde{c}}_{\beta }}=-\frac{2}{\beta }\mathop{\prod }\limits_{i=1}^{d}{{(-{{y}_{i}})}^{-\beta /2}},\end{eqnarray} \tag{ 2.35 }$$

where $\mathcal{C}$ is as in (2.32). Substituting (2.35) in (2.34), interchanging the order of integration and noting

$$\begin{eqnarray*}&&\int _{0}^{\infty }{{({\rm log} \lambda )}^{2}}{{{\rm e}}^{-\beta \lambda z/2}}\;{\rm d}\lambda =\frac{2}{\beta z}({{({\rm log} \frac{2}{\beta z})}^{2}}-2\gamma ({\rm log} \frac{2}{\beta z})+{{\gamma }^{2}}+\frac{{{\pi }^{2}}}{6})\end{eqnarray*}$$

shows that

$$\begin{eqnarray}&&\int _{0}^{\infty }{{({\rm log} \lambda )}^{2}}{{p}_{\beta }}(\lambda )\;{\rm d}\lambda =({{\gamma }^{2}}+\frac{{{\pi }^{2}}}{6}){{I}_{0}}-2\gamma {{I}_{1}}+{{I}_{2}},\end{eqnarray} \tag{ 2.36 }$$

where

$$\begin{eqnarray}&&{{I}_{p}}\frac{1}{{{\tilde{c}}_{\beta }}}\frac{1}{2\pi i}{{\int }_{\mathcal{C}}}(\frac{2}{\beta z}){{({\rm log} \frac{2}{\beta z})}^{p}}\mathop{\prod }\limits_{i=1}^{d}{{(z-{{y}_{i}})}^{-\beta /2}}\;{\rm d}z.\end{eqnarray} \tag{ 2.37 }$$

By deforming the contour to go around the pole at z = 0 we deduce that ${{I}_{0}}=1$. We read off from [18, equation (13)] that

$$\begin{eqnarray}&&{{I}_{1}}=2{{\mu }_{1}}+\gamma .\end{eqnarray} \tag{ 2.38 }$$

Also, simple manipulation shows

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{I}_{2}} & = & {{({\rm log} \frac{2}{\beta })}^{2}}{{I}_{0}}+2({\rm log} \frac{2}{\beta })({{I}_{1}}-{\rm log} \frac{2}{\beta })-{{J}_{2}} \\ {} & = & {{({\rm log} \frac{2}{\beta })}^{2}}+2({\rm log} \frac{2}{\beta })(2{{\mu }_{1}}+\gamma -{\rm log} \frac{2}{\beta })-{{J}_{2}}, \\ \end{array}\end{eqnarray} \tag{ 2.39 }$$

where J_p is specified by (2.32). We substitute (2.38) and (2.39) in (2.36) then subtract ${{(2{{\mu }_{1}})}^{2}}$, using the fact that [18, equation (14)]

$$\begin{eqnarray}&&2{{\mu }_{1}}=-\gamma +{\rm log} \frac{2}{\beta }-{{J}_{1}}\end{eqnarray} \tag{ 2.40 }$$

to arrive at (2.33).□

In the case $\beta =2$ the integrals in (2.33) can be evaluated using the residue theorem to reclaim (2.27), noting too that ${{\Psi }^{{\prime} }}(1)={{\pi }^{2}}/6$. In the case that ${{y}_{1}}=\cdots ={{y}_{d}}=1$ and thus $\Sigma ={{\mathbb{I}}_{d}}$, provided $\beta d/2\in {{\mathbb{Z}}^{+}}$ use of the residue theorem shows that

$$\begin{eqnarray}&&-{{J}_{1}}=\mathop{\sum }\limits_{s=1}^{d\beta /2-1}\frac{1}{s},\qquad -{{J}_{2}}=\mathop{\sum }\limits_{s=2}^{d\beta /2-1}{{\alpha }_{s}}\quad {\rm where}\quad {{\alpha }_{s}}=\frac{2}{s}\mathop{\sum }\limits_{p=1}^{s-1}\frac{1}{p},\end{eqnarray} \tag{ 2.41 }$$

which when substituted in (2.33) can be verified to give the same value of $\sigma _{1}^{2}$ as (2.8) in the case i = 1. It is shown in [18] that

$$\begin{eqnarray}&&{{J}_{1}}=-\int _{0}^{\infty }({{\chi }_{x\in (0,1)}}-\mathop{\prod }\limits_{i=1}^{d}{{(1+{{\frac{x}{{{y}_{i}}}}^{{}}})}^{-\beta /2}})\;\frac{{\rm d}x}{x},\end{eqnarray} \tag{ 2.42 }$$

as is consistent with the equality between (2.31) and (2.40). Analogous working shows

$$\begin{eqnarray}&&{{J}_{2}}=2\int _{0}^{\infty }({{\chi }_{x\in (0,1)}}-\mathop{\prod }\limits_{i=1}^{d}{{(1+{{\frac{x}{{{y}_{i}}}}^{{}}})}^{-\beta /2}})\frac{{\rm log} x}{x}\;{\rm d}x+\frac{{{\pi }^{2}}}{3},\end{eqnarray} \tag{ 2.43 }$$

so it is also possible to express (2.33) in terms of real integrals.

As remarked in the opening paragraph of section 2.1, the tractability of the computation of the Lyapunov spectrum for the Gaussian random matrices G is due to Gaussian matrices belonging to the class of isotropic matrices, meaning that any vector $\vec{x}$, $G\vec{x}/|\vec{x}|$ has distribution independent of $\vec{x}$. This property holds for any class of random matrices that are unitary invariant, and thus their distribution in unchanged by multiplication by a unitary matrix U, which is furthermore required to have entries of the same type (real, complex or real quaternion) as the random matrices.

Let ${{Z}^{(\beta ,d+n)}}$ be a random unitary matrix of size $(d+n)\times (d+n)$ with real ($\beta =1$, complex ($\beta =2$) or real quaternion ($\beta =4$) elements, and chosen with Haar measure. The latter has the defining feature that it is unchanged by multiplication of the matrices by fixed unitary matrices of the appropriate type. Suppose the fixed unitary matrices have the block structure

$$\begin{eqnarray*}U=\left[ \begin{array}{ccccccccccccccc} {{V}_{d\times d}} & {} \\ {} & {{\mathbb{I}}_{n}} \\ \end{array} \right],\end{eqnarray*}$$

where ${{V}_{d\times d}}$ is a $d\times d$ unitary matrix. Such matrices act only on the top $d\times d$ sub-matrix of ${{Z}^{(\beta ,d+n)}}$, so it must be that square submatrices of Haar distributed unitary matrices (also referred to as truncated unitary matrices [26]) are themselves unitary invariant, and are thus examples of isotropic matrices. With the matrices A_i in (1.1) chosen from the ensemble of truncated unitary matrices, our interest in this section is to compute the Lyapunov spectrum as well as the associated variances. Crucial to our computations is the fact that a $d\times k$ ($d\geqslant k$) sub-block $Z_{\cdot ,k}^{(\beta ,d)}$ is distributed with a probability density function proportional to

$$\begin{eqnarray}&&{{({\rm det} ({{\mathbb{I}}_{k}}-Z_{\cdot ,k}^{(\beta ,d)\,\dagger }Z_{\cdot ,k}^{(\beta ,d)}))}^{{{c}_{2}}}},\qquad {{c}_{2}}(\beta /2)(n-k+1-2/\beta ).\end{eqnarray} \tag{ 3.1 }$$

We remark that the singular values of products of truncated complex unitary matrices for finite N are the topic of the recent work [19].

According to (2.5) we have

$$\begin{eqnarray}&&{{\mu }_{1}}+\cdots +{{\mu }_{k}}=\langle {\rm log} \;{\rm det} {{(Z_{\cdot ,k}^{(\beta ,d)\,\dagger }Z_{\cdot ,k}^{(\beta ,d)})}^{1/2}}\rangle ,\end{eqnarray} \tag{ 3.2 }$$

where the average is with respect to (3.1). The rhs of (3.2) is easy to compute, allowing for a simple formula for the individual ${{\mu }_{i}}$.

Proposition 3.1. Let the matrices A_i in (1.1) be chosen as the top $d\times d$ sub-block of a $(d+n)\times (d+n)$ Haar distributed random unitary matrix with real $(\beta =1)$, complex $(\beta =2)$ or real quaternion $(\beta =4)$ entries. For $n\geqslant d$ at least we have

$$\begin{eqnarray}&&{{\mu }_{i}}=\frac{1}{2}(\Psi ((\beta /2)(d-i+1))-\Psi ((\beta /2)(n+d-i+1))).\end{eqnarray} \tag{ 3.3 }$$

Proof. As in the Gaussian case, an essential point is that the average (3.2) is a function only of the positive definite matrix ${{W}^{(\beta ,k)}}=Z_{\cdot ,k}^{(\beta ,d)\,\dagger }Z_{\cdot ,k}^{(\beta ,d)}$. Since $Z_{\cdot ,k}^{(\beta ,d)}$ is a sub-matrix of a unitary matrix, ${{W}^{(\beta ,k)}}$ is further constrained to have eigenvalues less than 1, which we denote by writing ${{W}^{(\beta ,k)}}\lt 1$. Regarding this as a change of variables, the corresponding Jacobian is proportional to (2.13). Hence our task becomes that of evaluating

$$\begin{eqnarray}&&\frac{1}{2}\frac{{\rm d}}{{\rm d}\mu }{{\langle {{({\rm det} {{W}^{(\beta ,k)}})}^{\mu }}\rangle }_{0\lt {{W}^{(\beta ,k)}}\lt 1}}{{|}_{\mu =0}},\end{eqnarray} \tag{ 3.4 }$$

where ${{W}^{(\beta ,k)}}$ has distribution with a probability density function proportional to

$$\begin{eqnarray}&&{{({\rm det} {{W}^{(\beta ,k)}})}^{{{c}_{1}}}}({\rm det} {{({{\mathbb{I}}_{k}}-{{W}^{(\beta ,k)}})}^{{{c}_{2}}}}\end{eqnarray} \tag{ 3.5 }$$

with c₂ as in (3.1) and

$$\begin{eqnarray}&&{{c}_{1}}=(\beta /2)(d-k+1-2/\beta ).\end{eqnarray} \tag{ 3.6 }$$

Both (3.4) and (3.5) are functions of the eigenvalues of ${{W}^{(\beta ,k)}}$ only. Introducing the appropriate Jacobian, as can be found in e.g. [7, proposition 1.3.4], gives for (3.4) the expression in terms of multiple integrals

$$\begin{eqnarray}&&\frac{1}{2}\frac{{\rm d}}{{\rm d}\mu }\frac{I_{k,{{c}_{1}}+\mu ,{{c}_{2}}}^{(\beta )}}{I_{k,{{c}_{1}},{{c}_{2}}}^{(\beta }}{{|}_{\mu =0}},\end{eqnarray} \tag{ 3.7 }$$

where

$$\begin{eqnarray}&&I_{k,{{\alpha }_{1}},{{\alpha }_{2}}}^{(\beta )}=\int _{0}^{1}{\rm d}{{x}_{1}}\cdots \int _{0}^{1}{\rm d}{{x}_{k}}\;\mathop{\prod }\limits_{l=1}^{k}x_{l}^{{{\alpha }_{1}}}{{(1-{{x}_{l}})}^{{{\alpha }_{2}}}}\mathop{\prod }\limits_{1\leqslant j\lt l\leqslant k}|{{x}_{l}}-{{x}_{j}}{{|}^{\beta }}.\end{eqnarray} \tag{ 3.8 }$$

This is precisely the Selberg integral (see e.g. [7, chapter 4]), already seen in a limiting form in (2.17). Inserting its gamma functions evaluation reduces (3.7) to

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \frac{1}{2}\frac{{\rm d}}{{\rm d}\mu }\mathop{\prod }\limits_{j=0}^{k-1}\frac{\Gamma ({{c}_{1}}+\mu +1+j\beta /2)\Gamma ({{c}_{1}}+{{c}_{1}}+{{c}_{2}}+2+(k+j-1)\beta /2)}{\Gamma ({{c}_{1}}+1+j\beta /2)\Gamma ({{c}_{1}}+\mu +{{c}_{2}}+2+(k+j-1)\beta /2)}{{|}_{\mu =0}} \\ \quad =\frac{1}{2}(\mathop{\sum }\limits_{j=0}^{k-1}\Psi ({{c}_{1}}+1+j\beta /2)-\mathop{\sum }\limits_{j=0}^{k-1}\Psi ({{c}_{1}}+{{c}_{2}}+2+(k+j-1)\beta /2)). \\ \end{array}\end{eqnarray} \tag{ 3.9 }$$

Equating this with the lhs of (3.2) and inserting the values of c₁ and c₂ from (3.6) and (3.1) we arrive at (3.3).□

We remark that even though our derivation of (3.3) requires $n\geqslant d$, used for example in the convergence of the integrals in (3.7), the final expression is well defined for $n\geqslant 0$, and we expect it to remain valid. A similar effect holds in other calculations relating to truncated unitary matrices which are based on (3.5) [11]. As evidence for our claim, note that for n = 0, (3.3) gives that each Lyapunov exponent is equal to 0. This is consistent with the matrices in the product (1.1) then themselves being unitary matrices, as the $d\times d$ sub-matrix of ${{Z}^{(\beta ,d+n)}}$ is for n = 0 the matrix itself. Furthermore, for specific $n\lt d$ the predictions given by (3.3) are found to be consistent with numerical computations, using the methods of section 4.1 below (for example with n = 1, d = 2 a Monte-Carlo estimate with $N={{10}^{6}}$ trials gave ${{\mu }_{1}}\approx -0.2509$ while (3.3) gives ${{\mu }_{1}}=-\frac{1}{4}$).

Similar working gives the variances associated with the Lyapunov exponents (3.3).

Proposition 3.2. The analogue of the formula (2.8) for the variances associated with the Lyapunov exponents in the case of products of the truncated unitary matrices of proposition 3.1 is given by

$$\begin{eqnarray}&&\sigma _{i}^{2}=\frac{1}{4N}({{\Psi }^{{\prime} }}((\beta /2)(d-i+1))-{{\Psi }^{{\prime} }}((\beta /2)(n+d-i+1)).\end{eqnarray} \tag{ 3.10 }$$

Proof. We use the analogue of the formula (2.9) with the lhs replaced by its leading large N form (2.10). Thus our main task is to compute

$$\begin{eqnarray}&&\langle {{({\rm log} \;{\rm det} {{(Z_{\cdot ,k}^{(\beta ,d)\,\dagger }Z_{\cdot ,k}^{(\beta ,d)})}^{1/2}})}^{2}}\rangle .\end{eqnarray} \tag{ 3.11 }$$

Following the working in the proof of proposition 3.1 we see that this average is equal to

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {} & {} & \frac{1}{4}\frac{{{{\rm d}}^{2}}}{{\rm d}{{\mu }^{2}}}\frac{I_{k,{{c}_{1}}+\mu ,{{c}_{2}}}^{(\beta )}}{I_{k,{{c}_{1}},{{c}_{2}}}^{(\beta }}{{|}_{\mu =0}} \\ {} & = & \frac{1}{4}\frac{{{{\rm d}}^{2}}}{{\rm d}{{\mu }^{2}}}\mathop{\prod }\limits_{j=0}^{k-1}\frac{\Gamma ({{c}_{1}}+\mu +1+j\beta /2)\Gamma ({{c}_{1}}+{{c}_{1}}+{{c}_{2}}+2+(k+j-1)\beta /2)}{\Gamma ({{c}_{1}}+1+j\beta /2)\Gamma ({{c}_{1}}+\mu +{{c}_{2}}+2+(k+j-1)\beta /2)}{{|}_{\mu =0}} \\ {} & = & \frac{1}{4}(\mathop{\sum }\limits_{j=0}^{k-1}{{\Psi }^{{\prime} }}({{c}_{1}}+1+j\beta /2)-\mathop{\sum }\limits_{j=0}^{k-1}{{\Psi }^{{\prime} }}({{c}_{1}}+{{c}_{2}}+2+(k+j-1)\beta /2)) \\ {} & {} & +\frac{1}{4}{{(\mathop{\sum }\limits_{j=0}^{k-1}\Psi ({{c}_{1}}+1+j\beta /2)-\mathop{\sum }\limits_{j=0}^{k-1}\Psi ({{c}_{1}}+{{c}_{2}}+2+(k+j-1)\beta /2))}^{2}}. \\ \end{array}\end{eqnarray*}$$

We recognize from (3.9) that the final line is equal to ${{(\sum _{i=1}^{d}{{\mu }_{i}})}^{2}}$, and so cancels out of the rhs of (2.9). Thus

$$\begin{eqnarray*}&&\mathop{\sum }\limits_{i=1}^{k}\sigma _{i}^{2}=\frac{1}{4}(\mathop{\sum }\limits_{j=0}^{k-1}{{\Psi }^{{\prime} }}({{c}_{1}}+1+j\beta /2)-\mathop{\sum }\limits_{j=0}^{k-1}{{\Psi }^{{\prime} }}({{c}_{1}}+{{c}_{2}}+2+(k+j-1)\beta /2)).\end{eqnarray*}$$

Recalling the values of c₁ and c₂ from (3.6) and (3.1) we see that this is equivalent to the stated formula (3.10). □

The formulas (3.3) and (3.10) have some similarity with the corresponding formulas (2.7) and (2.8) for products of rectangular Gaussian random matrices, one of size $d\times d$ (and thus square), the other of size $(d+n)\times d$. In particular the formulas agree if the $d\times d$ matrices were to occur in 'proportion' 1, and the $(d+n)\times d$ matrices in 'proportion' −1.

In our earlier paper on Lyapunov exponents [8] a discussion and references were given on the difficulty of accurately computing the Lyapunov exponents. A more robust, but less accurate approach is to use Monte-Carlo simulation based on (2.1). Numerically stable ways to carry out such computations have been reported in [6]. Thus, in relation to ${{\mu }_{1}}$, and with d = 2, ${{\vec{x}}_{0}}={{\left[ \begin{array}{ccccccccccccccc} 1 & 0 \\ \end{array} \right]}^{ \sf T}}$, we iteratively define ${{\vec{x}}_{i}}={{A}_{i}}{{\vec{y}}_{i-1}}$ $(i=1,2,...)$, where ${{\vec{y}}_{i}}={{\vec{x}}_{i}}/|{{\vec{x}}_{i}}|$. According to (2.1), a Monte-Carlo estimate of the largest Lyapunov exponent is then given by

$$\begin{eqnarray}&&{{\mu }_{1}}=\frac{1}{N}\mathop{\sum }\limits_{j=1}^{N}{\rm log} |{{\vec{x}}_{j}}|,\end{eqnarray} \tag{ 4.1 }$$

and similarly an estimate of the corresponding variance is given by

$$\begin{eqnarray}&&\sigma _{1}^{2}=\frac{1}{N}(\frac{1}{N}\mathop{\sum }\limits_{j=1}^{N}{{({\rm log} |{{\vec{x}}_{j}}|)}^{2}}-\mu _{1}^{2}).\end{eqnarray} \tag{ 4.2 }$$

To use (2.1) to compute ${{\mu }_{1}}+{{\mu }_{2}}$, again taking d = 2 for simplicity, let $\vec{x}_{0}^{(1)}={{\left[ \begin{array}{ccccccccccccccc} 1 & 0 \\ \end{array} \right]}^{ \sf T}}$ and $\vec{x}_{0}^{(2)}={{\left[ \begin{array}{ccccccccccccccc} 0 & 1 \\ \end{array} \right]}^{ \sf T}}$, then define $\vec{x}_{i}^{(p)}={{A}_{i}}\vec{y}_{i}^{(p)}$, ($p=1,2$, $i=1,2,...$) where $\{\vec{y}_{i}^{(1)},\vec{y}_{i}^{(2)}\}$ is obtained from $\{\vec{x}_{i-1}^{(1)},\vec{x}_{i-1}^{(2)}\}$ by the Gram–Schmidt orthonormalization procedure. It follows from (2.1) that a Monte-Carlo estimate of ${{\mu }_{1}}+{{\mu }_{2}}$ is

$$\begin{eqnarray}&&{{\mu }_{1}}+{{\mu }_{2}}=\frac{1}{N}\mathop{\sum }\limits_{j=1}^{N}{\rm log} \;{\rm det} {{(X_{j}^{\dagger }{{X}_{j}})}^{1/2}},\end{eqnarray} \tag{ 4.3 }$$

where X_i is the $2\times 2$ matrix with columns given by $\vec{x}_{i}^{(1)}$ and $\vec{x}_{i}^{(2)}$, while an estimate of the corresponding variance is

$$\begin{eqnarray}&&\sigma _{1}^{2}+\sigma _{2}^{2}=\frac{1}{N}(\frac{1}{N}\mathop{\sum }\limits_{j=1}^{N}{{({\rm log} \;{\rm det} {{(X_{j}^{\dagger }{{X}_{j}})}^{1/2}})}^{2}}-{{({{\mu }_{1}}+{{\mu }_{2}})}^{2}}).\end{eqnarray} \tag{ 4.4 }$$

As an example, let's consider the case of (1.1) with ${{A}_{i}}={{\Sigma }^{1/2}}G_{i}^{(2,2)}$. We choose the eigenvalues of ${{\Sigma }^{-1}}$ to equal $1,1/4$ so that

$$\begin{eqnarray*}{{A}_{i}}=\left[ \begin{array}{ccccccccccccccc} 1 & 0 \\ 0 & 2 \\ \end{array} \right]G_{i}^{(2,2)}.\end{eqnarray*}$$

The formula (2.27) gives $N\sigma _{1}^{2}=(\frac{{{\pi }^{2}}}{24}-\frac{4}{9}{{({\rm log} 2)}^{2}})=0.197699...$, whereas a Monte-Carlo calculation based on (4.2) with $N={{10}^{6}}$ gave $N\sigma _{1}^{2}\approx 0.19686$.

To perform Monte-Carlo simulations in relation to products of truncated unitary matrices, we must first generate random unitary matrices with Haar measure. One way to do this is to apply the Gram–Schmidt orthogonalization producure to a standard Gaussian matrix. In the complex case, with $d=n=2$ we did this and with $N={{10}^{5}}$ found the estimates ${{\mu }_{1}}\approx -0.41672$ and ${{\mu }_{1}}+{{\mu }_{2}}\approx -1.6685$ whereas according to proposition 3.1 ${{\mu }_{1}}=-\frac{5}{12}=-0.4166...$ and ${{\mu }_{1}}+{{\mu }_{2}}=-\frac{7}{6}=-1.666...$. For the corresponding variances we found $N\sigma _{1}^{2}\approx 0.0894$ and $N(\sigma _{1}^{2}+\sigma _{2}^{2})\approx 0.3952$ whereas proposition 3.2 gives $N\sigma _{1}^{2}=13/144=0.0902...$ and $N(\sigma _{1}^{2}+\sigma _{2}^{2})=29/72\approx 0.4027...$.

As revised in the first paragraph of the Introduction, the Lyapunov exponents are defined in terms of the eigenvalues of the positive definite matrix $P_{N}^{\dagger }{{P}_{N}}$, where P_N is the random matrix product (1.1). In the recent works [3, 17], for the particular cases that the A_i in (1.1) are standard Gaussian matrices with real, complex or real quaternion entries, explicit asymptotic analysis of the joint distribution of the spectrum of P_N has revealed that with the eigenvalues written in the form ${{z}_{k}}={{{\rm e}}^{N{{\lambda }_{k}}+i{{\theta }_{k}}}}$, the $\{{{\lambda }_{k}}\}$—referred to as the stability exponents—are identical to the Lyapunov exponents. Moreover in [3] a proof of the equality between the Lyapunov and stability exponents in the case of $2\times 2$ random matrices from an isotropic ensemble has been given. It was then conjectured that all matrix products formed from an isotropic ensemble have this property.

Actually the question of the asymptotic equality between the Lyapunov and stability exponents has been investigated at least as far back as 1987 [13], albeit in the setting of the linearization of sets of differential equations about fixed points. In this work it was conjectured that the exponents agree whenever the spectrum of P_N is non-degenerate. The reference book [6] cites [13], and furthermore addresses directly the case of products of random matrices, writing 'if the matrix P_N can be written in diagonal form, then (in our notation) ${{\lambda }_{k}}={{\mu }_{k}}$'. One immediate consequence of this conjecture is the prediction of the Lyapunov spectrum for a product of upper triangular random matrices. There, up to ordering $\{{{\lambda }_{k}}\}=\{{\rm log} \langle |{{A}_{kk}}|\rangle \}$ where A_kk denotes the entry $(kk)$ in the matrices A_i of (1.1), so the conjecture gives this same set for the Lyapunov exponents. In the case of ${{\mu }_{1}}$ this agrees with a known theorem [23].

Following [14], we emphasize that the conjecture ${{\lambda }_{k}}={{\mu }_{k}}$ appears to be an intrinsic property of the asymptotic limit of random matrix products, and is not shared by other asymptotic limits of spectra in random matrix theory. For example, in the case of $N\times N$ standard real, complex or real quaternion standard Gaussian random matrices

$$\begin{eqnarray*}&&\mathop{{\rm lim} }\limits_{N\to \infty }\frac{{{\nu }_{1}}({{({{X}^{\dagger }}X)}^{1/2}})}{|{{\nu }_{1}}(X)|}=\sqrt{2},\end{eqnarray*}$$

where ${{\nu }_{1}}(Y)$ denotes the eigenvalue of largest modulus of Y. That this ratio is at least unity follows from the general inequality $|{{\nu }_{1}}(Y)|\leqslant {{\nu }_{1}}({{({{Y}^{\dagger }}Y)}^{1/2}})$, which in turn is a consequence of the sub-multiplicity of the rhs viewed as a matrix norm.