arXiv:yymm.nnnn

 Hammerstein equations for sparse random matrices

Pawat Akara-pipattana^a and Oleg Evnin^b,c
a Université Paris-Saclay, CNRS, LPTMS, 91405 Orsay, France ^b High Energy Physics Research Unit, Department of Physics, Faculty of Science, Chulalongkorn University, 10330 Bangkok, Thailand ^b Theoretische Natuurkunde, Vrije Universiteit Brussel and
International Solvay Institutes, 1050 Brussels, Belgium
pawat.akarapipattana@universite-paris-saclay.fr, oleg.evnin@gmail.com

ABSTRACT

Finding eigenvalue distributions for a number of sparse random matrix ensembles can be reduced to solving nonlinear integral equations of the Hammerstein type. While a systematic mathematical theory of such equations exists, it has not been previously applied to sparse matrix problems. We close this gap in the literature by showing how one can employ numerical solutions of Hammerstein equations to accurately recover the spectra of adjacency matrices and Laplacians of random graphs. While our treatment focuses on random graphs for concreteness, the methodology has broad applications to more general sparse random matrices.

1 Introduction

The spectra of random matrices have found far-reaching applications in situations where large complex matrices are involved, including fields as diverse as linear [1, 2] and nonlinear [3] dynamics in the presence of random couplings, quantum chaos [4], localization phenomena in disordered media [5], vibrational properties of amorphous solids [6], the physics of heavy nuclei [7], explorations of discrete geometry [8] and mathematical neurobiology [9]. One particular topic is spectral problems related to graphs, and hence to network theory [10], since graphs are connected to matrices via the adjacency matrix representation; see [11] for a textbook exposition.

A number of spectral problems for sparse random matrices, that is, $N\times N$ matrices with a large $N$ and a fixed average number $c$ of nonzero entries per row, can be reduced [13, 12, 14, 15, 16, 17, 18, 19, 20, 21] to solving nonlinear integral equations of the form

$$g(x)=s(x)+\int dx^{\prime}\,K(x,x^{\prime})F(x^{\prime},g(x^{\prime})),$$ (1)

where the ‘source’ $s(x)$, the kernel $K(x,x^{\prime})$ and the nonlinearity $F$ are prescribed, the latter typically being an exponential nonlinearity with respect to the unknown $g(x)$ in the examples we will consider. Such equations are known as Hammerstein equations, after [22]. A typical reaction of a theoretical physicist to an equation like (1) is that of bewilderment. We thus read in the influential paper [23] on sparse matrix spectra: “This integral equation… has, however, so far resisted exhaustive analysis or full numerical solution.” An early numerical treatment can be seen in [16], and it approximates the desired solution via sophisticated sums of Gaussians adapted to the concrete equation at hand, without evoking a systematic theory of equations of the form (1).

Despite having remained unincorporated in the context of sparse random matrix considerations, a systematic theory of Hammerstein equations has been developed in the mathematical literature, starting with [22], see [24] for a comprehensive review. (Curiously, the physical application of such equations mentioned in [24] is in Chandrasekhar’s theory of radiative transfer [25], very far from our present perspective.) In fact, a splash of research activity on numerical convergence properties in concrete cases is seen in the recent years. Our purpose in this article is to demonstrate that this systematic approach to solving Hammerstein equations numerically provides an efficient procedure for recovering sparse matrix spectra in a way that has been previously overlooked in the physics literature. The concrete examples we shall consider are the spectra of adjacency matrices, ordinary graph Laplacians and normalized graph Laplacians of sparse Erdős-Rényi random graphs. These examples are chosen as concrete random matrix ensembles that have received considerable attention in the literature. The technology we discuss should be applicable with minimal modifications to many other sparse random matrices.

The paper is organized as follows: In section 2, we review the emergence of integral equations from the problem of computing the spectrum of adjacency matrices of sparse Erdős-Rényi graphs. In section 3, we describe the process for solving Hammerstein equations numerically. In section 4, we present various tweaks necessary for handling the concrete Hammerstein equations arising from random matrices, as well as the practical results of our numerical work for the spectra of adjacency matrices, ordinary graph Laplacians and normalized graph Laplacians of sparse Erdős-Rényi random graphs. We conclude with a discussion and a wishlist for possible improvements in the numerical algorithms.

2 An integral equation for sparse random graph spectra

We start with a review of the emergence of integral equations of Hammerstein type from the application of statistical field theory methods [26, 27] to spectral problems for sparse random graphs. These methods have been used for a wide range of problems involving random graphs and random matrices [12, 14, 15, 23, 13, 28, 29, 30, 31, 32, 33, 35, 34, 36]. Our purpose here is not to provide a systematic technical exposition of these methods, but rather to give the readers an impression of the sort of structures and steps involved. The readers have the option of fast-forwarding through this section, as the topic of deriving the Hammerstein equations for sparse random matrix spectra is logically independent from the topic of solving them.

The specific variation of the statistical field theory methods we use relies on supersymmetry and auxiliary integrals over anticommuting variables, see [37, 38] for textbook expositions. In application to random matrices, some essential ingredients we use go back to Fyodorov and Mirlin [14, 15]. We have given a pedagogical account of this technology in [19], and refer the readers to that paper for further details.

The example we choose for our demonstration is the eigenvalue density of the adjacency matrix of an Erdős-Rényi random graph. This graph is defined as a set of $N$ vertices with $N\gg 1$, and each pair of vertices is connected randomly and independently with probability $c/N$, so that at $N\to\infty$, each vertex has on average $c$ connections. These data are conveniently coded into the adjacency matrix $\mathbf{A}$, which is a symmetric zero-diagonal matrix whose entry $A_{ij}$ equals 1 if vertices $i$ and $j$ are connected, and 0 otherwise. The above specification of the Erdős-Rényi ensemble is translated into the following probability distribution for the entries of $\mathbf{A}$:

$$P(\mathbf{A})=\prod_{i}\delta(A_{ii})\,\,\prod_{i<j}\left[\left(1-\frac{c}{N}\right)\delta(A_{ij})+\frac{c}{N}\delta(A_{ij}-1)\right]\delta(A_{ji}-A_{ij}),$$ (2)

where $\delta$ stands for the standard $\delta$-function. Finding the eigenvalue density of the adjacency matrices described by this ensemble is a standard problem in random matrix theory. Our treatment follows the original works [15, 14]. We choose this example for its relative simplicity, and also for the fact that it is not directly covered by the pedagogical exposition of [19]. At large $c$, the eigenvalue distribution tends to the semicircle familiar from many incarnations of random matrix theory [39]. At smaller $c$, the shape gets deformed in a complicated way, and it is this deformed shape that will be captured by integral equations we consider in this paper.

We start with the so-called resolvent representation of the eigenvalue density $p(\lambda)$. Let the eigenvalues of $\mathbf{A}$ be $\lambda_{i}$. Then, by definition:

$$p(\lambda)\equiv\frac{1}{N}\sum_{k=1}^{N}\delta(\lambda-\lambda_{k})=\frac{1}{\pi N}\operatorname{Im}\sum_{k}\frac{1}{\lambda-\lambda_{k}-i0}=-\frac{1}{\pi N}\operatorname{Im}\mathrm{Tr}[\mathbf{A}-z\mathbf{I}]^{-1}\Big{|}_{z=\lambda-i0},$$ (3)

where we have used the Sokhotski–Plemelj formula and the spectral decomposition of the matrix inverse, while $z=\lambda-i0$ implies evaluating the relevant expression just below the real axis. The resolvent $\mathrm{Tr}[\mathbf{A}-z\mathbf{I}]^{-1}$ is much simpler to handle within statistical averages than the eigenvalues themselves, since it can be effectively represented using Gaussian integrals, making the dependence on the components of $\mathbf{A}$ simplify further. Note first that, due to the symmetry of the Erdős-Rényi ensemble under vertex permutations, after the averaging, one has

$$\frac{1}{N}\langle\mathrm{Tr}[\mathbf{A}-z\mathbf{I}]^{-1}\rangle\equiv\frac{1}{N}\sum_{k=1}^{N}\langle[\mathbf{A}-z\mathbf{I}]^{-1}_{kk}\rangle=\langle[\mathbf{A}-z\mathbf{I}]^{-1}_{11}\rangle.$$ (4)

The next step is to represent the resolvent (4) through Gaussian integrals, which results in factorization over the entries of $\mathbf{A}$ and hence in a straightforward $\mathbf{A}$-averaging. To do so, it is customary [15, 14, 37, 38] to introduce a 4-component ‘supervector’ at every vertex $k$, given by $\Psi_{k}\equiv(\phi_{k},\chi_{k},\xi_{k},\eta_{k})$, where $\phi_{k}$ and $\chi_{k}$ are ordinary numbers, while $\xi_{k}$ and $\eta_{k}$ are anticommuting numbers satisfying

$$\xi_{i}\xi_{j}=-\xi_{j}\xi_{i},\qquad\xi_{i}^{2}=0,\qquad\int d\xi_{i}\xi_{j}=\delta_{ij},\qquad\int d\xi_{i}=0,$$ (5)

