ASYMPTOTIC PRODUCTS OF INDEPENDENT GAUSSIAN RANDOM MATRICES WITH CORRELATED ENTRIES

In this work we address the problem of determining the asymptotic spectral measure of the product of independent, Gaussian random matrices with correlated entries, as the dimension and the number of multiplicative terms goes to inﬁnity. More speciﬁcally, let { X p ( N ) } ∞ p = 1 be a sequence of N × N independent random matrices with independent and identically distributed Gaussian entries of zero mean and variance 1 (cid:112) N . Let { Σ( N ) } ∞ N = 1 be a sequence of N × N deterministic and Hermitian matrices such that the sequence converges in moments to a compactly supported probability measure σ . Deﬁne the random matrix Y p ( N ) as Y p ( N ) = X p ( N )Σ( N ) . This is a random matrix with correlated Gaussian entries and covariance matrix E ( Y p ( N ) ∗ Y p ( N )) = Σ( N ) 2 for every p ≥ 1. The positive deﬁnite N × N matrix converges in distribution to a compactly supported measure in [ 0, ∞ ) as the dimension of the matrices N → ∞ . We show that the sequence of measures ν n converges in distribution to a compactly supported measure ν n → ν as n → ∞ . The measures ν n and ν only depend on the measure σ . Moreover, we deduce an exact closed-form expression for the measure ν as a function of the measure σ


Introduction
Considerable effort has been invested over the last century in determining the spectral properties of ensembles of matrices with randomly chosen elements and in discovering the remarkably broad applicability of these results to systems of physical interest. In the last decades, a considerable amount of work has emerged in the communications and information theory on the fundamental limits of communication channels that makes use of results in random matrix theory. In spite of a similarly rich set of potential applications e.g., in the statistical theory of Markov processes and in various chaotic dynamical systems in classical physics, the properties of products of random 353 DOI: 10.1214/ECP.v16-1635 matrices have received considerably less attention. See [9] for a survey of products of random matrices in statistics and [7] for a review of physics applications.
In this work we consider the problem of determining the asymptotic spectral measure of the product of random matrices. More specifically, let {X p (N )} ∞ p=1 be a sequence of N ×N independent random matrices with independent, and identically distributed Gaussian entries of zero mean and variance 1 N . Let {Σ(N )} ∞ N =1 be a sequence of N × N deterministic and Hermitian matrices such that the sequence converges in moments to a compactly supported probability measure σ. More precisely, for every k ≥ 1 the limit lim N →∞ tr N (Σ k (N )) = R t k dσ(t) (1) where tr N (·) is the normalized trace (sum of the diagonal elements divided by N ). Define the random matrix Y p (N ) as Y p (N ) = X p (N )Σ(N ). This is a random N × N matrix with correlated Gaussian entries with covariance matrix E(Y p (N ) * Y p (N )) = Σ(N ) 2 for every p ≥ 1. We will sometimes drop the index N for notation simplicity. We will show that the positive definite N × N matrix converges in distribution to a compactly supported measure in [0, ∞) as the dimension of the matrices N → ∞. Moreover, the sequence of measures ν n converges weakly to a compactly supported measure ν n → ν.
The measures ν n and ν only depend on the measure σ. Moreover, we deduce a exact closed-form expression of the measure ν as a function of the measure σ. In particular, this gives us a map ∆ : → + from the compactly supported measure in the real line to the compactly supported measures in [0, ∞). We would like also to mention that the normalization 1 2n in the matrix B n is necessary for convergence, and it is indeed the appropriate one. We can think of this result as a multiplicative version of the central limit theorem for random matrices. The case where the matrices Σ(N ) change with p and the corresponding limit laws of eigenvalues σ p are allowed to change from the different values of p is an interesting case to study. However, in this case the matrices in the product are not identically distributed and the problem is a little bit more involved. This question is out of the scope of this paper and we leave it for a subsequent work.
The Lyapunov exponents play an important role in a number of different contexts including the study of the Ising model, the Hausdorff dimension of measures, probability theory and dynamical systems. Recently there has been renewed interest because of their usefulness in the study of the entropy rates of Markov models. It is a fundamental problem to find an explicit expression for the Lyapunov exponents. Unfortunately, there are very few analytical techniques available for their study. Traditionally, they have been approximated using methods such as Monte Carlo approximations and others. The Lyapunov exponents of a sequence of random matrices was investigated in the pioneering paper of Furstenberg and Kesten [8] and by Oseledec in [18]. Ruelle [21] developed the theory of Lyapunov exponents for random compact linear operators acting on a Hilbert space. Newman in [15] and [16] and later Isopi and Newman in [12] studied Lyapunov exponents for random N × N matrices as N → ∞. Later on, Vladislav Kargin [14] investigated how the concept of Lyapunov exponents can be extended to free linear operators (see [14] for a more detailed exposition). The probability distribution of the Lyapunov exponents associated to the sequence {Y p } ∞ p=1 , is the spectral probability distribution γ of the Hermitian operator with spectral measure given by γ := ln(ν). Moreover, γ is absolutely continuous with respect to Lebesgue measure and from our work we can obtain a closed-form expression for its Radon-Nikodym derivative with respect to Lebesgue measure.
This problem is not only interesting from the mathematical point of view but it is also important for the information theory community. For example, it was studied in [5] that if one is interested in exploring the performance in a layered relay network having a single source destination pair, then the message is passed from one relay layer to the next till it reaches the destination. Assume that there are n + 1 layers of relay nodes between the source and the destination, with each layer having N relay nodes. Each relay node has a single antenna which can transmit and receive simultaneously. Thus, the complete state of this network is fully characterized by the n channel matrices denoting the channels between adjacent layers. The matrix Y m denotes the channel between layer m and m + 1, i.e., Y m (i, j) is the value of the channel gain between the i-th node in layer L m+1 and the j-th relay node in layer L m . Thus, the amplify-and-forward scheme converts the network into a point-to-point MIMO system, where the effective channel is Y n Y n−1 . . . Y 1 the matrix product of Gaussian random matrices. Under the appropriate hypothesis the capacity of this channel is given by where snr is the signal to noise ratio of the system. Now we will describe the content of this paper. In Section §2, we recall some necessary preliminaries as well as some known results. In Section §3, we prove our main Theorem and present some examples and simulations. In Section §4, we derive the probability distribution of the Lyapunov exponents of the sequence {Y p } ∞ p=1 . Finally, in Section §5 we provide the proofs.

