Universal distribution of Lyapunov exponents for products of Ginibre matrices
Starting from exact analytical results on singular values and complex eigenvalues of products of independent Gaussian complex random matrices also called Ginibre ensemble we rederive the Lyapunov exponents for an infinite product. We show that for a large number of product matrices the distribution of each Lyapunov exponent is normal and compute its -dependent variance as well as corrections in a large- expansion. Originally Lyapunov exponents are defined for the singular values of the product matrix that represents a linear time evolution. Surprisingly a similar construction for the moduli of the complex eigenvalues yields the very same exponents and normal distributions to leading order. We discuss a general mechanism for matrices why the singular values and the radii of complex eigenvalues collapse onto the same value in the large- limit. Thereby we rederive Newman’s triangular law which has a simple interpretation as the radial density of complex eigenvalues in the circular law and study the commutativity of the two limits and on the global and the local scale. As a mathematical byproduct we show that a particular asymptotic expansion of a Meijer G-function with large index leads to a Gaussian.
Lyapunov exponents are useful to study the stability of dynamical systems, but they also play an important role in statistical mechanics of disordered systems, localization theory, hidden Markov models and many others areas of physics and engineering.
The problem of the determination of Lyapunov exponents is intimately related to the asymptotic properties of products of random matrices in the limit when the number of factors tends to infinity. The randomness encoded in these matrices depends on the details of the problem in question and it is usually very difficult to find the exact values of the exponents. There are however some general theorems that guide the calculations. For example it is known that the largest Lyapunov exponent of the product of a random sequence of matrices generated by a stochastic process converges almost surely to a limiting deterministic value in the limit of infinite sequence length. For large but finite sequences the largest Lyapunov exponent is a normally distributed random variable with the variance inversely proportional to the sequence length .
The relevance of products of random matrices to dynamical systems and ergodic theory was realized in the sixties  and since then the study of matrix products has been an active field of research in probability theory , condensed matter physics, and statistical mechanics [4, 5, 6].
It was noticed long time ago [7, 8] that products of random matrices naturally arise in the analysis of disordered systems in statistical mechanics. As an example one can think of the transfer matrix formulation of random Ising chains [9, 10]. In this case the transfer matrices are random matrices. In the thermodynamic limit the free energy density is given by the largest Lyapunov exponent of the product of transfer matrices. Another important physical example is the localization phenomenon in electronic systems . In this case the leading Lyapunov exponent is related to the inverse localization length [12, 13, 14]. Other solvable physical models can be found in Yang-Mills theories . In this field unitary transfer matrices in the group find applications in calculations of the Wilson loop operator for .
Products of random matrices have many practical applications in other fields as well. For instance they arise in calculations of the capacity of a sequence of multiple-input-multiple-output arrays in wireless telecommunication [17, 18, 19] and in hidden Markov models applied in stochastic inference , in time series analysis, speech recognition, biological sequence analysis. In hidden Markov models the Lyapunov exponents correspond to the entropy rates [21, 22]. Also in image processing  product matrices play an important role.
The spectrum of Lyapunov exponents gives important information on the stability and the complexity of dynamical systems  and their effective information dimension . For this reason a great effort has been made to develop computational methods to determine Lyapunov exponents for given theoretical models or to estimate them from experimental data. Numerical methods are directly based on the analysis of the equation of motion or measurements of the expansion rates of phase space [25, 26]. Algorithms have been developed for the Lyapunov spectrum from sample time series . Also analytical approximations include methods based on the weak disorder expansion  or properties of determinants associated with transfer matrices [29, 30].
There are only a few models where Lyapunov exponents can be calculated exactly. They usually involve products of matrices with randomness controlled by a single random parameter where the exact expressions result from some model specific simplifications which occur during calculations. The examples include classical disordered harmonic chains [7, 31], the tight-binding Anderson model [32, 33], quantum spin chains [34, 35, 36] and random Schrödinger operators , see also [4, 5, 6] for reviews. Recently a general method has been worked out to derive a scaling form for the Lyapunov exponents in the continuum limit for products of matrices close to the identity  based on the Iwasawa decomposition of SL(2,R) .
An important solvable case where one can calculate the Lyapunov exponents exactly is the product of identically distributed Gaussian random matrices with independent, identically distributed (i.i.d.) centered real entries . Such matrices are usually called real Ginibre matrices. This is a special case, first of all because one can analytically derive the whole spectrum of Lyapunov exponents for any system size . Second, the calculation uncovers a deep connection between the spectrum and the law of large numbers . The exponents are exclusively shaped by the statistics of matrix elements and not by the matrix structure. In other words the two effects do not mix. A second much more recent example where all Lyapunov exponents have been calculated are products of independent Ginibre matrices, where each factor is multiplied by a fixed positive definite matrix [41, 42]. When these constant matrices are equal to the identity the results for the real, complex, and quaternion Ginibre ensembles agree up to a scaling factor where is the Dyson index.
The fact that one can derive the whole spectrum is very useful for practical purposes since the spectrum can be used to test numerical algorithms [25, 26, 27]. Moreover one can analytically calculate the limiting law for the distribution of Lyapunov exponents in the limit . For the numbers constructed from Lyapunov exponents, that we call in this paper incremental singular values, , , the distribution is given by the triangular law .
In the present work we further elaborate on the Lyapunov spectrum for the product of complex Ginibre matrices. We consider complex Ginibre matrices that are Gaussian matrices with i.i.d. complex elements. We derive an exact form of finite corrections to the Lyapunov spectrum, where is the number of matrices in the product. For finite the Lyapunov exponents are random variables. We calculate the joint probability distribution for these variables. For large it is asymptotically given by a permanent of the product of independent Gaussian functions centered at the limiting values. Thereby we determine the widths of the distributions. We also improve this Gaussian approximation by considering another approximation based on the saddle point approximation. The latter approach works even better for a product of a small number of matrices since it still incorporates asymmetric parts of the individual eigenvalue distributions and to a small extent the original level repulsion.
In addition to the Lyapunov exponents , which are related to the singular values of the product matrix, one can define the corresponding exponents for the moduli of the complex eigenvalues. The complex eigenvalue distribution of the product of Ginibre matrices is rotationally invariant in the complex plane [43, 44]. We find that the moduli of the eigenvalues become uncorrelated random variables in the large- limit and we determine the form of their joint probability distribution. Surprisingly, the spectrum and the joint probability distribution of these exponents is identical to that of the Lyapunov exponents, for .
A further consequence of this observation is discussed in Section 5. The triangular law for Lyapunov exponents corresponding to the singular values found by Isopi and Newman  can be understood as the radial distribution of eigenvalues of the Ginibre matrix. The fundamental reason behind this interpretation is twofold. First, our insight says that the Lyapunov exponents constructed from the singular values and from the moduli of the eigenvalues agree with each other. Second, Ginibre ensembles belong to the class of isotropic random matrix ensembles. For those ensembles the sometimes called self-averaging property of the product of isotropic matrices [46, 47] and the Haagerup-Larsen theorem  are known. These two properties imply that the spectral statistics of a product of independent random matrices is equal to the statistics of the power of a single matrix in the limit of large matrix dimension . After taking the root of the product matrix the level density is the one of an ordinary Ginibre matrix which is the circular law for the complex eigenvalues and is equal to the triangular law for the moduli of the eigenvalues.
The paper is organized as follows. In Section 2 we define the linear evolution given by the product of Ginibre matrices and define the corresponding Lyapunov exponents. In Section 3 we derive their joint probability density based on the singular value distribution of the product matrix for finite and large , keeping finite. In Section 4 we compute the joint probability density for exponents based on the moduli of complex eigenvalues for finite and large . In Section 5 we discuss the limit for Lyapunov exponents and show that this limit commutes with the limit on the global scale while it does not commute on the local scale of the mean level spacing. In Section 6 we conjecture the collapse of singular and eigenvalues for general isotropic ensembles and exemplify this for . We conclude the paper in Section 7. In the appendices we recall some identities of Meijer G-functions, compute a particular kind of a Hankel determinant and present some further details of our calculations.
2 Linear time evolution with Ginibre matrices
Let us consider a linear discrete-time evolution of an -dimensional system described by complex degrees of freedom. The state of the system at time is given by an -dimensional vector . The state at is related to the state at time by the following linear equation
with the evolution operator represented by an matrix. The total evolution from the initial state
is effectively driven by the product matrix
Here we are interested in ’s being i.i.d. complex non-Hermitian random matrices. In particular we consider the case of Ginibre matrices which centered and Gaussian distributed,
for all . The differential denotes the product of the differential of all independent matrix elements. Towards the end of the paper we comment on the evolution for general isotropic random matrices which are defined by the invariance of the probability measure where are arbitrary unitary matrices. Isotropic matrices are sometimes called bi-unitarily invariant or rotational invariant. Ginibre matrices belong to this class.
We are interested in the large behavior of the system, approximating a continuous time evolution. This behavior is controlled by the Lyapunov exponents which are related to the singular values of . Let us denote the real eigenvalues of the positive matrix
by . Their square roots correspond to singular values of . Then the Lyapunov exponents are defined as
where are the ordered eigenvalues of : . Throughout this paper we denote ordered (increasing) sequences like or by a hat.
In many physical situations the number of time steps in the evolution is large but finite. Hence it is interesting to study finite size corrections to the limiting values, and the rate of convergence to these values. Thus we want to address the question how this limit is realized when tends to infinity (). Our focus lies on the corresponding quantities for finite
which we call finite Lyapunov exponents, for . In the limit , after ordering, they become the standard Lyapunov exponents: . We look for a probabilistic law that governs the distribution of the finite Lyapunov exponents, or equivalently their joint probability density for finite and . Given the recent progress on the joint distribution of singular values (and complex eigenvalues) for a finite product of Ginibre matrices for finite and this can be easily calculated, and the limits and subsequently can be taken.
3 Lyapunov exponents from singular values
where is the Vandermonde determinant
The function is a particular case of the Meijer G-function (1.1) whose properties and definition are recalled in A. As any special function, it possesses many helpful properties which facilitate calculations. For simplicity we drop the explicit -dependence of the singular values and of the Lyapunov exponents in the ensuing discussions as it will be clear from the context if is finite or infinite.
The road map to find the large asymptotics is the following. In subsection 3.1 we find a determinantal representation of the joint probability distribution of Lyapunov exponents made of one-point probability distributions. We calculate the moments of these one-point distributions. The cumulant expansion yields an asymptotic expansion to any order in . This result is discussed in detail for large , in subsection 3.2. Moreover, we compare the cumulant expansion with a saddle point approximation which also incorporates a residual level repulsion as well as an asymmetric part of the distributions of the individual Lyapunov exponents. In subsection 3.3 we come back to the discussion of the corresponding singular values which we call incremental singular values since they are the average contribution to the total singular value of each single random matrix in the product .
3.1 Reduction to “decoupled” random variables
The joint probability distribution for Lyapunov exponents can be directly read off from eq. (LABEL:Ps) by the change of variables ,
The change of variables introduces a Jacobian which yields for each variable the exponential factor . These factors have been absorbed in the last equation in the Vandermonde determinant by replacing . The first determinant in eq. (3.1) can be expanded as
where denotes the group of permutations of elements and “” is the sign function which is for even permutations and for odd ones. The factors can be absorbed into the second determinant
By virtue of eq. (1.3) the last expression can be cast into the form
The skew-symmetry of the determinant under permutations of its rows and columns allows us to absorb the prefactor into the determinant via rearranging the rows. Hence we end up with
Thus the problem is reduced to the analysis of the function . By construction this function is positive semi-definite. With help of the integral identity (1.2), can be normalized such that the function
can be interpreted as a probability density for a single random variable. Replacing with its normalized version the joint probability distribution reads
In passing from eq. (3.7) to eq. (3.10) we have pulled the factor out of the determinant. This factor cancels the corresponding prefactor in eq. (3.7) leaving the product of the second powers in front of the determinant in eq. (3.10).
Using the cumulant expansion we argue in the next subsection that the probability densities can be approximated by Gaussian functions in the limit . Therefore let us define the moment generating function
where are the moments. This moment generating function can be calculated with help of eq. (1.2),
The expansion in at yields the moments . The logarithm of the moment generating function is the cumulant generating function
The coefficients of the corresponding Taylor series of at are the cumulants ,
Hereby we employed the definition of the digamma function and its derivatives,
The first cumulant (=first moment) corresponds to the mean value and is equal to
The second cumulant corresponds to the variance and takes the value
We emphasize that so far all results are exact for finite .
3.2 Large limit
We apply the standard argument based on the analysis of the large- behavior of cumulants to show that can be approximated by a Gaussian function for large . Thereby we have first to center the distribution and normalize its second moment. The exact limit will yield a Gaussian. This limit justifies to replace by a Gaussian centered at and with the standard deviation .
For this purpose we define the standardized random variable . Thereby we denote standardized quantities by in this and the next section. The random variable is distributed as . The same notation is applied for cumulants. By construction, the standardized mean is and the standardized variance is . The higher standardized cumulants are
They tend to zero when goes to infinity. Therefore the standardized cumulant generating function is in the limit ,
By analytic continuation to imaginary values we get and hence . The inverse Fourier transform for the moment generating function yields the limit
Inverting the process of standardization we get the following asymptotic expansion
with and given by eqs. (3.16) and (3.17). In other words, for large we can replace in (3.10) by the Gaussian function eq. (3.21). Here we have also reinserted the definition of from eqs. (3.9) and (3.8) in order to stress that this is the first main result of this section, namely the asymptotic expansion of a Meijer G-function in the double scaling limit of large argument and large index. We are not aware of such a result in the literature. In particular it is different from the well-known large argument expansion, cf. .
and hence with
Here the sum over permutations without signs is equal to the definition of the permanent, . The prefactor simplifies to since
as recalled in B.
Let us state the main result of this section in its explicit form which is the joint probability distribution for large ,
The limiting joint probability distribution sustains its invariance under permutations of the indices, . More explicitly, the joint probability density is a symmetrized product of one-point functions or densities, which means in physical language that it describes a system of independent, non-interacting, indistinguishable bosons. Starting from the determinantal process of the singular values the appearance of a permanent is somewhat surprising, whereas it quite naturally arises for complex eigenvalues after integrating over the angles, see e.g. in [50, 51]. We will come back to this point at the end of section 4.
Note that the dependence of on appears only through the widths of the Gaussian peaks. Their positions are independent of in this approximation.
The density defined as
is in our case
When increases the peaks become more narrow and, eventually in the limit , the Gaussian peaks turn into Dirac delta functions and we recover the deterministic laws [40, 41] for the Lyapunov exponents ,
Employing Newman’s argument  one can show that the positions of the peaks for general Dyson index are given by with . Thus the positions we calculated fit into the results obtained for products of real Ginibre matrices by Newman  and agree with the more general recent result by Forrester  who considered complex Ginibre matrices multiplied by a fixed positive definitive matrix. Forrester’s work was extended by Kargin  to . Let us emphasize that our result (3.28) gives finite- corrections to this deterministic law. Moreover we stress that the same limit has a corresponding consequence for the Meijer G-functions for the individual peaks, when taking the limit ,
Already for finite but sufficiently large when the peaks cease to overlap, each Gaussian peak , see eq. (3.23), can be identified as a finite size distribution of the -th largest Lyapunov exponent . Due to the recursion the distance between neighboring peaks is and the sum of their widths is . So the peaks separate when implying . Thus, for the system with degrees of freedom all peaks get separated for . Note that the positions and the widths are independent of . When increases, just new peaks appear in the distribution while the old ones neither change in shape nor shift their positions.
Let us study the quality of the approximation that has led us to eq. (3.26). In the derivation of the asymptotic form (3.26) for large we used the fact that the functions can be approximated by Gaussian functions (3.21) and that their mean values and their variances asymptotically only depend on a single index , see eq. (3.22), if one neglects terms. The 1/t terms have a twofold effect on the shape of the density. First, the positions and widths of the peaks solely resulting from the single random variable distributions weakly dependent on . Second a repulsion between peaks is introduced due to the determinant in eq. (3.7). We illustrate these two effects in Fig. 1 for the level density where we compare the asymptotic formula (3.28) and a saddle point approximation of for the inverse Fourier transform of the moment generating function (3.11),
and being the Heaviside function. This approximation is derived in C. Note that in the large limit the distribution indeed becomes the Gaussian (3.23) and independent of the index . The level density in the approximation (3.32) is
Hereby we integrated over all but one Lyapunov exponents, , and we expanded the determinant (3.10) in the columns and rows where the remaining distribution stands. Note that as well as are normalized. The cofactor of the Hankel determinant (3.25) is calculated in B.
The main conclusion from the comparison in Fig. 1 is that the corrections do not have any significant effect on the shape of the distribution when the peaks are separated. In particular for the smallest singular values this requirement is often satisfied. Nevertheless the corrections can become quite important for up to in which case the saddle point approximation (3.32) is better suited. For the largest eigenvalues the effect of these corrections is the strongest.
3.3 Incremental singular values
We close this section by going back to the singular values because in some physical situations it is more convenient to use them rather than Lyapunov exponents. Consider the -th root of the matrix ,
in contrast to eq. (2.5). We define incremental singular values as
which correspond to the real positive eigenvalues of the matrix . Intuitively, the incremental singular values give the typical incremental contraction or expansion factors for the configuration space under a single average time step of the evolution. Of course they contain exactly the same information as the Lyapunov exponents. Their joint probability distribution is obtained from that for the Lyapunov exponents by the simple change of variables in eq. (3.36) inserted in eq. (3.1). Using eq. (3.10) this gives
For large when is approximated by normal distributions, can be approximated by log-normal distributions. Otherwise everything works exactly in the same way as for Lyapunov exponents. In particular, when is large enough to neglect the corrections, we obtain the counterpart of eq. (3.26)
and , are given by eq. (3.22). The functions have maxima at . The density of incremental singular values is given by the normalized sum
in analogy to eq. (3.28). Again this turns into a sum of delta functions in the limit ,
We have tested this prediction against Monte-Carlo simulations for finite size systems. In Fig. 2 we show histograms of incremental singular values calculated analytically and numerically. We see that the log-normal functions provide a very good approximation to the actual shapes.
4 Lyapunov exponents from the moduli of complex eigenvalues
Rather than using singular values, the complex eigenvalues, , , are an alternative way to characterize the spectral properties of the matrix , see eq. (2.3). In general the singular values and the moduli of complex eigenvalues are unrelated, apart from their product which is equal to and bounds on their Euclidean norm which result from the trace (see eq. (6.23)), respectively. However in the large limit, the moduli of the complex eigenvalues will behave exactly in the same way as the singular values . In fact repeating the same construction as in section 3, taking the -th root of will lead to the very same normal distribution, frozen at identical positions as the limiting singular values. For that reason we will use the same term Lyapunov exponent which is otherwise reserved for the singular values, only.
We pursue a calculation similar to section 3. Thereby we first show that all complex eigenvalues can be traced back to decoupled random variables apart from a trivial determinantal coupling, see subsection 4.1. In the second step we employ the cumulant expansion to find Dirac delta functions in the leading order and Gaussian (for the corresponding Lyapunov exponents) and log-normal (for the moduli of eigenvalues) distributions in the next-to-leading order, see subsection 4.2. In subsection 4.3, we present an alternative approach by first integrating over the angles and then taking the limit . This alternative construction is also applied to the case since the analytical result for the joint probability density of the complex eigenvalues is known [53, 54, 51] for this case as well.
4.1 Reduction to “decoupled” random variables
The definition (2.6) of Lyapunov exponents requires to take the -th root and the logarithm of the positive singular values. However, for complex variables this is not a unique procedure. If one takes for example the root , the question arises which of the roots we have to take. When choosing the primary root the resulting spectrum will be mapped onto a circular sector of the angle which eventually shrinks to the positive semi-axis in the limit . Another alternative choice is to take the root of the moduli of the eigenvalues only, i.e.
Indeed this choice seems to be a more natural construction. When multiplying the product by new matrices, the angular parts of the eigenvalues will run around on circles while the radial part will either exponentially contract or expand. Thus it is not the angular part we have to worry about in the large limit since it stays in a compact set. It is the radial part of the eigenvalues which has to be rescaled such that the support stays fixed. Therefore we decide for the rooting (4.1). We emphasize that the kind of rooting is crucial to find our results which may change for other constructions.
The definition of the Lyapunov exponents at finite and infinite starting from the moduli of complex eigenvalues are
These definitions are the analog of those for the Lyapunov exponents corresponding to the singular values, see Eqs. (2.6) and (2.7). Hereby recall that the variables are the squared singular values which results in an additional prefactor .
where is the flat measure in the complex plane. As in the previous section we again drop the explicit -dependence of all quantities. We change to polar coordinates and employ the variables (4.2) such that the joint-probability distribution reads
We extend the first product by the identity . With help of the identity we get