Real eigenvalues of non-Gaussian random matrices and their products

# Real eigenvalues of non-Gaussian random matrices and their products

Sajna Hameed 1 School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455. U.S.A..    Kavita Jain 2 Theoretical Sciences Unit, Jawaharlal Nehru Centre for Advanced Scientific Research,
Jakkur P.O., Bangalore 560064, India.
Arul Lakshminarayan 3 Department of Physics, Indian Institute of Technology Madras, Chennai, 600036, India Max Planck Institute for the Physics of Complex Systems, 01187 Dresden, Germany
11e-mail:hamee007@ umn.edu
22e-mail: jain@jncasr.ac.in
33e-mail: arul@physics.iitm.ac.in
September 23, 2019
###### Abstract

We study the properties of the eigenvalues of real random matrices and their products. It is known that when the matrix elements are Gaussian-distributed independent random variables, the fraction of real eigenvalues tends to unity as the number of matrices in the product increases. Here we present numerical evidence that this phenomenon is robust with respect to the probability distribution of matrix elements, and is therefore a general property that merits detailed investigation. Since the elements of the product matrix are no longer distributed as those of the single matrix nor they remain independent random variables, we study the role of these two factors in detail. We study numerically the properties of the Hadamard (or Schur) product of matrices and also the product of matrices whose entries are independent but have the same marginal distribution as that of normal products of matrices, and find that under repeated multiplication, the probability of all eigenvalues to be real increases in both cases, but saturates to a constant below unity showing that the correlations amongst the matrix elements are responsible for the approach to one. To investigate the role of the non-normal nature of the probability distributions, we present a thorough analytical treatment of the single matrix for several standard distributions. Within the class of smooth distributions with zero mean and finite variance, our results indicate that the Gaussian distribution has the maximum probability of real eigenvalues, but the Cauchy distribution characterised by infinite variance is found to have a larger probability of real eigenvalues than the normal. We also find that for the two-dimensional single matrices, the probability of real eigenvalues lies in the range .

## I Introduction

The problem of the number of real roots of algebraic equations has a long history, but continues to attract attention, from the early works of Littlewood and Offord Littlewood and Offord (1939) to more recent developments Edelman and Kostlan (1995); Nguyen et al. (2014) (see the latter article for more related history and references). Mark Kac, in a seminal paper, proved that the expected number of real zeros of polynomials (of order ) whose coefficients are chosen from a normal distribution of zero mean is Kac (1943). Others substantially extended these results and showed universality in that this is the leading order behavior, irrespective of the underlying distribution as long as they are zero-centered and have a finite variance Stevens (1969). Logan and Shepp Logan and Shepp (1966) studied the same when the coefficients are Cauchy distributed (and hence have infinite variance, and undefined average) the leading behavior goes as , with which is larger than , and hence there are more real zeros in the Cauchy case than when the variance is finite. Thus, the number of real roots of a random polynomial are generally quite small.

More recently, several works have explored the fraction of real eigenvalues for matrices drawn from the real Ginibre ensemble. This ensemble consists of matrices whose elements are independently drawn from a normal distribution such as . For such matrices, it has been shown analytically that the expected number of real eigenvalues and the probability that all eigenvalues are real are given by Edelman et al. (1994); Edelman (1997); Kanzieper and Akemann (2005)

 limn→∞En√n=√2π,Pn,n=2−n(n−1)/4. (1)

Thus the probability that all eigenvalues are real tends to zero as the matrix dimension increases, although the expected number of them increases algebraically.

Products of random matrices have also been studied, at least, since the classic work of Furstenberg and Kesten Furstenberg and Kesten (1960) (see, Bougerol and Lacroix (1985); Crisanti et al. (1993) for a discussion of subsequent work and some applications). The study of the spectra and singular values of the products of random matrices is currently a very active area of research Roga et al. (2011); Burda et al. (2012), and recent progress has been reviewed in Akemann and Ipsen (2015). When studying a problem related to the measure of “optimally entangled” states of two qubits Shuddhodan et al. (2011), the problem of the number of real eigenvalues of a product of matrices came up. It was shown in Lakshminarayan (2013) that the probability of all eigenvalues being real for the product of two, matrices from the real Ginibre ensemble is . Somewhat surprisingly then there is a lesser probability that all eigenvalues are real for a single matrix (, from Eq. (1)) than for a product of two. This observation motivated numerical explorations in Lakshminarayan (2013) which showed that this probability tends to as the number of matrices in the product increases. Numerical results showing that this is true for higher dimensions was also discussed therein. This result has been proven analytically by Forrester Forrester (2014), who calculated explicitly in terms of Meijer G-functions for square matrices, and generalised by Ipsen Ipsen (2015) to rectangular matrices.

Most of the previous works on the properties of the eigenvalues have studied Gaussian matrices with independent, identically distributed (i.i.d.) entries. Exceptionally, Edelman, Kostlan and Shub Edelman et al. (1994) presented some numerical results indicating universality for the expected value of real eigenvalues of real matrices. However, to the best of our knowledge, there is no study exploring such extensions to the products of matrices. Thus, two conditions are relaxed: first, non-Gaussian i.i.d. elements are taken as entries of the matrices, and second, we consider products of such independent matrices. In this case, the entries of the product matrix are naturally correlated in a complex manner and their distribution is, in general, anyway non-Gaussian.