the latter relations known as Berezin integration rules. (All components of $\xi$ anticommute with all components of $\eta$ as well. Note that Berezin integrals should not be visualized as any generalization of Riemann sums. They are linear maps from functions of anticommuting numbers to real numbers that have been developed historically for the purpose of convenient path-integral representation of fermionic systems.) We furthermore introduce a ‘scalar product’ of supervectors defined as

$$\Psi^{\dagger}_{k}\Psi_{l}\equiv\phi_{k}\phi_{l}+\chi_{k}\chi_{l}+\frac{\xi_{k}\eta_{l}-\eta_{k}\xi_{l}}{2},\qquad(\Psi^{\dagger}_{k}\Psi_{k}\equiv\phi_{k}^{2}+\chi_{k}^{2}+\xi_{k}\eta_{k}).$$ (6)

With these preliminaries, we can write

$$(\mathbf{A}-z\mathbf{I})^{-1}_{11}=\int d\mathbf{\Psi}\,(2i\phi_{1}^{2})\,e^{i\sum_{kl}A_{kl}\Psi^{\dagger}_{k}\Psi_{l}-iz\sum_{k}\Psi^{\dagger}_{k}\Psi_{k}}.$$ (7)

Note that, if substituting (6) into (7), one obtains independent Gaussian integrals involving $\exp(i\sum A_{kl}\phi_{k}\phi_{l}-iz\sum\phi_{k}^{2})$, $\exp(i\sum A_{kl}\chi_{k}\chi_{l}-iz\sum\chi_{k}^{2})$ and $\exp(i\sum A_{kl}\xi_{k}\eta_{l}-iz\sum\xi_{k}\eta_{k})$. Absorbing the powers of $\pi$ arising from the standard Gaussian integrals over $\phi$ and $\chi$ into the definition of the measure, the first integral evaluates to $(\mathbf{A}-z\mathbf{I})^{-1}_{11}/\sqrt{\det(\mathbf{A}-z\mathbf{I})}$ due to the insertion of $\phi_{1}^{2}$, the second integral evaluates to $1/\sqrt{\det(\mathbf{A}-z\mathbf{I})}$, while the last integral evaluates to $\det(\mathbf{A}-z\mathbf{I})$ by the standard formula for Gaussian integrals in Berezin integration theory. As a result, the powers of determinants exactly cancel, leaving behind (7).

The purpose of inserting integrals over anticommuting variables in (7) was precisely to make the powers of determinant arising from Gaussian integration exactly cancel out, yielding an expression factorized over the components of $\mathbf{A}$. An alternative strategy, often seen in statistical physics literature, is the replica method [12, 13, 23, 35, 40], where one introduces an arbitrary controllable number of extra Gaussian integrals over ordinary commuting variables, formally extrapolating the total number of Gaussian integrals to 0 in the final formulas at the end of the computation. While, in practical terms, this strategy is largely interchangeable with the supersymmetry-based approach we are employing here, the latter dispenses of the need for the ambiguous operation of extrapolating a positive integer to 0.

Averaging (7) over the random matrix ensemble (2) is straightforward and yields

$$\left<(\mathbf{A}-z\mathbf{I})^{-1}_{11}\right> =\int d\mathbf{\Psi}\,(2i\phi_{1}^{2})\,e^{-iz\sum_{k}\Psi^{\dagger}_{k}\Psi_{k}}\prod_{k<l}\left[\left(1-\frac{c}{N}\right)+\frac{c}{N}e^{2i\Psi^{\dagger}_{k}\Psi_{l}}\right]$$

$$=\int d\mathbf{\Psi}\,(2i\phi_{1}^{2})\,e^{-\frac{1}{2N}\sum_{kl}C(\Psi_{k},\Psi_{l})-iz\sum_{k}\Psi^{\dagger}_{k}\Psi_{k}},$$ (8)

where we have used $c/N\ll 1$ and then introduced

$$C(\Psi,\Psi^{\prime})=c\left(1-e^{2i\Psi^{\dagger}_{k}\Psi_{l}}\right).$$ (9)

To solve the supervector model (8) at large $N$, we employ, after Fyodorov and Mirlin [14, 15, 19], the following Gaussian functional integral over an auxiliary function $g(\Psi)$:

$$\int\mathcal{D}g\exp\left[-\frac{N}{2}\int d\Psi\,d\Psi^{\prime}\,g(\Psi)\,C^{-1}(\Psi,\Psi^{\prime})\,g(\Psi^{\prime})+i\sum_{k}g(\Psi_{k})\right]=e^{-\frac{1}{2N}\sum_{kl}C(\Psi_{k},\Psi_{l})},$$ (10)

where $C^{-1}$ is the inverse of $C$ in the sense of integral convolution:

$$\int d\Psi\,C(\Psi_{1},\Psi)\,C^{-1}(\Psi,\Psi_{2})=\delta(\Psi_{1}-\Psi_{2}).$$ (11)

The explicit form of $C^{-1}$ will not be necessary for our derivation. Once this transformation has been implemented, the integrals over $\Psi_{k}$ with different $k$ completely factorize and are performed independently, yielding an explicit saddle point problem:

$$\begin{split}&\langle(\mathbf{A}-z\mathbf{I})^{-1}_{11}\rangle=\int\mathcal{D}g\,e^{-N\hskip 0.28453ptS[g]}\,\frac{\int d\Psi\,(2i\phi^{2})\,e^{ig(\Psi)-iz\Psi^{\dagger}\Psi}}{\int d\Psi\,e^{ig(\Psi)-iz\Psi^{\dagger}\Psi}},\\ &S[g]\equiv\frac{1}{2}\int d\Psi\,d\Psi^{\prime}\,g(\Psi)\,C^{-1}(\Psi,\Psi^{\prime})\,g(\Psi^{\prime})-\ln\left(\int d\Psi\,e^{ig(\Psi)-iz\Psi^{\dagger}\Psi}\right).\end{split}$$ (12)

At large $N$, the integral is dominated by the stationary points of the functional $S[g]$ defined by $\delta S/\delta g(\Psi)=0$, which is

$$g(\Psi)=i\,\frac{\int d\Psi^{\prime}\,C(\Psi,\Psi^{\prime})\,e^{ig(\Psi^{\prime})-iz\Psi^{\prime\dagger}\Psi^{\prime}}}{\int d\Psi^{\prime}\,e^{ig(\Psi^{\prime})-iz\Psi^{\prime\dagger}\Psi^{\prime}}}.$$ (13)

This equation is closely reminiscent of the analogous Bray-Rodgers equation in the replica approach [13, 12, 23].

The final step is to assume that the dominant saddle point described by (13) is invariant under ‘superrotations’ that preserve the scalar product $\Psi^{\dagger}\Psi$ (an alternative would have been continuous families of saddle points related by superrotations, since the equation is symmetric under applying a superrotation to the argument of $g$). This implies that $g(\Psi)=g(\Psi^{\dagger}\Psi)$, and with this simplification, the integral equation (13) can be reduced to a single one-dimensional integral equation. To do so, we express the commuting components of $\Psi$ in polar coordinates as

$$\rho=\phi^{2}+\chi^{2},\qquad\phi=\sqrt{\rho}\cos\alpha,\qquad\chi=\sqrt{\rho}\sin\alpha,\qquad d\phi\,d\chi=\frac{1}{2}d\rho\,d\alpha.$$ (14)

With this,

$$\Psi^{\dagger}\Psi=\rho+\xi\eta,\qquad g(\Psi^{\dagger}\Psi)=g(\rho)+\partial_{\rho}g(\rho)\,\xi\eta,$$ (15)

where the Taylor expansion of $g(\rho+\xi\eta)$ truncates at the $\xi\eta$-term, since all higher powers of anticommuting variables vanish by (5). Inserting this representation in the $\Psi^{\prime}$-integral in (13) and performing the integration over all variables except for $\rho^{\prime}$, we obtain, after some integration by parts in $\rho^{\prime}$:

$$g(\rho)=ic\rho\int_{0}^{\infty}\!\!\!d\rho^{\prime}\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}})}{\sqrt{\rho\rho^{\prime}}}\,e^{ig(\rho^{\prime})-iz\rho^{\prime}},$$ (16)

where $J_{1}$ is the usual Bessel function. This equation, which dates back to [12] where it was derived by replica methods, and its analogs for graph Laplacians, will be the main subjects of our numerical work in what follows. Similar integral equations with exponential nonlinearities and Bessel kernels arise also in other sparse random matrix problems.

Once a solution of (16) has been found, one has to substitute it in (12) and perform a saddle point evaluation of the integral. This, in general, would involve complicated determinants from Gaussian integrations, but in our case, these determinants must cancel out. An elementary way to see it is as follows: consider, tautologically, (8) with $(\mathbf{A}-z\mathbf{I})^{-1}_{11}$ replaced by 1 (so that the average is guaranteed to equal 1). Then, our entire derivation will go through, arriving at (12), but without the insertion ${\int d\Psi\,(2i\phi^{2})\,e^{\{\cdots\}}}/{\int d\Psi\,e^{\{\cdots\}}}$. Neither the saddle point equation (16) nor the determinants arising from Gaussian integrations depend on the presence of this insertion. On the other hand, by our tautological starting point, the result without the insertion must equal 1. The result with the insertion, as in (12), must then be equal to the insertion evaluated at the saddle point, that is,