Preliminaries and Notation
We begin with an analytic method for the calculation of multiplicative free convolution discovered by Voiculescu. Denote C the complex plane and set for z ∈ C\[0, ∞). The measure ν is completely determined by ψ ν . The function ψ ν is univalent in the half-plane iC + , and ψ ν (iC + ) is a region contained in the circle with center at −1/2 and radius 1/2. Moreover, , the function ψ ν has an inverse with respect to composition Finally, define the S-transform of ν to be See [3] for a more detailed exposition.
Denote by the family of all compactly supported probability measures defined in the real line R. We denote by + the set of all measures in which are supported on [0, ∞). On the set there are defined two associative composition laws denoted by * and . The measure µ * ν is the classical convolution of µ and ν. In probabilistic terms, µ * ν is the probability distribution of a + b, where a and b are commuting independent random variables with distributions µ and ν, respectively. The measure µ ν is the free additive convolution of µ and ν introduced by Voiculescu [24]. Thus, µ ν is the probability distribution of a + b, where a and b are free random variables with distribution µ and ν, respectively. There is a free analogue of multiplicative convolution also. More precisely, if µ and ν are measures in + we can define µ ν the multiplicative free convolution by the probability distribution of a 1/2 ba 1/2 , where a and b are free random variables with distribution µ and ν, respectively. The following is a classical Theorem originally proved by Voiculescu and generalized by Bercovici and Voiculescu in [4] for measures with unbounded support.
for every z in the connected component of the common domain of S µ and S ν .
Analogously, the R-transform is an integral transform of probability measures on R. Its main property is that it linearises the additive free convolution. More precisely, the following result proved by Voiculescu [24] holds.
for every z in the connected component of the common domain of R µ and R ν .
Hence in the analogy between the free convolution and the classical one, the R-transform plays the role of the log-Laplace transform.