Motivations for studying products of random matrices are well known in physics, and range from the study of localization in random media where transfer matrices are multiplied to dynamical systems where products of local stability matrices determine the Lyapunov exponents Crisanti et al. (1993). Applications for counting number of real eigenvalues of random matrices, especially of products, is of more recent provenance. The real eigenvalues of a class of random matrices indicate topologically protected level crossings at the Fermi energy in the bound states of a Josephson junction Beenakker et al. (2013). While the relevant object is a single matrix and not a product, it is conceivable that this provides a context for further investigation. In the context of entanglement of two qubits Lakshminarayan (2013), the real eigenvalues of products of two matrices were relevant in finding the measure of optimally entangled states Shuddhodan et al. (2011). Different distributions of matrix elements would result in different sampling of states on the Hilbert space. Finally, it is worth pointing out that the fact that real eigenvalues dominate the spectrum of products of matrices implies that generic orbits of dynamical systems are not of the complex unstable variety. Complex instability can occur in Hamiltonian systems with more than two degrees of freedom Contopoulos et al. (1994), when exponential instability is combined with rotational action in phase space. However, as a consequence of discussions in this paper, it seems unlikely that they would generically arise for long orbits, as they depend on eigenvalues of products of Jacobian matrices being complex.

To begin with, we ask: how do non-normal distributions affect the probability of real eigenvalues of a single matrix? In fact, a major part of this paper is concerned simply with matrices and various families of the distributions of matrix elements, and the probability that the eigenvalues are real is calculated exactly in many cases in Sec. II. In particular, we find the lower and upper bounds on the probability that both eigenvalues are real. In the family of symmetrized Gamma distributions, the range of probabilities is shown to be . The upper bound is associated with the case when the matrix elements have a large weight at the origin and the lower bound in the opposite case when the maximum is away from the origin. In fact, these same bounds are obtained in a class of truncated distributions thus indicating that the probability of real eigenvalues depends crucially on whether the distributions are concentrated near the origin or away from it. The lower bound of is shown to be a tight one by proving this in general. The upper bound of still remains specific to these distributions, however we do believe that this is widely applicable as well.

We also find that the other features of the distributions of the matrix elements that play a significant role in determining the probability of real eigenvalues are the degree of smoothness and the existence of finite moments. It is noteworthy that these (apart from zero mean) are also the crucial features for universality in the case of random polynomials Edelman and Kostlan (1995). Given that the distributions are smooth and have finite moments, we tentatively propose that the normal distribution is the one with the maximum probability of finding real eigenvalues. If true, this provides a further unique characterization of the normal.

Besides the bounds, we also observe a hierarchy for the fraction of real eigenvalues for matrices constructed from these distributions. Specifically, for some commonly occurring distributions, we find that the eigenvalues are more real when the matrix elements are chosen from Cauchy distribution than Laplace, Gaussian or uniform distributions, which are arranged in the increasing order of the decay. However, we must mention that such hierarchy patterns are not simply determined by the tail behavior of the distribution, and the trends appear to be more complex as explained in the following section.

In Sec. III, we numerically study the properties of the eigenvalues of the matrix obtained after taking product of several matrices, and find that the increase in the probability to one holds for several other distributions. In particular, we explore the zero-centered uniform distribution, the symmetric exponential (Laplace) distribution and Cauchy distribution, as simple representative ones. The hierarchy mentioned above is seen to hold even after the multiplication of independent random matrices. The expected number of real eigenvalues is also studied which tends to the matrix dimension exponentially fast (with the number of matrices in the product), although the rate of approach is distribution-dependent.

It is unclear as to why the eigenvalues tend to become real when multiplying random structureless matrices. In an attempt at seeing how much this has to do with the act of multiplication, we study numerically the probability of real eigenvalues in the case of Hadamard (or Schur) products where the elements are simply the products of the corresponding elements of the multiplying matrices. Of course, in this case, the matrix elements remain uncorrelated. We find numerically that the probability of real eigenvalues increases with the number of matrices in the product, but tends to saturate at a value smaller than . Numerical evidence that this approach is a power law is also provided. This then highlights that the correlations built up in the process of (usual) matrix multiplication are responsible for the phenomenon that the fraction of real eigenvalues tends to one.

## Ii Real eigenvalues of a single 2×2 matrix

Let denote the probability that a product of random matrices has real eigenvalues. In this section, we study the simplest case, viz., the probability that all the eigenvalues of a matrix are real. We assume that the matrix elements are i.i.d. random variables chosen from the distribution with support on the interval , where is finite for bounded distributions and infinity for unbounded ones, and that the distribution is symmetric about the origin as a result of which the mean of the probability distribution is guaranteed to be zero. We first give analytical results for the probability , which we study for many probability distributions, and then present some numerical results for the more general quantity in the following section.

Consider a matrix with i.i.d. elements defined as

As the discriminant must be nonnegative for real eigenvalues, the probability that all eigenvalues are real is

 P(1)2,2=∫u−u∫u−u∫u−u∫u−uΘ[(v−w)2+4xy] p(v)p(w)p(x)p(y) dvdwdxdy. (2)

Let the probability distribution of be with . Observe then that the probability of real eigenvalues can be written as

 P(1)2,2=12+4∫0−udx∫u0dy∫2u0dz q(z)p(x)p(y)Θ[z2+4xy]. (3)

Here the even symmetry of the distributions is used, and also that when the signs of both and are the same, the constraint imposed by the Heaviside function is trivially satisfied. In the above equation, the distribution is given by the convolution

 q(z)=∫u−u p(v) p(w) δ(z−(v−w)) dvdw=∫u−u p(v)p(v−z) dv. (4)

