Spectral properties of the generalized diluted Wishart ensemble

The celebrated Marčenko–Pastur law, that considers the asymptotic spectral density of random covariance matrices, has found a great number of applications in physics, biology, economics, engineering, among others. Here, using techniques from statistical mechanics of spin glasses, we derive simple formulas concerning the spectral density of generalized diluted Wishart matrices. These are defined as F≡12dXYT+YXT , where X and Y are diluted N × P rectangular matrices, whose entries correspond to the links of doubly-weighted random bipartite Poissonian graphs following the distribution P(xiμ,yiμ)=dNϱ(xiμ,yiμ)+1−dNδxiμ,0δyiμ,0 , with the probability density ϱ(x, y) controlling the correlation between the matrices entries of X and Y . Our results cover several interesting cases by varying the parameters of the matrix ensemble, namely, the dilution of the graph d, the rectangularity of the matrices α = N/P, and the degree of correlation of the matrix entries via the density ϱ(x, y). Finally, we compare our findings to numerical diagonalisation showing excellent agreement.


Introduction
The Wishart distribution was introduced in 1928 by John Wishart [1] and it can be thought of as a generalization to a higher dimensional gamma distribution. Along the years it has become one of the central matrix ensembles in random matrix theory, the so-called Wishart ensemble. It naturally arises as the probability function of sample covariance matrices for a sample of multivariate normal distributions and, as such, has found a great number of applications in the context of multivariate statistical analysis.
From the several statistical properties derived from the Wishart ensemble, the limiting distribution of its spectral density, the celebrated Marčenko-Pastur [2] law, can be ubiquitously found in multiple applications ranging from physics, to biology, to economics, to engineering, among others. More precisely, let us consider a set {x μ i } of independent and identically distributed random variables with zero mean and unit variance, with i = 1, . . . , N and μ = 1, . . . , P. This may be understood as N times series each of size P. The entries of the N × N sample covariance matrix E are defined as: where we have introduced an N × P matrix X whose elements are precisely (X) iμ = x μ i . Notice that in the limit P → ∞, while keeping N fixed, the covariance matrix becomes the identity matrix and its spectral density is a Dirac delta distribution centered at one. However, in the limiting case of P → ∞ and N → ∞, while keeping the ratio α = N/P fixed, the spectral density ρ(λ) is described the Marčenko-Pastur law [2]: This ensemble was recently generalized in [3] by considering the symmetric cross-correlation ensemble, composed of matrices of the following form: where X and Y are two N × P rectangular matrices, each with independent and identically distributed zero mean and unit variance entries, but such that E[x μ i y μ j ] = cδ ij , so that c ∈ [−1, 1] controls the correlation between the elements of both matrices. Obviously, for c = 1 one recovers the classical Wishart ensemble. The ensemble, as pointed out in [3], can be seen, for instance, as a null hypothesis for detecting correlations between different multidimensional time series, or as the statistical properties of the anticommutator of random matrices for N = P. Its spectral density was derived in [3] using free probability theory.
On the other hand, random graphs constitute the main tool to model the complex behaviour of large empirical networks observed in social, technological, and biological systems [4,5]. In random graph models a network is typically represented by nodes that interact through edges. Random graph theory leads to important insights into the structure of networks as well as on the dynamical processes occurring on them, such as the spreading of diseases [6,7], the stability of ecosystems to perturbations [8], the dynamics and equilibrium properties of sparsely connected neurons [9][10][11], models of formation of opinion [12], and so on. Being the spectral properties of random graphs a central property to understand these systems, in a series of seminal works the spectral density of standard, fully connected, random matrix ensembles was generalised to diluted random graphs, first for symmetric matrices in [13], subsequently for asymmetric ones [14], which was later followed by correlated networks [15]. Moreover, recently a set of powerful mathematical methods have also been introduced to obtain the large deviation properties of diluted random matrix ensembles [16][17][18].
The main goal of the present work is to introduce the diluted symmetric cross-correlation ensemble and to study its spectral density using mathematical tools originated in statistical mechanics of spin glasses [19]. This work is organized as follows. In section 2 we start by defining this ensemble, discuss how by changing the parameters of the ensemble we cover previous cases, and go through the main derivations to obtain the spectral density by using the cavity method. In section 3 we summarize our theoretical results, exploring some particular cases, and compare them to numerical diagonalization. We conclude with a summary and discuss future lines of research.