$$\langle(\mathbf{A}-z\mathbf{I})^{-1}_{11}\rangle=\frac{\int d\Psi\,(2i\phi^{2})\,e^{ig(\Psi^{\dagger}\Psi)-iz\Psi^{\dagger}\Psi}}{\int d\Psi\,e^{ig(\Psi^{\dagger}\Psi)-iz\Psi^{\dagger}\Psi}},$$ (17)

where $g$ is now understood as the solution of (16). The integrals can again be simplified using (14-15) and assuming that $g(\rho)$ decays at infinity, and then plugged back into (3-4) to yield the following result for the eigenvalue density:

$$p(\lambda)=\frac{1}{\pi}\,\operatorname{Re}\int_{0}^{\infty}\!\!\!d\rho\,e^{ig(\rho)-iz\rho}\Big{|}_{z=\lambda}.$$ (18)

By construction, $p(\lambda)$ is normalized as a probability density $\int d\lambda\,p(\lambda)=1$. It is rather striking, at the face value, that the solution of the formidable equation (16) processed with (18) is automatically normalized in this way. It is nonetheless correct, as follows from our arguments, and as will be verified by numerical simulations. (In [19], we ignored all normalization at the intermediate stages of computation, since it is always easy to normalize the end result. Here, however, we restore all normalizations for comparisons with numerics, which can be done, as it turns out, with minimal effort.)

3 Handling Hammerstein equations

Equation (16) evidently matches the general structure of Hammerstein equations given by (1). Before proceeding with its solution, we will review the general underlying principles, following [24].

In the random matrix literature, the systematic theory of Hammerstein equations appears to have never been applied. Indeed, in [23], the apparent intractability of integral equations like (16) is used to motivate an alternative approach, where one goes back to an analog of (13) and then develops a sampling scheme in the spirit of ‘population dynamics’ to approximate its solutions. In [16], equations analogous to (16) are approached directly, and their solutions are approximated by elaborate sums of Gaussians. This approach is not apparently connected, however, to the strategies for solving Hammerstein equations seen in the mathematical literature.

What the mathematical literature instructs us to do is to take (1) and approximate its solutions by functions within some finite-dimensional linear subspace spanned by $\phi_{j}(x)$ with $j=1..J$, that is $g(x)=\sum_{j}\beta_{j}\phi_{j}(x)$. One can then replace the original equation (1) by the following modification

$$\sum_{j}\beta_{j}\phi_{j}(x)=\mathbb{P}\left[s(x)+\int dx^{\prime}\,K(x,x^{\prime})F(x^{\prime},{\textstyle\sum_{j}}\beta_{j}\phi_{j}(x^{\prime}))\right],$$ (19)

where $\mathbb{P}$ is a projector on the subspace spanned by $\phi_{j}(x)$. There are different possible choices for this projector, leading to different numerical solution methods, as we shall discuss shortly. One can view (19) as a system of $J$ nonlinear algebraic equations for $J$ unknowns $\beta_{j}$ that can be solved by any standard numerical methods. One expects that, as the dimension $J$ of the subspace spanned by $\phi_{j}(x)$ increases, the resulting $g(x)=\sum_{j}\beta_{j}\phi_{j}(x)$ will approximate the true solution of (1) better and better. The details of this convergence have been studied by mathematicians, but we shall be pragmatic about it since rather small values of $J$ around 10 will be sufficient to reproduce sparse matrix spectra with a good precision.

There is a technical issue with equation (19) as it stands. Because $\beta_{j}$ appear inside the integral over $x^{\prime}$ nonlinearly, the integral will have to be recomputed for each assignment of $\beta_{j}$, for example, when solving the algebraic system (19) by means of iterative methods. These repeated integrations are numerically costly. Kumar and Sloan proposed in [41] a simple trick to bypass this issue. Instead of $g(x)$, one introduces

$$f(x)\equiv F(x,g(x)),$$ (20)

so that (1) is equivalently rewritten as

$$f(x)=F\left(x,s(x)+\int dx^{\prime}\,K(x,x^{\prime})f(x^{\prime})\right).$$ (21)

We then expand $f(x)$, rather than $g(x)$, in the truncated basis $\phi_{j}(x)$ with $j=1..J$:

$$f(x)=\sum_{j=1}^{J}\beta_{j}\phi_{j}(x),$$ (22)

and project (21), rather than (1), on this basis using a suitable projector $\mathbb{P}$:

$$\sum_{j}\beta_{j}\phi_{j}(x)=\mathbb{P}F\left(x,s(x)+{\textstyle\sum_{j}}\beta_{j}\int dx^{\prime}\,K(x,x^{\prime})\phi_{j}(x^{\prime})\right).$$ (23)

This is, again, solved for $\beta_{j}$ as a system of nonlinear algebraic equations. An advantage over (19) is that the integrals $\int dx^{\prime}\,K(x,x^{\prime})\phi_{j}(x^{\prime})$ are computed once and for all, and are not affected by the values of $\beta_{j}$. After $\beta_{j}$ has been found in this way, providing an approximation to $f(x)$, $g$ is best reconstructed by applying one more integration:

$$g(x)=s(x)+\int dx^{\prime}\,K(x,x^{\prime})f(x^{\prime}).$$ (24)

Regarding the choice of the projector $\mathbb{P}$, a natural first thought is to use an orthogonal projector with a suitably defined inner product on the space of functions. This choice leads to Galerkin (or spectral) methods for solving Hammerstein equations. These methods are extensively covered in the literature, but a more economical setup in the context of (23) is given by the collocation method that uses the so-called interpolating projector [42]. For a given input $h(x)$, the interpolating projector provides $\mathbb{P}h(x)$ as a function of the form $\sum_{j}\beta_{j}\phi_{j}(x)$ that equals $h(x)$ at prescribed collocation points $x_{k}$ with $k=1..J$. This provides exactly $J$ linear conditions for $J$ unknowns $\beta_{j}$. With such a prescription for the projector $\mathbb{P}$, (21) becomes

$$\sum_{j}\beta_{j}\phi_{j}(x_{k})=F\left(x_{k},s(x_{k})+{\textstyle\sum_{j}}C_{kj}\beta_{j}\right),\qquad C_{kj}\equiv\int dx\,K(x_{k},x)\phi_{j}(x).$$ (25)

Note that the $J\times J$ matrix $C_{kj}$ only has to be evaluated once. With this matrix computed, the $J$ equations for $\beta_{j}$, indexed by $k$, are solved using the standard numerical routines for multidimensional root search.

One still has to decide on how to choose the basis $\phi_{j}$ and the collocation points $x_{k}$. Mathematical literature consistently advises us to work with orthogonal polynomials. In [42], integral equations on the interval $[-1,1]$ are considered, and then it is natural to choose $\phi_{j}$ as the Legendre polynomials. Some considerations of half-infinite integration ranges, as in (16), can be seen in [43]. Since the integrals are from $0$ to $\infty$, it is natural to use Laguerre polynomials but we will make some different choices in the concrete examples below, depending on the concrete form of the integrand appearing in the Hammerstein equation. Note that since no orthogonal projection is involved in the collocation equation (25), the orthogonality property of the polynomial family is never used explicitly, and any other family of linearly independent functions could be used, at least in principle. Using orthogonal polynomials is, however, a way to ascertain a ‘good degree’ of linear independence. For example we have checked that constructing $\phi_{j}$ from plain polynomials $1$, $x$, $x^{2}$, etc without orthogonalization results in a poorly conditioned system, and the numerical search for the solutions of (25) does not reliably converge. Orthogonal polynomials work well, by contrast. In [43], we are instructed that a good choice of collocation points is the roots of the lowest-degree orthogonal polynomial not included in the chosen set $\phi_{j}$. We will generally follow this prescription, though according to our simulations reported below, reasonable variations in the positions of the collocation points do not immediately compromise the operation of the algorithm.

4 Numerical results for sparse Erdős-Rényi graphs

We proceed to report how numerical solutions of Hammerstein equations via the collocation method (25) work for determining eigenvalue spectra of adjacency matrices and Laplacians of sparse random graphs. For each concrete equation at hand, we will have to indicate how to effectively partition the functional dependences in the integral convolutions on the right-hand side to match the notation in (1), and how to choose the sets of trial functions $\phi_{j}$ and collocation points.

4.1 Adjacency matrices

Before applying the technology of section 3 to equation (16), it is wise to implement the following rescalings:

$$\rho=\tilde{\rho}/\sqrt{c},\hskip 14.22636ptz=\tilde{z}\sqrt{c},\hskip 14.22636ptg(\rho(\tilde{\rho}))=\tilde{g}(\tilde{\rho}).$$ (26)

Dropping the tildes for convenience, we obtain the rescaled equation

$$g(\rho)=i\rho\int_{0}^{\infty}d\rho^{\prime}\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}/c})}{\sqrt{\rho\rho^{\prime}/c\,}}\,e^{ig(\rho^{\prime})-iz\rho^{\prime}}.$$ (27)