We note that this is indeed the usual convolution with due to the symmetry of . If the above expression for is valid, else it is zero.

Equation (3) shows that the probability of both eigenvalues being real is at least one half for any probability distribution. The integral in Eq. (3) can be simplified if we use the convolution as

 P(1)2,2=12+4∫u0dx∫u0dy∫2u2√xydz q(z)p(x)p(y)=1−4∫u0dx∫u0dy∫2√xy0dz q(z)p(x)p(y). (5)

The last equality follows from the fact that the convolution is itself a symmetric normalized distribution as . Thus the evaluation of the probability of real eigenvalues reduces to evaluating the triple integral above.

It is also useful to write the distribution as

 q(z)=12π∫∞−∞dk eikz |~p(k)|2,~p(k)=∫u−up(x) e−ikx dx. (6)

Here is the characteristic function (or Fourier transform) of the probability distribution . Next consider the integral defined as

 I(α)=∫0−udx∫u0dy∫2u0dz q(z)p(x)p(y)Θ[z2+αxy]. (7)

The value of interest is . Towards this end, differentiating with respect to (w.r.t.) , and performing the integral over the resulting delta function leads to

 ∂I∂α=∫∞−∞dω∫u0dx∫u0dyxy2√αxy|~p(ω)|22πeiω√αxyp(x)p(y), (8)

where is given by Eq. (6). Noting that the above equation is valid only for , we integrate the last expression w.r.t. from to , and use to get and finally the probability that both the eigenvalues are real as

 P(1)2,2=1−4π∫∞0dω∫u0dx∫u0dy|~p(ω)|2sin(2ω√xy)ωp(x)p(y). (9)

Below we will apply the result in either Eq. (5) or Eq. (9) to various choices of distribution . Before turning to explicit calculations, we note that the unbounded distributions have the following scale invariance property. Consider a division of the random matrix elements by a nonzero constant . If , then where and . From Eq. (2), it is seen that the probability is not affected by this scaling for unbounded functions, but for bounded ones, the limits should be redefined. While the above is an elementary analysis, we are not aware of it being discussed before.

Under conditions of smoothness (in the sense that all the derivatives exist at least in intervals) and symmetry of , it is interesting to enquire about the maximum value or upper bounds of , for instance. For the lower bound, we have already stated that the probability of all real eigenvalues is at least one half. A tighter lower bound is however , a proof of which was suggested to us by an anonymous referee which we discuss now. The probability of eigenvalues being real is

 P[(v−w)2+4xy>0]=12P[(v−w)2+4xy>0|xy>0]+12P[(v−w)2+4xy>0|xy<0], (10)