Definitions and main theoretical derivations
We start by introducing the diluted symmetric cross-correlation ensemble, which we will also call as the generalized diluted Wishart ensemble. To do this, let us take two N × P rectangular matrices X and Y and use them to introduce an N × N matrix F given by Next, we assume that the matrix entries of X and Y weigh the double links of a Poisson bipartite graph with probability distribution: where the distribution function (x μ i , y μ i ) controls the correlation between the matrix entries x μ i and y μ i between a pair of connected nodes i and μ. Note that if we take d = N we recover the fully connected or dense limit studied in [3]. If the distribution function ρ is such that the variables x μ i and y μ i are perfectly correlated, one recovers the original Wishart ensemble. Let be the spectrum of a matrix F and recall that its empirical spectral density is given by Then, following [13,20], we recast this as a spin-glass-type problem in which the ρ F (λ) is given by where Z F (z) can be thought of as a partition function of a spin glass model: where we have introduced the following Hamiltonian: with In this context, the empirical spectral density is given by the statistical average of the w-variables via the formula: As, according to equation (11), the empirical spectral density depends on the expectation value on single nodes, it seems natural to apply the cavity method to find a closed set of equations for single-node marginals. To do so, we notice that the Hamiltonian H F (w; z) models interacting variables on an underlying bipartite graph as shown in figure 1. The graph has two types of nodes: factor nodes or μ-nodes, represented by squares and labelled by Greek letters; graph nodes or i-nodes, represented by circles and labelled by Latin letters. Furthermore, on the i-nodes we have dynamical variables, denoted as w = (w 1 , . . . , w N ), while on the factor nodes we have a pair of dynamical variables {m, n} = {(m 1 , . . . , m P ), (n 1 , . . . , n P )}. The graphs nodes are connected to the factor nodes, and vice versa, by pairs of links which correspond to the matrix entries of X and Y. As usual the neighbourhood of a graph node, let us say i, is denoted as ∂i, and similarly, the neighbourhood of a factor node, say μ, is denoted as ∂μ. Finally, given a subset A of nodes, w A ≡ {w } ∈A and dw A ≡ ∈A dw . Before carrying on, it is important to remember that the main pillar the cavity method rests on to be a useful (e.g. exact) method is to look for a physical quantity that becomes uncorrelated after performing an operation on the graph. What that physical quantity is, how it uncorrelates and which operation to perform on the graph depends on the system under study. For instance, in fully connected disordered systems, like the paradigmatic SK model, decorrelation between fields (physical quantity) is achieved via the statistical independence of the couplings if an extra site is removed or added (operation on the system) [21]. In disordered systems on treelike graphs, decorrelation of fields or spins (physical quantity) is achieved via the topological independence induced on the graph after a node has been removed (operation on the graph). If the graph were a tree-like graph at the level of small clusters the operation to perform to achieve topological decorrelation would be to remove those small clusters.
In our case, due the tree-likeness of the underlying graph, the exactness of cavity method rests on topological decorrelation which is achieved by removing single nodes. Note that the density ρ(x, y) contains correlations on the pair of links between two connected nodes, and thus it does not affect the topological decorrelation.
So, after performing the conventional derivation using the cavity method (see, for instance, [13]), we find that single-node marginals P i (w i ) are given by: for all i = 1, . . . , N, where the so-called cavity marginals Q (i) ν (m ν , n ν ) obey the following closed set of equations As it was highlighted in [13], the closed set of equations given by equations (14) and (15) can be further simplified by noticing that the set of Gaussian distributions is a fixed point. With this in mind, and with a modest amount of foresight, we write After plugging equations (16) and (17) into the set of equations (14) and (15), and a bit of algebra, we arrive to the following set of equations Once we have found a solution to the set of equations (18) and (19), numerically or otherwise, the spectral density is given by: with The set of equations (18)- (21) are the first main result of this paper and correspond the empirical spectral density based on the cavity equations. Note that when applying the cavity method to derive these equations one assumes that the underlying graph is a typical large locally tree-like graph, being Poissonian graphs one of several possibilities. For Poissonian graphs the average ensemble of the cavity equations acquires a simpler form. Indeed, following [22] to carry out the ensemble average of the cavity equations, we arrive to the following set of self-consistency equations: where w(g, γ, δ) and v(Δ) are densities. In this case the spectral density is given by instead by the following expression:

Analysis of the theoretical results and comparison with numerical estimates
Let us start first by considering the dense limit consisting on taking the limit d → ∞ while keeping α fixed. There are two ways to do this, both instructive and complementary: either using the set of cavity equations (18) and (19) or the ensemble average equations (22) and (23). We take the first route. Consider for instance the relationship between the coefficients g (i) ν and the Δ (ν) . The Onsager's correction term [19] between g (i) ν and g ν , and similarly for the Δ variables, is order O(d −1 ) and vanishes in the dense limit (see, for instance, [13]). Moreover the coefficients Δ (ν) cannot depend on the {x ν } random variables since for that coefficient the factor node ν has been removed. Thus we can graciously write that We similarly arrive at the following limits This yields the following equations for Δ: from where the spectral density is obtained by ρ(λ) = lim →0 + Im Δ(λ − i )/π, in agreement with the results found in [3]. Notice that in the limit c → 1 this equation reduces to which is the propagator equation for the Marčenko-Pastur spectral density. For the case c = 0 we have instead These cases were thoroughly analysed in [3] and therefore we will not dwell around them any longer.
To obtain the spectral density for the general case of finite d, the set of equations (22) and (23) is solved by using population dynamics, a method that combines Monte Carlo integration and fixed point iteration method. To apply it, each density w(g, γ, δ) and v(Δ) is represented by a population of N random variables {g a , γ a , δ a } N a=1 and {Δ a } N a=1 , respectively, whose histograms are precisely estimates of the aforementioned densities and they become more and more accurate as the population size N increases. More precisely, Then, starting with a random population {g a , γ a , δ a } N a=1 and {Δ a } N a=1 , these are updated according to the following steps: Convergence is usually monitored by using any stopping criterion commonly found in numerical iteration methods.
In figure 2, we compare the theoretical results using population dynamics (represented with markers) with numerical diagonalisation (represented by the histograms that fill each curve to the bottom). For the former we use a population size of N = 10 5 while for the numerics we diagonalise a set of 10 3 matrices of size N = 10 3 . For the comparison we take d = 12 and α ∈ {3, 1, 0.3} (from left to right figures, respectively) and, in each plot, for varying values of the correlation c = E[x μ i y μ i ]. As we can see, the agreement between theory and numerics is excellent. Generally speaking the spectral density does not have hard edges at the end points of its domain. When the rectangularity parameter α > 1, the spectral density ρ(λ) is bounded, while for α 1 it develops a gap and a singularity at λ = 0.

Conclusions
Being random graphs one of the main tools to model complex phenomena, in this paper we have introduced the diluted version of the symmetric cross-correlation matrix ensemble and have obtained its limiting spectral density exactly by using the cavity method. We have checked that our results are correct by thoroughly comparing them to estimates of the empirical density obtained by numerical diagonalization. This work opens the door to study other spectral properties of this ensemble as, for instance, the distribution of extreme eigenvalues, rate functions using large deviation theory, and to explore asymmetric versions of this matrix ensemble. These and some other studies are currently under way.

Data availability statement
No new data were created or analysed in this study.