An advantage of this representation is that the relevant range of $\rho$ becomes approximately independent of $c$. Indeed, at large $c$, we can expand $J_{1}(2x)=x+\cdots$ and the equation turns into $g(\rho)=i\rho\int d\rho^{\prime}\,e^{ig(\rho^{\prime})-iz\rho^{\prime}}$, solved by a fixed-size Wigner semicircle [19]. At smaller $c$, the functional dependences evidently get deformed, but they occupy roughly the same region with respect to the new $\rho$ variable, rather than getting scaled $\sim\sqrt{c}$ with respect to the old $\rho$ variable. This makes the numerical implementation more neat.

There is some ambiguity in how to split the various $\rho$-dependences in (27) to match the structure of (1). Evidently, $s(\rho)=0$. We find it convenient to identify

$$K(\rho,\rho^{\prime})=i\rho\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}/c})}{\sqrt{\rho\rho^{\prime}/c\,}},\qquad F(\rho,g(\rho))=e^{ig(\rho)-iz\rho}.$$ (28)

An advantage of including the $z$-dependence into the nonlinearity $F$ and not into the kernel $K$ is that, as we vary $z$ to derive the eigenvalue density at different points, the integrals $C_{kj}$ in the Kumar-Sloan collocation problem (25) will not have to be recomputed. For the basis $\phi_{j}$ we choose Laguerre polynomials multiplied by exponentials. Namely, we write

$$e^{ig(\rho)-iz\rho}=\sum_{j=0}^{J}\beta_{j}L_{j}(\gamma\rho)\,e^{-\gamma\rho},$$ (29)

where $L_{j}$ are Laguerre polynomials orthogonal with respect to $\int_{0}^{\infty}dx\,e^{-x}L_{i}(x)L_{j}(x)$. We will treat $\gamma$ as a tunable real parameters and adjust it to optimize the numerical performance. We remark that we could have used different scalings for the arguments of the Laguerre polynomials and the exponential, which gives more control over the numerical performance, but the practical gain is not significant and the formulas become more bulky. We give a summary of that scheme in Appendix A, and proceed here with (29) as given.

We can evaluate the integral in (27) analytically using the identity [44]

$$\sum_{n=0}^{\infty}\frac{w^{n}\,L_{n}^{(k)}(x)}{(n+k)!}=e^{w}(xw)^{-k/2}J_{k}(2\sqrt{xw}),$$ (30)

where $L^{(k)}$ are the associated Laguerre polynomials, and choosing $x=\gamma\rho^{\prime}$, $w=\rho/\gamma c$ and $k=1$. From (27), (29) and (30), together with $L^{(1)}_{n}=\sum_{m=0}^{n}L_{m}$, we obtain

$$g(\rho) =i\rho\,e^{-\rho/\gamma c}\sum_{n=0}^{\infty}\frac{(\rho/\gamma c)^{n}}{\gamma(n+1)!}\sum_{m=0}^{\min(n,J)}\beta_{m}$$

$$=i\rho\,e^{-\rho/\gamma c}\left[\sum_{m=0}^{J}\beta_{m}\sum_{n=0}^{\infty}\frac{(\rho/\gamma c)^{n}}{\gamma(n+1)!}-\sum_{n=0}^{J-1}\frac{(\rho/\gamma c)^{n}}{\gamma(n+1)!}\sum_{m=n+1}^{J}\beta_{m}\right]$$

$$=ic\left[(1-e^{-\rho/\gamma c})\sum_{m=0}^{J}\beta_{m}-e^{-\rho/\gamma c}\sum_{n=1}^{J}\frac{\rho^{n}}{n!(\gamma c)^{n}}\sum_{m=n}^{J}\beta_{m}\right]$$

$$=ic\sum_{m=0}^{J}\beta_{m}\left(1-e^{-\rho/\gamma c}\sum_{n=0}^{m}\frac{(\rho/\gamma c)^{n}}{n!}\right).$$ (31)

Note that the $\rho^{\prime}$-integral is completely gone due to the orthogonality of Laguerre polynomials. Recombined with (29), this generally results in the following problem:

$$\sum_{j=0}^{J}\beta_{j}L_{j}(x)=\mathbb{P}\exp\left[x-\frac{izx}{\gamma}-c\sum_{j=0}^{J}\beta_{j}\left(1-e^{-x/\gamma^{2}c}\sum_{n=0}^{j}\frac{(x/\gamma^{2}c)^{n}}{n!}\right)\right],$$ (32)

where we defined $x\equiv\gamma\rho$, and $\mathbb{P}$ is the projector of choice on the space of degree $J$ polynomials in terms of $x$. This should be treated as a system of $J+1$ equations, from identifying the polynomial coefficients on the two sides, for $J+1$ unknowns $\beta_{j}$. By orthogonality of Laguerre polynomials, (18) then yields for the following eigenvalue density, keeping in mind that we must undo the rescalings (26):

$$p_{\mathrm{adj}}(\lambda)=\frac{\operatorname{Re}\beta_{0}\big{|}_{z=\lambda/\sqrt{c}}}{\gamma\pi\sqrt{c}}.$$ (33)

For the simplest choice of $\mathbb{P}$ we take the interpolating projector, equating the polynomial and the exponential in (32) at the collocation points $x_{k}$ with $k=0..J$, chosen as the roots of $L_{J+1}$:

$$\sum_{j=0}^{J}\beta_{j}L_{j}(x_{k})=\exp\left[x_{k}-\frac{izx_{k}}{\gamma}-c\sum_{j=0}^{J}\beta_{j}\left(1-e^{-x_{k}/\gamma^{2}c}\sum_{n=0}^{j}\frac{(x_{k}/\gamma^{2}c)^{n}}{n!}\right)\right].$$ (34)

Figure 1: Solutions of the collocation problem (34) for Erdős-Rényi graph adjacency matrices at $J=10$, $\gamma=1$, converted to eigenvalue distribution estimates using (33) and plotted as solid black lines for (left) $c=15$ and (right) $c=8$. As the distributions are reflection-symmetric, only the $\lambda>0$ part is plotted explicitly. The grey shaded areas represent empirical eigenvalue density histograms obtained from a sample of 1000 Erdős-Rényi graphs with 10000 vertices each at the corresponding values of $c$. Tiny undulations are seen near the center of the $c=8$ plot. These numerical artifacts become more problematic at smaller values of $c$.

To check the performance, we program this in Python (a sample script is provided in Appendix B) using the root function of the SciPy library [45] to solve (34) numerically with $\gamma=1$ and the initial seed for the root search taken as $\beta_{0}=1$, $\beta_{j>0}=0$ at the first $z$-point, with the previous solutions reused as the seeds for each subsequent $z$-point.

The results of our simulations are presented in Fig. 1. It can be seen that the performance is perfect at $c=15$, which is shared by any higher values of $c$. At $c=8$ tiny undulations are seen near the center. These undulations become much more dramatic at lower $c$, disrupting the precision of the algorithm performance, though the overall shape of the distribution remains qualitatively correct. We shall comment in the conclusions on possible strategies for stabilizing the numerical performance at lower values of $c$.

4.2 Ordinary graph Laplacians

Graph Laplacians are of central importance for transport phenomena on networks. We start with the ordinary graph Laplacian that controls, in particular, random walks on graphs where there is a fixed transition rate per edge per unit time [19]. In terms of the adjacency matrix $\mathbf{A}$, drawn for Erdős-Rényi graphs from the distribution (2), this Laplacian is defined as follows: we first construct the diagonal degree matrix $\mathbf{D}$ whose diagonal entries are $D_{ii}\equiv\sum_{j}A_{ij}$, and then the ordinary Laplacian $\mathbf{L}$ is defined by

$$\mathbf{L}=\mathbf{D}-\mathbf{A}.$$ (35)

Unlike the adjacency matrix, this matrix is positive semi-definite.

An analytic theory of the ordinary Laplacian along the lines of section 2 was developed in [19]. One has to solve the Hammerstein equation

$$g(\rho)=-ic\left\{e^{i\rho}-1-\rho\,e^{i\rho}\int_{0}^{\infty}d\rho^{\prime}\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}})}{\sqrt{\rho\rho^{\prime}}}e^{ig(\rho^{\prime})-i(z-1)\rho^{\prime}}\right\}$$ (36)

that replaces (16) derived for adjacency matrices. The eigenvalue density of (35) is recovered from the solution as

$$p_{L}(\lambda)=\frac{1}{\pi}\,\operatorname{Re}\int_{0}^{\infty}\!\!\!d\rho\,e^{ig(\rho)-iz\rho}\Big{|}_{z=\lambda}=\frac{1}{\pi c}\operatorname{Im}[\partial_{\rho}g(\rho,z)]\Big{|}_{\rho=0,\,\,z=\lambda+1},$$ (37)

where we have restored, based on the principles outlined in section 2, the exact probability density normalization, ignored in the analytic formulas of [19]. To get to the last representation, one uses the $\rho$-derivative of (36) at $\rho=0$. Equations of the form (36) date back to [13] where the derivations were based on replica methods. They have also been derived, together with equations like (16) for adjacency matrices, in [17] where the analysis based on the methods of moments. The latter strategy has the advantage of being mathematically rigorous, but requires a genuine tour de force of combinatorial accounting and resummations.