Main Results
In this Section we prove our main results. Let us first fix some notation. We say that two N × N random matrices A and B have the same * -distribution if and only if The proof of this result is in Section 5.1.
Since the random matrices {X p } ∞ p=1 are Gaussian and independent then the random matrices {Y p } ∞ p=1 are asymptotically free as N → ∞ (see [23,24,17]). Denote by µ the limit distribution measure of the sequence Y * 1 Y 1 as the dimension N → ∞. More precisely, µ is the compactly supported probability measure such that for every k ≥ 1. The measure µ depends only on σ and their relation is well known. Moreover, this topic is a focus of a lot of work in the information theory community, since it gives information on the capacity of the correlated Gaussian channel (see [23,24,17,1,2,20,19,22] for more on this). If we consider the random matrix B 2 defined as is not difficult to see that its limit measure is µ µ the multiplicative free convolution of µ with itself. This is essentially because the moments of B 2 are the same as the moments of Y * 2 Y 2 Y * 1 Y 1 by the trace property. Analogously, the random matrices B n converge in distribution to as N → ∞. The relationship between the measure µ and µ n is given by Voiculescu's S-transform as explained in the previous Section. Our interest is in the normalized version of µ n . More specifically, in the measure ν n defined as the limit distribution of B 1/2n n as N → ∞. The relationship between the moments of µ n and ν n is given by for every k ≥ 1.
Since the measure µ n is compactly supported then it is clear that ν n is compactly supported as well. Our interest is in the asymptotic behavior of these measure as n → ∞. Now we are ready to state our main Theorem.

Theorem 3.2.
Let {Y k } k be a sequence of random matrices as before. Let µ in + and B n be as before. The sequence of measures ν n converges in distribution to a compactly supported measure ν.

Examples
In this Section we present an example and some simulations.
Example 3.4. The simplest case is when Σ(N ) = I N is the identity N × N matrix. Therefore, the measure σ = δ 1 and µ is the limit spectral measure of X * 1 X 1 which is known to be the Marchenko-Pastur distribution with parameter one. Its density is given by A simple computation shows that the S transform S µ (z) = 1 z + 1 .

Hence, by Theorem 3.2 we see that
In Figure 1 we see the spectral measure of X * 1 X 1 for N = 300 whose limit as N → ∞ is the well known Marchenko-Pastur distribution of parameter 1. In Figure 1 we also see the spectral measure of B 1/2 1 = (X * 1 X 1 ) 1/2 for N = 300 (whose limit is the well known quarter-circular law). Analogously, in Figure 2 we see the spectral measure of B 1/4 2 = (X * 1 X * 2 X 2 X 1 ) 1/4 for N = 300. Finally, in Figure 2 we also show the spectral measure of B 1/12 6 for N = 500 also. We can appreciate that as n increases the spectral measures of the operators converge to the ramp measure described in the previous example. Further simulations show that this convergence is relatively slow.

Lyapunov Exponents of Random Matrices
{Y k } ∞ k=1 be the sequence of random matrices as before. Let µ be the spectral probability measure of Y * 1 Y 1 and assume that µ({0}) = 0. Using Theorem 3.2 we know that for every fixed n the sequence of random matrices converges in distribution to a compactly supported measure in [0, ∞) as the dimension of the matrices N → ∞. Moreover, the sequence of measures ν n converge weakly to a compactly supported measure This distribution is absolutely continuous with respect to the Lebesgue measure and has Radon-Nikodym derivative (t) and F µ (t) = S µ (t − 1) −1/2 . Let Λ be a random variable with probability distribution ν and let L be the possibly unbounded random variable defined by L := ln(Λ), and let  γ be the spectral probability distribution of L. It is a direct calculation to see that γ is absolutely continuous with respect to Lebesgue measure and has Radon-Nikodym derivative The probability distribution γ of L is what is called the distribution of the Lyapunov exponents (see [15], [16] and [21] and [14] for a more detailed exposition on Lyapunov exponents in the classical and non-classical case). Theorem 4.1. Let {Y k } k be a sequence of random matrices as before. Let µ in + and B n be as before. Let γ be probability distribution of the Lyapunov exponents associated to this sequence. Then γ is absolutely continuous with respect to Lebesgue measure and has Radon-Nikodym derivative