where the factors of arise as probability that or otherwise. However the first conditional probability is , and therefore

 P[(v−w)2+4xy>0] = 12+12P[(v−w)2−4xy>0|xy>0] = 12+14P[(v+w)2−4vw−4xy>0|xy>0,vw>0] + 14P[(v+w)2−4vw−4xy>0|xy>0,vw<0] = 12+14P[(v+w)2−4vw−4xy>0|xy>0,vw>0] + 18P[(v+w)2−4vw−4xy>0|xy>0,vw<0,|vw|

Here, in the second equality, further conditions on the sign of are used with equal probabilities for either occurrence (from the symmetry and independence of the distributions). In the third equality, a further certainty is carved out by using the possibility that when (occurence probability being ), the discriminant is certainly positive. However our attempts at a similar approach for the upper bound were not successful.

Although one can obtain get some insight into the bounds, for what kind of distributions these bounds actually occur is discussed below.

### ii.1 Bounded distributions

#### ii.1.1 Symmetric Beta distribution

Consider the symmetric Beta distribution with zero mean defined as

 pμ,ν(x)=Γ(μ+ν+2)2Γ(1+μ)Γ(1+ν) (1−|x|)μ |x|ν , μ,ν>−1 (15)

with support on the interval . This is a two-parameter symmetric family which can be nonsmooth at . For , the above distribution reduces to a uniform distribution, while and correspond to tent-shaped or V-shaped distribution respectively.

Case:

For the uniform distribution, the convolution on and zero elsewhere. Using Eq. (5), we find that

 P(1)2,2(ν=0)=12+∫10dx∫10dy∫2√4xydz q(z)=4972=0.680556. (16)

Thus the uniform distribution results in a smaller fraction of real eigenvalues when compared to the normal distribution Edelman et al. (1994), but not by very much.

When and is nonzero and positive, the distribution is zero at the origin. As tends to infinity, the weight of the distribution gets increasingly concentrated at and as Table 1 shows, the probability of real eigenvalues decreases. The distribution is smooth when is an even integer. In this case, an analytical expression is possible for arbitrary , and given by Eq. (38) of Appendix A. Numerical analysis of Eq. (38) for large strongly suggests that

 P(1)2,2(0,ν)=58+O(ν−1),ν≫1. (17)

An insight into the above result can be gained by the following heuristic arguments: If we consider the limiting distribution as the Bernoulli ensemble with only two choices of entries with equal probabilities, we get only matrices in the ensemble of which have real eigenvalues. Thus it would seem that the ratio should have converged to . However, there are cases where both the eigenvalues are exactly zero. If we believe that in the case of continuous distributions these are modified into cases with real and complex eigenvalues and equally, we get cases of complex values and the fraction is consistent with for real eigenvalues. Another related heuristic argument, closer to the continuous distributions under consideration, starts with the distribution

 p(x)=12(δ(x+1)+δ(x−1)), (18)

and applies Eq. (9) with , leading to

 P(1)2,2=1−2π∫∞0sinωcos3ωωdω=58.

A proof of the assertion in Eq. (17) by analysing the double sums in Eq. (38) seems rather difficult to obtain, but it is possible to tackle the Gamma distribution given by Eq. (24), for which also the probability is zero at the origin and increases algebraically away from it, and show that indeed the limit probability is , see Sec. II.2.2. Thus the lower bound derived above is a tight one.

We now turn to the case when where the probability distribution piles up at the origin. As shown in Table 1, the probability of real eigenvalues increases as decreases towards . It is obvious that if the matrix elements are chosen from a Dirac-delta distribution centred about zero, both the eigenvalues are definitely real. But whether the probability limits to something less than one as approaches is of natural interest. As shown in Appendix A, we find that

 P(1)2,2(0,ν)=78−O(1+ν),  ν→−1. (19)

Thus it is interesting that the distribution restricted to an interval spans a range of behaviors for the probability of real eigenvalues. As the matrix elements probability gets increasingly piled up at the ends (), the probability of real eigenvalues decreases and tends to , while in the opposite case when the elements are piled up around origin, the probability approaches the maximum value of .

Case:

For , the probability distribution and we find that with increasing , the probability that the eigenvalues are real increases, refer Table 2. Partial analytical results are possible. For example, for the case when (“tent” distribution), the convolution is

 q(z)=⎧⎨⎩16(2−|z|)31≤|z|≤216(4−6|z|2+3|z|3)0≤|z|≤1.. (20)

The resulting integral in Eq. (5) can now be done using hyperbolic coordinates. This results in

 P(1)2,2(μ=1,ν=0)=1614322400−23ln21008≈0.704854. (21)

It is reasonable to expect that as increases (and ), the probability increases to that when the elements are distributed according to the Laplace distribution which is obtained as . We will deal with the Laplace distribution below, but state here that the probability of real eigenvalues in this case is . However, for negative where the matrix elements have a tendency to be close to , from the discussion in the last subsection, we expect the probability to approach as . Thus for the distribution , the probability of real eigenvalues lies in the range .

#### ii.1.2 A smooth family

Consider another class of smooth zero mean distributions bounded on the interval defined as

 pη(x)=Γ(η+32)√πΓ(η+1)(1−x2)η , η>−1. (22)

Of course, for , the distribution is continuous but not smooth at , but this does not seem crucial. When , is uniform on the said interval, and as the distribution approaches the normal distribution with the variance scaling as . As discussed above, the probability of real eigenvalues is independent of the variance and therefore, we expect that the large value for this probability will coincide with the known result for the normal distribution, namely, Edelman et al. (1994) (also, see Sec. II.2.1). The probability of real eigenvalues is calculated for some values of integer in Appendix B, and we find that the probability indeed approaches with increasing .

When is negative, the distribution diverges at , and from the discussion in the preceding subsection, we expect the probability of having real eigenvalues to approach as . The case of is that of the arcsine distribution. The convolution with itself can be found and is a complete Elliptic integral. However, even if further analytic results seem to be hard, this enables a more accurate numerical estimate of the probability of real eigenvalues which is . Further decreasing makes the numerical evaluation unstable, but results indicate a monotonic decrease in the probability.

### ii.2 Distributions with infinite support

#### ii.2.1 Gaussian

The case of normal or Gaussian distribution is the most studied and there are general results Edelman et al. (1994); Lakshminarayan (2013); Forrester (2014). We consider a zero centered, unit variance Gaussian distribution for the matrix elements given by

 p(x)=e−x2/2√2π;x∈(−∞,∞). (23)

The probability of both eigenvalues being real has already been shown, using Eq. (9), to be exactly equal to in Lakshminarayan (2013), and alternative derivations naturally exist, see the works of Edelman (Edelman et al., 1994) and Forrester Forrester (2014). However to place it in the context of this paper, we rederive the probability of real eigenvalues in this case in Appendix C.

It is well known that the normal or Gaussian distribution is singled out in numerous ways, for example, as one that maximizes entropy for a given mean and variance, or as the limit of sums of random variables. In the context of the present work, it seems plausible that the normal distribution is once again to be singled out as the distribution of matrix elements that maximizes the probability of finding real eigenvalues among the class of smooth, symmetric, and finite moments distributions. To test this further, we have looked at distributions such as , ( etc. and verified in all these cases that the probability of real eigenvalues is indeed less than . We advance this proposition tentatively based on evidence gathered so far, and from results to be presented below.

#### ii.2.2 Gamma distribution

We next consider the symmetrized Gamma distribution defined as

 pγ(x)=12Γ(γ)|x|γ−1e−|x| , γ>0. (24)

where the scale for the exponential decay has been set to one. The special case of Laplace distribution for which and the general case using the convolution route are discussed in Appendix  D. The convolution itself is given by

 q(z)=14Γ(2γ)e−|z||z|2γ−1+12γ+12√πΓ(γ)|z|γ−12Kγ−12(|z|), (25)

where is the modified Bessel function of the second kind DLMF . The results obtained by analysing the resulting integrals are shown in Table 3, and we find that with increasing , as the weight of the distribution at the origin decreases, the probability also decreases.

The behavior of the probability of real eigenvalues as approaches zero can be understood as follows. As approaches zero from the positive side, the distribution and both the terms in the convolution diverge at the origin as . The latter can be seen as , and , as follows from the duplication formula for the Gamma function Abramowitz and Stegun (1964). Since we expect that the divergence at the origin is all that matters, the probability for real eigenvalues will converge to the case of the bounded distribution already studied above, and therefore the probability for real eigenvalues increases to as .

For increasing , the distribution itself is peaked away from the origin and has two symmetric maxima at . The probability now decreases from the value and monotonically seems to approach , the lower limit for the bounded distributions as increased. Indeed the piling up of the probability of the matrix elements at two symmetric “walls” makes this plausible. To analyze this further, we restrict attention to integer values of and give an exact evaluation in terms of finite sums as

 P(1)2,2(γ)=12+A1(γ)+A2(γ), (26)

where and are given by Eq. (62) and originate from the first and second terms of the convolution in Eq. (25) respectively. To find the limiting value for the probability of real eigenvalues for large , we note that like the distribution , the first term in the convolution is also peaked around , but the second term peaks at so that the two terms in the convolution have practically disjoint supports. As a consequence, for large , the dominant contribution to the probability in Eq. (26) comes from which represents the overlap between the distribution and the convolution, and can be neglected. As discussed in Appendix D, on analysing the sum , we get

 P(1)2,2(γ)=58+116√2πγ+O(γ−1),γ≫1. (27)

In summary, for the symmetrized Gamma distributions, the probability of real eigenvalues also seems to be in the range . The distribution having a nonanalyticity at the origin leads to larger probability for real eigenvalues compared to the normal distribution and it seems comparable to the bounded power law distributions studied above except in the rates of convergence.

#### ii.2.3 Power law distributions

A qualitatively different process is interesting to consider, and as in the case of random polynomials Logan and Shepp (1966), it will be interesting to study what happens when the underlying probability distributions have diverging moments. The Cauchy distribution which is the simplest and best studied of these and occurs in many contexts, is given by

 p(x)=1π(1+x2) (28)

for . Again the possible parameter in the distribution is rendered inoperative via scaling. Using and performing a change of variables in (9), the integral w.r.t. can be calculated. Then computing the integral w.r.t. using series expansion yields

 P(1)2,2=1−4π2∫∞011+y2tan−1⎛⎜ ⎜⎝11+√2y⎞⎟ ⎟⎠dy=34. (29)

Although the integral is written as if it is carried out, it is in fact a numerical evaluation which is almost certainly correct. Evaluations using the convolution path lead to interesting alternative integral forms of , but none of them (including the above) seem to be either in standard tables or calculable using symbolic mathematical packages.

In general, one may consider the distribution which possesses finite (and nonzero) moments up to order only. Numerical evaluation of Eq. (5) gives and for and respectively, both of which are smaller than the result obtained above for the Cauchy distribution which is the slowest decaying power law distribution with finite mean. We also note that compared to the Gaussian case, the eigenvalues are less likely be real for the power law-distributed matrix elements with , and therefore the probability is not merely determined by the decay behavior of the distribution of the matrix elements. Other distributions such as which is slower than Cauchy distribution but not smooth, or power laws corrections that vanish or diverge at the origin to the fat-tailed distributions have not been investigated.

To summarize, the case of a single matrix with elements drawn from various i.i.d. distributions have revealed an interesting phenomenology. It may seem like the weight of the distributions near the origin matters and if the values are clustered around zero, the probability of real eigenvalues increases. This statement must however be qualified: the normal distribution with however a small variance always have only probability of having real eigenvalues. Thus the differentiability of the underlying distributions play an important role. In the extreme case of delta distributions, we may have eigenvalues real. But given that the distributions be smooth to all orders, the normal distribution seems to be singled out as the one with the largest probability of real eigenvalues. Also, the probability of real eigenvalues seems to be in the range for the classes of distributions considered here. For the case of the Cauchy distribution which is smooth but has diverging variance, the probability is large at for real eigenvalues. Thus this rather simple problem of the probability of real eigenvalues of random real matrices seems to possess a multitude of interesting features that warrants further study and clarification.

## Iii Real eigenvalues of product of matrices

We now turn to the properties of an matrix obtained after taking a product of square matrices, and study how the probability that some or all of the eigenvalues are real behaves for . It is of interest to see if the results for single matrix described in the last section carry over to the product of matrices. Thus there is a two-fold generalization: (1) products of matrices are considered, and (2) their dimensionality can be more than two.

### iii.1 Asymptotic value and maintenance of hierarchy

We numerically studied the probability that eigenvalues are real for a product of random matrices. As the products of random matrices can have, in general, a positive Lyapunov exponent Furstenberg and Kesten (1960), numerical procedures for all cases of products “renormalizes” the matrices for each product by dividing the Frobenius norm. This however does not alter the quantities of interest here. For most of the discussion, the matrix elements are assumed to be i.i.d. random variables distributed according to one of the following symmetric probability distributions: uniform, Gaussian, Laplace and Cauchy. Note that all of these distributions are finite at the origin, but have different tail behavior. The data were averaged over independent realisations of the random matrices. For the matrix, there can be either zero or two real eigenvalues, but for higher dimensional matrix with that we consider here, the probability of real eigenvalues is nonzero for and only. Numerical results for the probability of all real eigenvalues for and are presented in Fig. (1) for various probability distributions as a function of the number of matrices in the product. We find that the probability of all eigenvalues being real increases to unity with monotonically for most distributions (see, however, the case of Gamma distribution defined by (24) with ). This effect has been previously observed by one of the authors Lakshminarayan (2013) for Gaussian-distributed matrix elements. Here we find that this result is quite general in that the eigenvalues tend to become real with increasing when matrix elements are distributed according to non-normal distributions as well.

The second important point illustrated by Fig. (1) is that the same hierarchy as observed for the case continues to hold for . Explicitly, the probability of all eigenvalues being real for the uniform distribution is the smallest followed by the Gaussian distribution, the Laplace distribution and finally the Cauchy distribution. Thus the slowly decaying distributions appear to have larger probability of having all real eigenvalues. This feature is seen to be valid for higher dimensions as well (results not presented). Thus the hierarchy apparent even with a single matrix continues to hold for larger dimensional matrices as well as the product of the random matrices. Furthermore, Fig. (2) shows that the ordering of the probability of real eigenvalues according to their tail behavior is not special to the probability of all real eigenvalues, but holds when as well. Note, however, while the probability of all real eigenvalues increases monotonically towards unity, the probability of real eigenvalues decays to zero, as the number of matrices in the product increases (see, the inset of Fig. (2)).

The increase in the probability of all real eigenvalues is reflected in the average number of real eigenvalues as well which is given by . Since, as discussed above, the probability as increases, the average approaches the dimension of the matrix. Numerical results in Fig. (3) for matrices strongly suggest an exponential approach of the expected number of real eigenvalues to the dimension of the matrix. The data also points to the continuance of the hierarchy even for the average number of real eigenvalues indicating a more microscopic adoption of the hierarchy to all .

### iii.2 Effect of correlations between matrix elements

We now consider the Hadamard (or Schur) product wherein the elements of the product are simply the products of the corresponding elements. Thus if is used to denote the Hadamard product of two matrices and , then . The Hadamard product of two positive semidefinite matrices is also positive semidefinite, unlike the usual product, and is one of the reasons it is important Horn and Johnson (1987). However the reason why we choose to study Hadamard products is that although the matrix elements of the product matrix are products of random numbers, there is no correlation amongst the matrix elements. Thus one can hope to disentangle two possible mechanisms that may be responsible for the phenomenon that all the eigenvalues of a random product matrix tend to be real: the act of simple multiplication which is presumably causing the matrix elements to have large weight near the origin and the various addition of such products that are leading to correlations between the matrix elements.

To this end, we consider the eigenvalue properties of a Hadamard product of random matrices with elements distributed according to uniform, Laplace, Gaussian and Cauchy distributions, each with zero mean. Figure (4) shows the case of Hadamard products of matrices, and we observe that the hierarchy effect is maintained at any (although there are more fluctuations in this case as compared to that of the ordinary matrix product). But, importantly, we find that the probability does not approach unity with increasing .

We can get some insight into this result by calculating the distribution of the matrix elements of the Hadamard product matrix for some distributions, and appealing to the results for the matrices obtained earlier in Sec. II. Let us first consider the case of Hadamard product of matrices with elements drawn from the bounded distribution defined in Eq. (15). Let represent a random variable formed by the product of i.i.d. random variables. As discussed in Appendix E, the distribution of the product for this case is given by

 pK(zK=x1x2...xK)=(ν+1)K|zK|ν2(K−1)! [ln(1|zK|)]K−1Θ(1−|zK|). (30)

For , the above distribution diverges at the origin while for positive , it vanishes at and and is a nonmonotonic function of . Thus, except for the uniform case, the behavior of the distribution is similar to that of near the origin. Using the above result, it is possible to numerically evaluate the integral in (5) to obtain for various values of and . For example, for , we find the probability of real eigenvalues to be , , , , , for respectively. For larger , numerical integration does not converge; however, as shown in Fig. 4 for some representative values of , the data obtained from direct sampling indicates the approach to a probability less than . For the case of uniform distribution, the probability seems to saturate around . But for negative and positive , the probability of real eigenvalues approaches a value higher and lower than respectively. This pattern is consistent with the results in Sec. II where the probability of real eigenvalues decreased with increasing .

We next consider unbounded distributions, and start with the case of Hadamard product of matrices with Gaussian distributed elements with zero mean and unit variance for which Springer and Thompson (1970); Forrester (2014)

 p2(z2)= 12π∫∞−∞∫∞−∞e−x21/2e−x22/2 δ(z2−x1x2) dx1dx2=12π∫∞−∞e−|z2|coshx1 dx1=K0(|z2|)π, z2∈(−∞,∞). (31)

which diverges logarithmically at the origin since as approaches zero. For this distribution, the integral in Eq. (5) can be numerically evaluated by using the fact that its convolution is given by a Laplace distribution (shown in Sec. III B 2) yielding , which is consistent with the value obtained in Fig. (4) by direct sampling. For , the probability distributions are given by the Meijer-G functions Springer and Thompson (1970); Forrester (2014) making it difficult to handle this using Eq. (5). We therefore give the results obtained via numerical sampling in Fig. 4 which again indicates an approach towards . In the case of Hadamard product of two matrices with Laplace-distributed elements with zero mean given by (24) with , we have Springer and Thompson (1970)

 p2(z2)= 14∫∞−∞∫∞−∞e−|x1|e−|x2| δ(z2−x1x2) dx1dx2=12∫∞−∞e−2√|z2|coshx1 dx1=K0(2√|z2|), z2∈(−∞,∞), (32)

where the last step in Eq. (32) has been done by a series of two transformations, followed by . Numerical evaluation of Eq. (5) using this yields , and the results obtained from direct sampling show a convergence towards a value close to . For the Cauchy-distributed matrix elements, the distribution of the product of two random variables is known to be Springer and Thompson (1966)

 p2(z2)=lnz22π2(z22−1), (33)

which diverges at the origin logarithmically.

Thus, as Fig. (4) shows, the probability of all real eigenvalues monotonically increases with and saturates to a value in the range . Note that these numbers are less than which, from our analysis in the previous section, is expected of distributions that diverge as a power law, with , at the origin. From the above calculations for , it is clear that the uniform distribution and the three unbounded distributions discussed above have Hadamard product elements distributed according to a probability distribution that diverges as at the origin. In fact, the elements of the Hadamard product of uniform distribution as well as Cauchy distributed random matrices have probability distributions that diverge as at arbitrary Springer and Thompson (1966). It is reasonable to expect that the Gaussian and Laplace cases also diverge in a similar manner. For these cases, as the divergence at the origin is somewhat “slower” than a power law, it seems reasonable that the probability of eigenvalues being real is a shade smaller than that of single matrices whose elements are distributed according to the power law discussed above.

Assuming that the probability of Hadamard products saturate to the same value for these four distributions in the limit, and taking this to be (in the case of ), our numerical results show a power law approach to the constant:

 P(K)2,2=P(∞)2,2−CKθ, (34)

where is a positive constant and the exponent for uniform, Gaussian, Laplace and Cauchy distribution respectively (see the inset of Fig. (4)). Thus, here the probability of real eigenvalues approaches the asymptotic value as a power law in comparison to the ordinary matrix product case, where this approach is seen to be exponentially fast Lakshminarayan (2013). The exponential behavior is also seen in the expected values of real eigenvalues as in Fig. 3. We also looked at the fraction of expected number of real eigenvalues for higher dimensional matrices () when each matrix in the Hadamard product has Gaussian-distributed elements with zero mean, and find that the fraction of real eigenvalues increases and saturates to a value less than one, unlike that observed for the case of usual matrix products. We also find that the asymptotic value of this fraction decreases with an increase in the dimensionality of the matrices.

The convergence of the probability that all eigenvalues are real as for the Hadamard product can also be understood using the Central Limit Theorem applied to the logarithm of the absolute value of the product random variable. The elements of the Hadamard product matrix would then admit a symmetrized log-normal probability distribution of the form Sornette (2000)

 pK(zK=x1x2...xK)=12|zK|√2πKσ2 exp[−(log|zK|−Kμ)22Kσ2], zK∈(−∞,∞), (35)

for large , with the parameters and being the mean and standard deviation of . Numerically, we find that the probability that all eigenvalues are real for a matrix with elements distributed according to Eq. (35) indeed converges to a value around at large , as for the Hadamard product.

#### iii.2.2 Usual matrix product

We now return to the properties of the matrix obtained after multiplying matrices, and ask why the asymptotic probability of real eigenvalues saturates to unity unlike that for the Hadamard product. We can think of two possible reasons: one, the (usual) product matrix elements have probability distributions that are significantly different from the Hadamard case. The other possibility is that the correlations between the matrix elements in the case of usual products lead to an increase in probability of all eigenvalues real to unity as opposed to the Hadamard case where this is not true. It is not obvious whether the difference in the probability distribution or the correlations between the matrix elements is responsible for a higher .

To gain some insight into this intriguing question, we looked at the matrices with i.i.d. random variables as elements, each distributed according to the probability distribution of the product matrix elements at various . Such matrices can be easily generated by taking independent samples of usual product matrices at each and forming new matrices by picking the corresponding matrix elements from the independent product matrix samples. These new matrices would then have elements with probability distributions corresponding to that of elements of a usual product matrix at product length , but the matrix elements would now be independent, having been selected from independent product matrix samples. Fig. (5) shows a comparison between for the Hadamard product and usual product and for matrices with probability distribution of elements corresponding to that of usual product at product length , but without any correlations between the elements. It is clear that the difference in probability distribution of matrix elements between the Hadamard and the usual product has an effect of lowering for the usual product case to a value below that of the Hadamard product case. However, the presence of correlations between the matrix elements of the usual product matrix is seen to have an effect of a substantial increase in , finally leading to an asymptotic value of one.

The effect of correlations can be seen analytically for the case of Gaussian-distributed matrices when . We have seen that the distribution of the product of two Gaussian-distributed random variables with zero mean and unit variance is given by Eq. (31) and the probability of real eigenvalues is . The distribution of the matrix elements obtained on taking the usual matrix product corresponds to that of the sum of two random variables, each distributed according to the distribution of the product of two random variables. The characteristic function of Eq. (31) is given by (due to Eq. (11.4.14) of Abramowitz and Stegun (1964)). It is easily seen that this is the square root of the characteristic function of Laplace distribution. Hence, the product matrix elements have a Laplace distribution, for which we know from Eq. (56) that would have been , had the elements been independent. But, the correlations between the matrix elements raise this probability to Lakshminarayan (2013); Forrester (2014) for the usual product, hence presenting a strong evidence that the correlations between the matrix elements are in fact leading to observed for usual products.

We also measured the correlations between the matrix elements to see how correlations increase with increasing for matrices with Gaussian-distributed matrix elements. More precisely, we consider

 Ci=(⟨x1xi⟩−⟨x1⟩⟨xi⟩K)1/K , i=1,...4 (36)

where are the matrix elements obtained after taking the product of matrices. For , obviously . But with increasing , all the four ’s seem to be of same order i.e. the variance of the new distribution and the correlations are similar. In particular, we find that for , the correlations for respectively, but these numbers increase to and for and respectively. These data suggest that similar to the probability , the correlations also tend to saturate with increasing .

## Iv Discussion

How does the probability of real eigenvalues for a random matrix depend on the probability distribution of the matrix? It is quite surprising that such an apparently elementary question presents many novel challenges, and throws up some surprises. Here we find that the probability of real eigenvalues depends on several detailed features of the probability distribution, some of which we have identified here as the finiteness and smoothness at the origin, existence of maxima away from origin and the finitness of the moments. We have shown that eigenvalues are most likely to be real for distributions limiting to which is the most divergent distribution at the origin (with zero mean) and occur with a probability . That the large weight at the origin correspond to high probability is perhaps not surprising since when the matrix elements are distributed according to , this probability is unity. However, interestingly, here we find the maximum probability to be less than one. The origin of this naturally rests in the smoothness of the probability distributions considered. We also mention that the values are not universal; for example, the probability turns out to be different for the Beta distribution and for Gamma distributed matrix elements with , although both have the same behavior near .

Our results suggest that the probability that all eigenvalues are real is larger for distributions with large weight at the origin and that decay slowly. More concretely, there exists a hierarchy between the probability values of the different distributions, which is Cauchy (0.750) Laplace (0.733) Gaussian (0.707) uniform (0.6805). Moreover, for class of distributions that decay in the same manner, the distributions with higher weight close to zero have a higher probability of real eigenvalues. For example, although the probability for distribution Eq. (24) for is larger than the uniform distribution as it decays slower than the bounded distribution, this probability for is found to be smaller than that for the uniform case. Thus the ‘reality’ is determined by both the tail and near-zero behavior of the probability distribution from which the matrix elements have been chosen. For the broad class of smooth distributions we have considered, we did not find the probability of real eigenvalues to lie outside the range , which appear to be natural boundaries for sufficiently smooth distributions.

Numerical results discussed above also shows that this remains valid even after taking a product of random matrices. It is quite surprising that the curves for probabilities of the different distributions do not cross each other. An intuitive explanation of this is yet to be figured out. We also measured the distribution of the matrix elements of the product matrix, and find that it decays slower (but faster than any power law) as the number in the product increases. This observation is consistent with the result for single matrices (with finite moments), namely, that faster the distribution decays, smaller is the probability that all eigenvalues are real. However this is only part of the explanation as the probability always stays below unity when correlations between matrix elements are ignored, and in order to approach unity, correlations are found to be crucial. It maybe stated that despite extensive and exact results obtained so far there is not much understanding about why the eigenvalues tend to become real under the usual matrix product, and more work in this direction is clearly desired.

###### Acknowledgements.
AL thanks JNCASR, Bangalore, for supporting a short visit facilitating discussions. It is a pleasure to thank an anonymous referee for suggesting the proof that is the lower-bound.

## Appendix A Detailed derivations for symmetric beta distribution

When and is a positive integer, the convolution can be written in the form of a finite series:

 q(z)=Θ(2−|z|)(2k+12)22k∑r=0(2kr)(−|z|)2k−r1−(|z|−1)2k+r+12k+r+1. (37)

Integrations can be done using the hyperbolic variables as the combination appears exclusively. To outline the method used, hyperbolic coordinates () where and are useful. Thus the variable does not appear in the transformed expression. The range of is , while for a given , the range of is . The integration over is easily carried out first, and noting that the Jacobian of the transformation is , leads finally to:

 P(1)2,2(μ=0,ν=2k)=1−(2k+1)4×2k∑r=0[(2kr)22k−r+12k+r+1((−1)r(2k−r+1)(6k−r+3)2+2k+r+1∑l=0(2k+r+1l)(−2)l(2k−r+l+1)(6k−r+l+3)2)]. (38)

As special cases:

 P(1)2,2(μ=0,ν=2)=890514112=0.631023,P(1)2,2(μ=0,ν=4)=4533248972144072≈0.628361. (39)

The case of odd integer is also accessible, however we only state the result when . The convolution is now piecewise continuous on and , and the probability of real eigenvalues is

 P(1)2,2(μ=0,ν=1)=36535760+ln2240≈0.63709. (40)

To understand the behavior of the probability of real eigenvalues for negative , it is convenient to write the convolution for the distribution defined by Eq. (15) as

 q(z)=∫+∞−∞p(x)p(z−x)dx=(ν+12)2∫+1−1|x|ν|z−x|νΘ(1−|z−x|)dx. (41)

Note that and implies that . From the symmetry it suffices to consider . Using the hyperbolic coordinates, and , we get

 P(1)2,2=1−(ν+1)2∫102vdvv2ν∫−lnvlnvdw∫2v0dzq(z)=1+(1+ν)4∫10v1+2νlnv∫2v0dz∫1−1|x|ν|x−z|νΘ(1−|x−z|)dxdv. (42)

Except when is an even integer, some care must be taken as is piecewise continuous in the two intervals and . For , we have

 P(1)2,2(μ=0,ν=−1/2)=148(41−π−2ln2)≈0.759836. (43)

To find the probability as , consider the two inner integrals over and in Eq. (42). On considering various cases arising due to the Heaviside theta function and performing the integral over first, we get

 P(1)2,2(μ=0,ν)=34−Γ(12)Γ(2+ν)Γ(32+ν) 1+4(1+ν)ln242(2+ν)−J1+J2, (44)

where

 J1=(1+ν)2(∫1/20dv v1+2νlnv (1−2v)1+ν−∫11/2dv v1+2νlnv (2v−1)1+ν), (45)

and

 J2=(1+ν)3(2∫1/20dv v2+2ν