As in the case of adjacency matrices, it is convenient to implement a $c$-dependent rescaling to make the relevant ranges of variables depend only weakly on $c$. We choose this rescaling as

$$z=c+1+{\tilde{z}}{\sqrt{c}},\qquad\rho=\frac{\tilde{\rho}}{\sqrt{c}},\qquad g(\rho(\tilde{\rho}),z(\tilde{z}))=\tilde{g}(\tilde{\rho},\tilde{z})+\tilde{\rho}\sqrt{c}.$$ (38)

The resulting equation, with tildes dropped, is

$$g(\rho)=-ic\left(e^{i\rho/\sqrt{c}}-1-\frac{i\rho}{\sqrt{c}}\right)+i\rho\,e^{i\rho/\sqrt{c}}\int_{0}^{\infty}d\rho^{\prime}\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}/c})}{\sqrt{\rho\rho^{\prime}/c\,}}e^{ig(\rho^{\prime})-iz\rho^{\prime}}.$$ (39)

At large $c$, this equation reduces to $g(\rho)=i\rho^{2}/2+i\rho\int d\rho^{\prime}\,e^{ig(\rho^{\prime})-iz\rho^{\prime}}$ [19, 46, 47], which describes a ‘free convolution of the Wigner semicircle and a Gaussian’ in the language of free probability theory [48]. To represent (39) in the form (1), we choose

$$s(\rho)=-ic\left(e^{i\rho/\sqrt{c}}-1-\frac{i\rho}{\sqrt{c}}\right),\quad K(\rho,\rho^{\prime})=i\rho\,e^{i\rho/\sqrt{c}}\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}/c})}{\sqrt{\rho\rho^{\prime}/c\,}},\quad F(\rho,g(\rho))=e^{ig(\rho)-iz\rho}.$$

We then follow the Kumar-Sloan strategy for implementing collocation described in section 3, and specifically apply the expansion (29).

Since the integral over $\rho^{\prime}$ in (39) is exactly the same as in (27), we can rely on the same strategy as under (30) to explicitly evaluate all the integrals under the ansatz (29). This yields, instead of (31),

$$g(\rho)=-ic\left(e^{i\rho/\sqrt{c}}-1-\frac{i\rho}{\sqrt{c}}\right)+ice^{i\rho/\sqrt{c}}\sum_{m=0}^{J}\beta_{m}\left(1-e^{-\rho/\gamma c}\sum_{n=0}^{m}\frac{(\rho/\gamma c)^{n}}{n!}\right),$$ (40)

and leads to the following collocation problem:

$$\sum_{j=0}^{J}\beta_{j}L_{j}(x_{k})=\exp\Bigg{[} x_{k}-\frac{izx_{k}}{\gamma}+c\left(e^{ix_{k}/\gamma\sqrt{c}}-1-\frac{ix_{k}}{\gamma\sqrt{c}}\right)$$ (41)

$$-ce^{ix_{k}/\gamma\sqrt{c}}\sum_{m=0}^{J}\beta_{m}\left(1-e^{-x_{k}/\gamma^{2}c}\sum_{n=0}^{m}\frac{(x_{k}/\gamma^{2}c)^{n}}{n!}\right)\Bigg{]},$$

where $x_{k}$ are the roots of $L_{J+1}$. From the solution of this algebraic system, one recovers via (37), (38) and (40)

$$p_{L}(c+1+{\tilde{z}}{\sqrt{c}})=\frac{1}{\pi\sqrt{c}}\operatorname{Re}\int_{0}^{\infty}\!\!\!d\tilde{\rho}\,e^{i\tilde{g}(\tilde{\rho})+i\tilde{\rho}\sqrt{c}-i(c+1+{\tilde{z}}{\sqrt{c}})\tilde{\rho}/\sqrt{c}}=\frac{1}{\gamma\pi\sqrt{c}}\operatorname{Re}\int_{0}^{\infty}\!\!\!dx$$

$$\hskip 42.67912pt\times\sum_{j=0}^{J}\beta_{j}L_{j}(x)\,e^{-(1+i/\gamma\sqrt{c})x}=\frac{1}{\gamma\pi\sqrt{c}}\operatorname{Re}\int_{0}^{\infty}\!\!\!dx\,\sum_{j=0}^{J}\beta_{j}\sum_{n=0}^{j}{j\choose n}\frac{(-x)^{n}}{n!}\,e^{-(1+i/\gamma\sqrt{c})x}$$

$$\hskip 42.67912pt=\frac{1}{\gamma\pi\sqrt{c}}\operatorname{Re}\sum_{j=0}^{J}\beta_{j}\sum_{n=0}^{j}{j\choose n}\frac{(-1)^{n}}{(1+i/\gamma\sqrt{c})^{n+1}}=\frac{1}{\pi}\operatorname{Im}\sum_{j=0}^{J}\frac{\beta_{j}(\tilde{z})}{\left(1-i\gamma\sqrt{c}\right)^{j+1}}.$$ (42)

We have performed numerical simulations based on (41) and (42) and observed that they reproduce very well the empirical distribution at $c\gtrsim 8$, though increasing the expansion parameter $\gamma$ becomes necessary towards the lower end of this range to ensure numerical convergence. At smaller values of $c$, the performance is less stable, especially in regions away from the main peak, where spurious oscillations develop. The performance at $c=8$ is seen on the left panel of Fig. 2. It captures the smooth descending part of the distribution perfectly, but is less accurate to the left of the main peak, where rough textures drvelop as $c$ tends to the ‘percolation threshold’ at $c=1$, even though the overall shape is reproduced correctly.

Figure 2: Results for the eigenvalue density of the ordinary graph Laplacian plotted as solid black lines: (left) at $c=8$ using the Laguerre collocation problem (41) with $J=10$, $\gamma=5$, from which the eigenvalue density is extracted via (42), (right) at $c=4$ using the general Kumar-Sloan collocation (25) with the half-range Hermite expansion (43) at $J=10$, followed by the integration in (37). Additionally, on the right panel, we add a dotted curve representing the result of solving (41) with $J=10$, $\gamma=5$ at $c=4$, showing more spurious undulations than the half-range Hermite method. The grey shaded areas represent empirical eigenvalue density histograms obtained from a sample of 1000 Erdős-Rényi graphs with 10000 vertices each at the corresponding values of $c$.

One can understand intuitively the deteriorating performance of the expansion scheme (29) away from the main peak. Indeed, there, the first derivative of $g(\rho)$ at $\rho=0$ becomes small and $e^{ig-iz\rho}$ is dominated by the Gaussian envelope $e^{-\rho^{2}/2}$, clearly visible [19] in the large $c$ limit of ($\ref{HammLscale})$. If we expand as in (29), this Gaussian envelope is approximated by a polynomial, which is clearly prone to instabilities due to the large $\rho$ growth of polynomials and large $\rho$ decay of $e^{-\rho^{2}/2}$. To remedy for this problem, we can replace (29) by an alternative expansion

$$e^{ig(\rho)-iz\rho}=\sum_{j=0}^{J}\beta_{j}H_{j}(\rho)\,e^{-\rho^{2}/2},$$ (43)

where $H_{j}$ are the half-range Hermite polynomials [49] orthogonal with respect to the inner product $\int_{0}^{\infty}d\rho\,e^{-\rho^{2}}H_{i}(\rho)H_{j}(\rho)$. The price to pay is that there are no longer simple and nice analytic formulas for $H_{j}$, but we can still process all the integrals numerically and run the generic formulation of Kumar-Sloan collocation given by (25). We have observed that this approach performs better at smaller $c$. (A sample script is given in Appendix C.) In the right panel of Fig. 2, the performance of this scheme at the low value $c=4$ is shown. The smooth heavy tail of the distribution is captured very accurately, while the rough part around the ascent and the peak is reproduces less well, though the overall shape is correct.

At larger values of $c$, both approaches quickly become very accurate, and the empirical distribution shape is effectively indistinguishable from our numerical solutions of the Hammerstein equation (36) already at c=15.

4.3 Normalized graph Laplacians

The normalized Laplacian is defined, in the notation of (35), by

$$\boldsymbol{\mathcal{L}}\equiv\mathbf{D}^{-1/2}(\mathbf{D}-\mathbf{A})\mathbf{D}^{-1/2},$$ (44)

This Laplacian controls, in particular, random walks where there is a fixed rate per unit time to leave the present vertex and jump with equal probability to one of its neighbors [19]. (In the above defintion, $\mathbf{D}^{-1/2}$ is understood as the square root of the pseudoinverse of $\mathbf{D}$.)

The spectral theory of normalized Laplacians of Erdős-Rényi graphs was developed in [19]. It results in the following Hammerstein equation:

$$g(\rho)=-ic\left(e^{i(1-z)\rho}-1-\rho\,e^{i(1-z)\rho}\int_{0}^{\infty}d\rho^{\prime}\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}})}{\sqrt{\rho\rho^{\prime}}}\,e^{ig(\rho^{\prime})+i(1-z)\rho^{\prime}}\right).$$ (45)