Remark 4.2.
Note that if the operator Y * 1 Y 1 is not invertible in the · 2 then the random variable L is unbounded.
The following is an example done previously in [14] using different techniques.
Therefore, we see that the probability measure of the Lyapunov exponents is γ with This law is the exponential law discovered by C. Newman as a scaling limit of Lyapunov exponents of large random matrices. (See [15], [16] and [12]). This law is often called the "triangle" law since it implies that the exponentials of Lyapunov exponents converge to the law whose density is in the form of a triangle.

Proof of Lemma 3.1
Proof. Let Y k = U k A k be the polar decomposition of the matrix Y k , where A k is positive definite and U k is a unitary matrix. We will proceed by induction on n. The case n = 1 is obvious since Y * 1 Y 1 = A 2 1 . Assume now that B k has the same distribution as b k for k < n. Then by the unitary invariance and the induction hypothesis Hence Since conjugating by a unitary does not alter the distribution we see that Since the random matrices {Y k } ∞ k=1 are independent then {{U k , A k }} ∞ k is also an independent family and A 1 ∼ d U 1 A 1 U * 1 and independent with respect to {A k } k≥2 . Then, concluding the proof.

Proof of Theorem 3.2
Before staring the proof let us review some necessary results.
In [17] Nica and Speicher introduced the class of R-diagonal operators in a non commutative C * -probability space. An operator T is R-diagonal if T has the same * -distribution as a product uh where u and h are * -free, u is a Haar unitary, and h is positive. The next Theorem and Corollary were proved by Uffe Haagerup and Flemming Larsen ( [10], Theorem 4.4 and the Corollary following it) where they completely characterized the Brown measure of an R-diagonal element. We will state their Theorem for completeness.

µ T is the only rotation symmetric probability measure satisfying (3).
Corollary 5.2. With the notation as in the last Theorem we have has an analytic continuation to a neighborhood of its domain and F > 0 on (µ h ({0}), 1).

µ T has a radial density function f on
Proof of Theorem 3.2: From the previous Lemma it is enough to prove the Theorem for A k = |Y k |. The sequence of random matrices {A k } ∞ k=1 converge in distribution to a sequence of free and identically distributed operators {a k } ∞ k=1 as the dimension N → ∞. Therefore, the measure ν n can we characterized as the spectral measure of the positive operator from h it is known (see [23], [24]) that the family of operators {u k h(u * ) k } ∞ k=0 is free. Therefore, defining c k = u k h(u * ) k we see that T * T ∼ d c 2 1 , (T * ) 2 T 2 ∼ d c 2 c 2 1 c 2 and it can be shown by induction that (T * ) n T n ∼ d c n c n−1 · · · c 2 1 · · · c n−1 c n . Therefore, since c k has the same distribution than a k , and both families are free, we conclude that the operators (T * ) n T n and b n have the same distribution. Moreover, by Theorem 2.2 in [11] the sequence (T * ) n T n 1 2n converges in distribution to a positive operator Λ. Let ν be the probability measure distribution of Λ. If the distribution of a 2 k is a Dirac delta, µ = δ λ , then h = λ and (T * ) n T n 1 2n = λ n (u * ) n u n 1 2n = λ.
Therefore, b 1 2n n has the Dirac delta distribution distribution δ λ and ν = δ λ . If the distribution of a k is not a Dirac delta, let µ T the Brown measure of the operator T . By Theorem 2.5 in [11] we know that We know by Theorem 5.1 and Corollary 5.2 that µ T = βδ 0 + ρ with dρ(r, θ ) = 1 2π f (r) 1 (F µ (β),F µ (1)] (r) d r dθ (15) where f (t) = F <−1> µ (t) and F µ (t) = S µ (t − 1) −1/2 . Hence, using equation (14) we see that r p f (r)d r for all p ≥ 1. Using the fact that if two compactly supported probability measures in + have the same moments then they are equal, we see that By Corollary 5.2, we know that