(Curiously, the same equation was recently found to control the probabilities of first return times of random walks on Erdős-Rényi graphs [21].) From the solution of this equation, the eigenvalue density can be recovered in the following form, with the normalization given explicitly:

$$p_{\mathcal{L}}(\lambda)=\frac{1}{\pi}\operatorname{Re}\int_{0}^{\infty}d\rho\,[c+ig(\rho)]\,e^{ig(\rho)}\Big{|}_{z=\lambda}.$$ (46)

As in our previous treatments, it is convenient to introduce rescalings that minimize the dependence of the contributing ranges of variables on $c$. This is accomplished with

$$z=1+\frac{\tilde{z}\sqrt{c}}{c+1},\qquad\rho=\frac{\tilde{\rho}}{\sqrt{c}},\qquad g(\rho(\tilde{\rho}),z(\tilde{z}))=\tilde{g}(\tilde{\rho},\tilde{z})-\tilde{z}\tilde{\rho}\frac{c}{c+1}.$$ (47)

Dropping all tildes, we then get

$$g(\rho)=-ic\left(e^{\frac{-iz\rho}{c+1}}-1+\frac{iz\rho}{c+1}\right)+i\rho\,e^{\frac{-iz\rho}{c+1}}\int_{0}^{\infty}d\rho^{\prime}\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}/c})}{\sqrt{\rho\rho^{\prime}/c}}\,e^{ig(\rho^{\prime})-iz\rho^{\prime}}.$$ (48)

The $c\to\infty$ limit $g(\rho)=i\rho\int d\rho^{\prime}\,e^{ig(\rho^{\prime})-iz\rho^{\prime}}$ is the same as for adjacency matrices, defining a Wigner semicircle [19]. For the Kumar-Sloan collocation, we identify $F(\rho,g(\rho))=e^{ig(\rho)-iz\rho}$ and

$$s(\rho)=-ic\left(e^{\frac{-iz\rho}{c+1}}-1+\frac{iz\rho}{c+1}\right),\quad K(\rho,\rho^{\prime})=i\rho\,e^{\frac{-iz\rho}{c+1}}\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}/c})}{\sqrt{\rho\rho^{\prime}/c\,}}.$$ (49)

The integral over $\rho^{\prime}$ in (48) is once again completely identical to (27), so we adopt the expansion (29) and apply the same processing as in section 4.1 to obtain

$$g(\rho)= -ic\left(e^{\frac{-iz\rho}{c+1}}-1+\frac{iz\rho}{c+1}\right)+ice^{\frac{-iz\rho}{c+1}}\sum_{m=0}^{J}\beta_{m}\left(1-e^{-\rho/\gamma c}\sum_{n=0}^{m}\frac{(\rho/\gamma c)^{n}}{n!}\right).$$ (50)

The collocation problem, analogous to (34) for adjacency matrices, then becomes

$$\sum_{j=0}^{J}\beta_{j}L_{j}(x_{k}) =\exp\Bigg{[}x_{k}-\frac{izx_{k}}{\gamma}+c\Bigg{(}e^{\frac{-izx_{k}}{\gamma(c+1)}}-1+\frac{izx_{k}}{\gamma(c+1)}\Bigg{)}$$ (51)

$$\hskip 85.35826pt-ce^{\frac{-izx_{k}}{\gamma(c+1)}}\sum_{m=0}^{J}\beta_{m}\Bigg{(}1-e^{-x_{k}/\gamma^{2}c}\sum_{n=0}^{m}\frac{(x_{k}/\gamma^{2}c)^{n}}{n!}\Bigg{)}\Bigg{]},$$

where the collocation points $x_{k}$ are the roots of $L_{J+1}$. After these equations have been solved numerically, one must reconstruct $e^{ig}$ from (29), $g$ from (50), undo the rescaling (47) and evaluate (46). In this way, we obtain:

$$p_{\mathcal{L}}(\lambda)=\frac{\sqrt{c}}{\pi\gamma}\operatorname{Re}\int_{0}^{\infty}dx\,\Bigg{[}1-\sum_{m=0}^{J}\beta_{m}\Big{(}1-e^{-x/\gamma^{2}c}\sum_{n=0}^{m}\frac{(x/\gamma^{2}c)^{n}}{n!}\Big{)}\Bigg{]}\sum_{j=0}^{J}\beta_{j}L_{j}(x)\,e^{-x}$$ (52)

$$=\frac{\sqrt{c}}{\pi\gamma}\operatorname{Re}\Bigg{[}\beta_{0}\Big{(}1-\sum_{m=0}^{J}\beta_{m}\Big{)}+\sum_{j,m=0}^{J}\beta_{m}\beta_{j}\sum_{n=0}^{m}\sum_{k=0}^{j}{j\choose k}\frac{(-1)^{k}}{n!k!(\gamma^{2}c)^{n}}\int_{0}^{\infty}dx\,x^{n+k}e^{-x(1+1/\gamma^{2}c)}\Bigg{]}$$

$$=\frac{\sqrt{c}}{\pi\gamma}\operatorname{Re}\Bigg{[}\beta_{0}\Big{(}1-\sum_{m=0}^{J}\beta_{m}\Big{)}-\sum_{j,m=0}^{J}\beta_{m}\beta_{j}\sum_{n=0}^{m}\sum_{k=0}^{j}{j\choose k}{n+k\choose k}\frac{(-\gamma^{2}c)^{k+1}}{(\gamma^{2}c+1)^{n+k+1}}\Bigg{]}\Bigg{|}_{z=\frac{(\lambda-1)(c+1)}{\sqrt{c}}}.$$

Figure 3: Solutions of the collocation problem (51) for the normalized Laplacians of Erdős-Rényi graphs at $J=9$, $\gamma=2$, converted to eigenvalue distribution estimates using (52) and plotted as solid black lines for (left) $c=8$ and (right) $c=4$. As the distributions are reflection-symmetric, only the $\lambda>1$ part is plotted explicitly. We excise small regions near $\lambda=1$ where the convergence of our collocation scheme is compromised by the presence of sharp spikes in the actual eigenvalue distribution. (There are generically $\delta$-function spikes at $\lambda=1$ that become more and more prominent at small $c$. We comment on them further in the conclusions.) The grey shaded areas represent empirical eigenvalue density histograms obtained from a sample of 1000 Erdős-Rényi graphs with 10000 vertices each at the corresponding values of $c$.

We display the output of solving (51), followed by an application of (52) in Fig. 3. The performance of this algorithm remains stable and successful down to rather low values $c=8$ and 4. Convergence is more problematic near the origin at low $c$, where a large abrupt spike develops, hence we excise these regions. In the bulk of the distribution, the numerics based on (51) and (52) gives very convincing results.

5 Outlook

We have applied a systematic treatment of nonlinear integral equations of the Hammerstein type to recover the spectra of sparse random matrices associated to Erdős-Rényi graphs with an asymptotically large number of vertices of average degree $c$. Our treatment is practically exact for $c$ larger than 15 or so, and remains rather accurate down to $c$ around 8. As $c$ is lowered further, approaching the percolation threshold at $c=1$, challenges in numerical performance are met, though for graph Laplacians, we still obtain very reasonable output for $c$ as low as 4. We hope that further improvements in numerical methods for solving Hammerstein equations will lead to better performance in the low-$c$ region.

The approach we have taken here is very complementary to the widely explored ‘cavity’ and ‘population dynamics’ methods for reconstructing sparse graph spectra [23, 50, 40]. In those approaches, stochastic simulations are designed whose stationary distributions contain information on eigenvalue densities. By contrast, our numerics is fully deterministic and consists in solving the integral equations (16), (36) and (45) using collocation methods. The implementation is very economical in its computational costs: numerical solutions of the Hammerstein equations in all of our examples (the solid black lines in the plots) take seconds to obtain on an ordinary personal computer. We feel that it would be very interesting to apply this technology to more sophisticated sparse matrix problems, for example, those arising in the theory of amorphous solids [6] or for non-Hermitian matrices [36, 51, 52].

The main challenges met by our numerical algorithms are either in the regions where the eigenvalue density becomes very small, or at low values of $c$, where it is known from empirical diagonalization of large matrices that the regularity of the eigenvalue density deteriorates dramatically. We provide below preliminary comments on the origin of these issues and potential pathways to their solution.

The poor performance outside the main support of the eigenvalue density can be easily understood qualitatively. The eigenvalue density is related to the imaginary part of the $\rho$-derivative of $g$ at the origin, as in (37). When the eigenvalue density is small, the imaginary part of $g$ grows slowly with $\rho$, and hence $e^{ig}$ decreases slowly with $\rho$. As a result, a large range of $\rho^{\prime}$ becomes relevant in the integrals on the right-hand side of the Hammerstein equations, with complicated oscillatory behaviors, and one should not expect that such integrals will be easily approximated using interpolating polynomials. The situation is significantly better for ordinary Laplacians, where $e^{ig}$ includes a Gaussian envelope $e^{-\rho^{2}/2}$, independent of the eigenvalue density, and this tames the integrals over $\rho^{\prime}$ when we use the half-range Hermite polynomial expansion adapted to this Gaussian envelope. (We note that the other cases also have similar Gaussian envelopes, but the width is of order $\sqrt{c}$ rather than of order 1. This may provide a viable alternative strategy at small $c$, again, moving away from simple Laguerre expansions.) We feel that a general scheme to handle the equations outside the main support of the eigenvalue density should exist, possibly by rearranging the integrals in the complex $\rho$-plane, but we have not been able to formulate a working recipe. This problem is of somewhat marginal importance for the eigenvalue distribution topic, since one is mostly interested in the shape of the eigenvalue distributions where they are nonvanishing, and not in the regions where they (almost) vanish. By contrast, for applications of Hammerstein equations to random walks on random graphs, as in [21], the behavior of $g$ near $z=\infty$, very far away from the main support of the eigenvalue density, is of primary importance.

Using Laguerre polynomial expansions, we have been able to evaluate all the integrals in the Hammerstein equations with Bessel kernels exactly, obtaining explicit systems of algebraic equations (34), (41) and (51). This both provides for very fast numerical performance and creates a comfortable environment for more detailed mathematical analysis of convergence, approximation accuracy, etc. Where it is advantageous to use a different kind of expansion, as in our half-range Hermite treatment of the ordinary graph Laplacians, the integrals can still be effectively handled numerically following the general Kumar-Sloan scheme (25).

As to improving the performance of our numerics at lower $c$, the issue seems to be that the shape of various functional dependences becomes more complicated at lower $c$, and our simple polynomial expansions fail with capturing it. One could hope that increasing the polynomial order $J$ in (29) would help, but it is not straightforward for two reasons. First, as we go to higher $J$, we would need to solve nonlinear systems of equations similar to (34) in higher and higher numbers of dimension, where naive root searches become less and less reliable. Perhaps this could be tackled by first running the algorithm at a lower value of $J$ and then reusing the output as seeds for the root search at higher values of $J$. More dramatically, however, the interpolating projector we have used with great efficiency at low $J$ is likely to be flawed at higher $J$. It is well-known that naively increasing the order of interpolation does not necessarily result in better approximation [53]. There could certainly be more refined choices for the projector $\mathbb{P}$ in (32) and other similar equations. A naive option is to rely on the infinite radius of convergence of the Taylor series for exponentials and employ a Taylor expansion on the right-hand side of (32). This, however, would require working at a very high arithmetic precision, since the convergence of Taylor expansions for exponentials requires massive numerical cancellations. There likely exist better choices for $\mathbb{P}$ that we leave for future investigations, also hoping for practical input from mathematicians specializing in numerical methods. Uniform approximations of functions by polynomials on infinite intervals in the presence of exponential suppression have been discussed systematically in the literature [54].

In our numerical work, we have treated the expansion parameters such as $\gamma$ in (29) and the initial seed for multidimensional root searches as free parameters chosen empirically to ensure satisfactory performance. In the cases we considered, simple straightforward choices sufficed to obtain the plots given in the paper. It could be good to develop a better approach where these parameters are bootstrapped toward good values, either by analyzing the integral equations at a given value of $z$, or by reusing the previously found solutions at other values of $z$ (the latter is actually the approach we took for the initial seed).

Nonlinear algebraic equations like (34), (41) and (51) are generally expected to have multiple solutions. Only one of these solutions corresponds to the actual eigenvalue density curve, while the others are spurious. We could see that deforming the parameters of our expansion schemes (for example, making $\gamma$ small) or the initial seed far from the optimal values we found leads to such spurious convergence, predicting a completely wrong eigenvalue density curve. For the concrete cases we considered, we indicated simple natural choices of the numerical implementation parameters that ensure convergence to the correct solution. It would be nice to have more of a theory of the solutions of (34), (41) and (51) and their basins of attraction relative to root search algorithm, especially given that the equations are provided explicitly, with all the integrals evaluated in a closed form. Uniqueness of solutions of Hammerstein equations (before converting them to approximate collocation schemes) has been discussed extensively in the mathematical literature.

One striking feature of the eigenvalue distributions of sparse matrices at lower values of $c$ is the emergence of sharp peaks [23, 55]. These peaks are understood qualitatively as coming from small disconnected components. For example, isolated vertices of degree 0 contribute zero eigenvalues to the adjacency matrix. Since there are on average $e^{-c}N$ vertices of degree 0, one expects a $e^{-c}\delta(\lambda)$ peak in the eigenvalue distribution. Similarly, graph Laplacians develop zero eigenvalues for each disconnected component of the graph. There are further peaks at nonzero values of $\lambda$. We suspect that they come from graph components that are only connected to the rest of the graph by one edge, and other peculiar arrangements (similar structures play a significant role in the distribution of resistance distances [34]). These sharp peaks are, incidentally, reflected differently in the empirical histograms and in the numerical solutions of the integral equations appearing in our plots. Numerical histograms, in the limit of large samples, always show the integrals of the eigenvalue density over single-bin intervals. The numerical solutions are evaluated at a single $\lambda$-point, except that the truncated approximations of the original integral equation, as in (29), are likely to resolve the infinitely thin $\delta$-function peaks somewhat.

The sharp peaks could be of significant interest physically, depending on the setting. For example, where random matrix spectra represent vibrational modes, the peaks in the eigenvalue density will signify strongly resonant points in the spectrum. The $\delta$-function singularities in the eigenvalue distribution should translate to singular behaviors of $g(z;\rho=0)$ at some special values of $z$. It would be very nice to develop a systematic theory of these sharp peaks, and use this analytic understanding to subtract that part within the equations of $g$, so that the target functions for numerical work become more smooth and the numerical methods become more stable.

Appendix A A more general expansion in terms of Laguerre polynomials

Instead of (29), we could have employed a more general expansion

$$e^{ig(\rho)-iz\rho}=\sum_{j=0}^{J}\beta_{j}L_{j}(a\rho)\,e^{-b\rho},$$ (53)

where the scalings of the argument of the Laguerre polynomial and the exponential are different. This does not immediately yield significant improvements over the case $a=b=\gamma$ treated in the main text, and the formulas become more bulky. These formulas, however, may become useful for other future applications, so we give a summary here on how to treat the adjacency matrix equation (27) with the expansion (53).

The integrals can in general be evaluated using the relation of the Bessel convolutions – that is, Hankel transforms – to 2d Fourier transforms. Specifically, we write $J_{1}(\rho)=-\partial_{\rho}J_{0}(\rho)=-\pi^{-1}\partial_{\rho}\int_{0}^{\pi}d\varphi\,e^{-i\rho\cos\varphi}$. Then, for any $f(\rho)$ that decays at infinity,

$$\rho\int_{0}^{\infty}d\rho^{\prime}\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}/c})}{\sqrt{\rho\rho^{\prime}/c\,}}f(\rho^{\prime})=-\frac{c}{\pi}\int_{0}^{\infty}d\rho^{\prime}\int_{0}^{\pi}d\varphi\,\partial_{\rho^{\prime}}e^{-2i\sqrt{\rho\rho^{\prime}/c}\cos\varphi}f(\rho^{\prime})=cf(0)$$ (54)

$$+\frac{c}{\pi}\int_{0}^{\infty}d\rho^{\prime}\int_{0}^{\pi}d\varphi\,e^{-2i\sqrt{\rho\rho^{\prime}/c}\cos\varphi}\partial_{\rho^{\prime}}f(\rho^{\prime})=\frac{2c}{\pi}\int_{-\infty}^{\infty}dx\int_{0}^{\infty}dy\,e^{-2ix\sqrt{\rho/c}}\partial_{\rho^{\prime}}f(\rho^{\prime})\Big{|}_{\rho^{\prime}=x^{2}+y^{2}}+cf(0),$$

where we have introduced the Cartesian coordinates $x\equiv\sqrt{\rho^{\prime}}\cos\varphi$, $y\equiv\sqrt{\rho^{\prime}}\sin\varphi$.

If we start with the definition of Laguerre polynomials

$$L_{j}(a\rho)e^{-b\rho}=e^{-b\rho}\sum_{n=0}^{j}{j\choose n}\frac{(-a\rho)^{n}}{n!},$$ (55)

substitute $\rho=x^{2}+y^{2}$ and expand everything in terms of monomials of the form $x^{2k}y^{2l}$, the integrals are evaluated term-by-term using the two elementary Gaussian quadratures:

$$\int_{-\infty}^{\infty}x^{2n}e^{2ix\sqrt{\rho/c}}e^{-bx^{2}}dx=e^{-\rho/bc}\int_{-\infty}^{\infty}(\tilde{x}+i\sqrt{\rho/c})^{2n}e^{-b\tilde{x}^{2}}d\tilde{x}$$ (56)

$$=e^{-\rho/bc}\sum_{k=0}^{n}{2n\choose 2k}\left(\frac{-\rho}{c}\right)^{\!\!n-k}\int_{-\infty}^{\infty}\tilde{x}^{2k}e^{-b\tilde{x}^{2}}d\tilde{x}=e^{-\rho/bc}\sqrt{\frac{\pi}{b}}\sum_{k=0}^{n}{2n\choose 2k}\left(\frac{-\rho}{c}\right)^{\!\!n-k}\frac{(2k-1)!!}{(2b)^{k}},$$

and

$$\int_{0}^{\infty}y^{2m}e^{-by^{2}}dy=\sqrt{\frac{\pi}{b}}\,\frac{(2m-1)!!}{2(2b)^{m}},$$ (57)

where one must define $(-1)!!=1$ to make the formulas work correctly. Then,

$$\int_{-\infty}^{\infty}dx\int_{0}^{\infty}dy\,(x^{2}+y^{2})^{n}e^{2ix\sqrt{\rho/c}}e^{-b(x^{2}+y^{2})}=\sum_{m=0}^{n}{n\choose m}\!\!\int_{-\infty}^{\infty}x^{2m}e^{2ix\sqrt{\rho/c}}e^{-bx^{2}}dx\!\int_{0}^{\infty}y^{2(n-m)}e^{-by^{2}}dy$$

$$\hskip 56.9055pt=\frac{\pi}{2b}e^{-\rho/bc}C^{\mathrm{aux}}_{n}(\rho),$$ (58)

where we have introduced the following auxiliary polynomials

$$C^{\mathrm{aux}}_{n}(\rho)\equiv\frac{1}{(2b)^{n}}\sum_{m=0}^{n}{n\choose m}(2(n-m)-1)!!\sum_{k=0}^{m}{2m\choose 2k}(2k-1)!!\left(\frac{-2b\rho}{c}\right)^{\!\!m-k}.$$ (59)

With these ingredients, we can write for the integrals of the individual terms in (55):

$$\rho\int_{0}^{\infty}d\rho^{\prime}\,\frac{J_{1}(2\sqrt{\rho\rho^{\prime}/c})}{\sqrt{\rho\rho^{\prime}/c\,}}\,\rho^{\prime n}e^{-b\rho^{\prime}}=\frac{2c}{\pi}\!\int_{-\infty}^{\infty}dx\!\int_{0}^{\infty}dy\,e^{2ix\sqrt{\rho/c}}\,[n(x^{2}+y^{2})^{n-1}\!-\!b(x^{2}+y^{2})^{n}]e^{-b(x^{2}+y^{2})}$$

$$\hskip 28.45274pt=ce^{-\rho/bc}\Big{[}\frac{n}{b}C^{\mathrm{aux}}_{n-1}(\rho)-C^{\mathrm{aux}}_{n}(\rho)\Big{]}.$$ (60)

To make this formula work correctly at $n=0$, we must add $c$ due to the last term in (54). If we now impose (27) at collocation points $\rho=\rho_{l}$ with $l=0..J$ under the ansatz (53) and use the above formulas to evaluate the Bessel integral, keeping in mind that $C^{\mathrm{aux}}_{0}(\rho)=1$, the output is the following collocation problem

$$\sum_{j=0}^{J}L_{jl}\beta_{j}=\exp\Bigg{[}(b-iz)\rho_{l}+c(e^{-\rho_{l}/bc}-1)\sum_{j=0}^{J}\beta_{j}+\sum_{j=1}^{J}C_{jl}\beta_{j}\Bigg{]},$$ (61)

with $\displaystyle L_{jl}\equiv\sum_{n=0}^{j}{j\choose n}\frac{(-a\rho_{l})^{n}}{n!}$ and

$$C_{jl}\equiv ce^{-\rho_{l}/bc}\sum_{n=1}^{j}{j\choose n}\frac{(-a)^{n}}{n!}\Big{[}C^{\mathrm{aux}}_{n}(\rho_{l})-\frac{n}{b}C^{\mathrm{aux}}_{n-1}(\rho_{l})\Big{]},$$ (62)

with $C^{\mathrm{aux}}$ given by (59). While the multiple sums may seem awkward, they merely encode the coefficients of explicit polynomials in $\rho$ and can be treated very effectively in numerical implementations. They also need to be evaluated only once. Note that no approximations have been made in (61-62) beyond truncating the expansion (53) at $j=J$.

For the collocation points we choose $\rho_{l}=x_{l}/a$. where $x_{l}$ are the roots $L_{J+1}$, so that $L_{J+1}(x_{l})=0$. We then use the above collocation problem to find $\beta_{j}$ and hence $e^{ig(\rho)-iz\rho}$, and then the eigenvalue density can be recovered from (18), which gives

$$p_{\mathrm{adj}}(\lambda)=\frac{1}{\pi b}\,\operatorname{Re}\sum_{j=0}^{J}\beta_{j}\sum_{n=0}^{j}{j\choose n}(-{a}/{b})^{n}=\frac{1}{\pi b}\sum_{j=0}^{J}(1-{a}/{b})^{j}\,\operatorname{Re}\beta_{j}\Big{|}_{z=\lambda/\sqrt{c}}.$$ (63)

Appendix B Python code for the adjacency spectrum

We provide a basic Python script that recovers the solution of the collocation problem (34) for adjacency matrices. Related collocation problems (41) and (51) for graph Laplacians can be treated by minimal modifications of this script.

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

def eqs(beta,z):
 betac=beta[:J]+1j*beta[J:]
 eqs=[sum([L[k,j]*betac[j] for j in range(J)])-np.exp(x[k]*(1-1j*z/gamma)
      -sum([C[k,j]*betac[j] for j in range(J)])) for k in range(J)]
 return np.append(np.real(eqs),np.imag(eqs))

c=15
sqrtc=np.sqrt(c)
J=11  #corresponds to J+1 in the notation of the paper
gamma=1
zs=np.linspace(0.01,2.1,100)
x=np.polynomial.laguerre.lagroots([0]*J+[1])

fac=np.ones((J))
for i in range(1,J-1):
 fac[i+1]=(i+1)*fac[i]

C=np.zeros((J,J))
L=np.zeros((J,J))
for k in range(J):
 xgc=x[k]/(gamma**2*c)
 expxgc=np.exp(-xgc)
 for j in range(J):
  L[k,j]=np.polynomial.laguerre.Laguerre([0]*j+[1])(x[k])
  C[k,j]=c*(1-expxgc*sum([xgc**n/fac[n] for n in range(j+1)]))

p=[]
beta=np.array([1]+[0]*(2*J-1))
for z in zs:
 sol=opt.root(eqs,beta,args=(z),tol=1e-12)
 beta=sol.x
 p.append(beta[0]/(gamma*np.pi*sqrtc))

plt.plot(zs*sqrtc,p)
plt.ylim(0)
plt.show()

Appendix C Python code for the ordinary Laplacian via half-range Hermite polynomials

We provide a simple Python script that implements the half-range Hermite collocation for ordinary graph Laplacians. While we could evaluate some of the integrals here analytically in principle, the numerical evaluation works well for practical purposes.

import numpy as np
from scipy.optimize import root
from scipy.integrate import quad
from scipy.special import jv
import matplotlib.pyplot as plt

def hrHermites(n):
 polys=[np.poly1d([(np.sqrt(np.pi)/2)**(-0.5)])]
 for i in range(1,n+1):
  this_poly=np.poly1d([1]+[0]*i)
  for last_poly in polys:
   this_poly-=last_poly*quad(lambda y: this_poly(y)*last_poly(y)*np.exp(-y**2),
                                             0,np.inf)[0]
  this_poly=this_poly/np.sqrt(quad(lambda y: this_poly(y)**2*np.exp(-y**2),
                                             0,np.inf)[0])
  polys.append(this_poly)
 return polys

def eqs(beta,z):
 betac=beta[:J]+1j*beta[J:]
 eqs=[]
 for k in range(J):
  xg=1j*x[k]/sqrtc
  expxg=np.exp(xg)
  eqs.append(sum([H[k,j]*betac[j] for j in range(J)])-np.exp(x[k]**2/2-1j*x[k]*z
      +c*(expxg-1-xg)-expxg*sum([C[k,j]*betac[j] for j in range(J)])))
 return np.append(np.real(eqs),np.imag(eqs))

c=5
sqrtc=np.sqrt(c)
J=12
zs=np.linspace(-3,4,100)
hrH=hrHermites(J)
x=hrH[-1].r

C=np.zeros((J,J))
H=np.zeros((J,J))
for k in range(J):
 for j in range(J):
  H[k,j]=hrH[j](x[k])
  C[k,j]=sqrtc*quad(lambda y:jv(1,2*np.sqrt(x[k]*y/c))*np.sqrt(x[k]/y)
                                       *hrH[j](y)*np.exp(-y**2/2),0,np.inf)[0]

p=[]
beta=np.zeros(2*J)
for z in zs:
 sol=root(eqs,beta,args=(z),tol=1e-12)
 beta=sol.x
 betac=beta[:J]+1j*beta[J:]
 f_real=lambda y: np.real(np.exp(-1j*y/np.sqrt(c)-y**2/2)*np.sum([betac[j]
                      *hrH[j](y) for j in range(J)]))
 f_imag=lambda y: np.imag(np.exp(-1j*y/np.sqrt(c)-y**2/2)*np.sum([betac[j]
                      *hrH[j](y) for j in range(J)]))
 p.append(np.real(quad(f_real,0,np.inf)[0]
                     +1j*quad(f_imag,0,np.inf)[0])/(np.pi*sqrtc))

plt.plot(c+1+zs*sqrtc,p)
plt.xlim(0)
plt.ylim(0)
plt.show()