Kac-Rice fixed point analysis for single- and multi-layered complex systems

Introduction

The stability of large complex systems has been of growing interest in the scientific literature ever since Robert May famously asked “Will a large complex be stable?” ; he considered criteria for a high-dimensional random linear model to be stable. At first sight, linearity appears strange in this context, since complex systems are almost always considered to be non-linear. In fact, non-linearity is typically taken as a defining characteristic of complex systems. The idea behind May’s model is to imagine the linear system as a leading order approximation near a fixed point, thus enabling stability analysis of the fixed point freed from the complications of non-linearity. Its beautiful simplicity is to a large extend what makes the linear model useful. Nonetheless, understanding the effect of non-linearity on the stability of large complex systems remains an outstanding open problem of high value.

Consider a general discrete-time dynamical system

where the map $\mathbf{f}$ is assumed to be highly non-linear. One of the simplest (yet interesting) questions we can ask about a dynamical system (1) is for the number of fixed points, i.e. solutions $\mathbf{x}_{*}$ to the system $\mathbf{x}_{*}=\mathbf{f}(\mathbf{x}_{*})$. Determining the fixed points would be the first step in a stability analysis of such systems.

In this paper, we will ask for number of fixed points of a general class of single- and multi-layered systems. The systems that we will consider are chosen to satisfy certain symmetry constraints but are otherwise chosen at random. In this way, we may view such systems as a null model for large complex systems. An analysis of mean number of fixed points in a continuous-time dynamical system was initiated in , but our model differ from the one considered in in several important aspects. Most notable, our system evolve in discrete-time and include the possibility of multi-layer constructions.

The paper is organised as follows. In the Section 2, we introduce the model that we are going to study and present a theorem (Theorem 2.1) which links the mean number of fixed points to a problem within Random Matrix Theory. The theorem is proved in the appendix. In Section 3, we exploit techniques from Random Matrix Theory (in particularly recently developed results for products of random matrices) to find the asymptotic behaviour of the mean number of fixed points for large (i.e. high dimensional) systems. We end the paper with a summary and brief discussion of some open problems in Section 4. The paper has been written in such a way so Section 2 as well as Appendix A can be read without any prior knowledge about Random Matrix Theory, while Section 3 requires no prior knowledge about random fields.

Our model and main results

In certain applications, it is natural to decompose the map $\mathbf{f}$ appearing in (1) into a sequence of iterated sub-maps $\mathbf{f}=\mathbf{f}^{(1)}\circ\mathbf{f}^{(2)}\circ\cdots\circ\mathbf{f}^{(D)}$ so that

In such cases, we say the system is multi-layered with depth $D$ and we refer to $\mathbf{f}^{(d)}$ ($d=1,\ldots,D$) as the $d$-th layer. For $D=1$, we say that the system is single-layered.

Since the layers are zero-mean Gaussians, their structure is completely determined by their (matrix-valued) correlation kernels,

It often useful also to have entry-wise notation, in which case we write

with $\mathbf{f}^{(d)}=(f^{(d)}_{i})_{i=1,\ldots,{N_{d-1}}}$ and $\mathbf{K}=(K_{ij})_{i,j=1,\ldots,N_{d-1}}$. In this paper, we will furthermore assume that the kernel has the form

It almost goes without saying that homogeneity corresponds to the (stochastic) symmetry of translation invariance, while domain- and codomain-isotropy correspond to rotation invariance (including parity inversions) in the domain and codomain, respectively. We note that there is no distinctions between (stochastic) symmetries of the wide or strict sense, since we are considering Gaussian maps. It is common in the probability literature to refer to the domain and codomain of a random function as time and space, respectively. Consequently, the stochastic symmetries described above are referred to as stationarity, time- and space-isotropy rather than homogeneity, domain- and codomain-isotropy (see e.g. ). However, our model (2) already contains a notion of time so such terminology is inappropriate in the present context. Moreover, to us, the notion of multi-dimensional time appears rather contrived, thus we prefer the above given terminology.

and has fast decay at infinity. These conditions are sufficient to ensure that we can choose the sample layers $\mathbf{f}^{(d)}$ to be regular enough for our purposes.

Under the above given regularity assumptions, our system (2) always has at least one fixed point (almost surely), since each layer $\mathbf{f}^{(d)}$, and consequently the iterated map $\mathbf{f}=\mathbf{f}^{(1)}\circ\mathbf{f}^{(2)}\circ\cdots\circ\mathbf{f}^{(D)}$, has zero-mean and fast decaying correlations for increasing $\lVert\mathbf{x}-\mathbf{y}\rVert$. Our main result (stated below and proven in Appendix A) gives a formula for the mean number of fixed points.

Let $\mathbf{J}_{d}$ be an $N_{d-1}\times N_{d}$ random matrix whose entries are i.i.d. centred Gaussian random variables with variance $\sigma_{d}^{2}$, and let $\mathbf{J}_{1},\ldots,\mathbf{J}_{D}$ be stochastically independent. Consider the multi-layered random dynamical system of depth $D$ as defined above with $\sigma_{d}:=(-\kappa_{d}^{\prime}(0))^{1/2}>0$ and denote by $\mathcal{N}^{(D)}_{\mathbf{f}}$ an integer-valued random variable which gives the number of fixed points. Then, we have

where the expectation on the right-hand side is with respect to the joint distribution of $\mathbf{J}_{1},\ldots,\mathbf{J}_{D}$, i.e.

As a prelude to the more involved analysis performed in Section 3, let us study the simplest possible case, namely $N=D=1$. This case requires no prior knowledge about techniques from Random Matrix Theory. However, we emphasise that while the $N=D=1$ problem may be trivial from the perspective of Theorem 2.1, it is a non-trivial problem in the sense that we are asking for solutions to a system

where $f$ is a highly non-linear function. Nonetheless, it follows immediately from Theorem 2.1 with $N=D=1$ that we have

with $\sigma=\sigma_{1}$. Asymptotically the mean number of fixed points (10) behave as

where we have used standard asymptotic notation in which $f\sim g$ means $f/g\to 1$.

Asymptotic behaviour for large systems

In the end of previous section, we saw that the one-dimensional single-layered system ($N=D=1$) has a plateau for small $\sigma$ where the mean number of fixed points is approximately equal to one, but that the number starts to increase for $\sigma\gtrapprox 1$, cf. Figure 1 (left panel). In this section, we will see that this behaviour is only intensified as $N$ grows larger. In fact, we will argue below that in the large-$N$ limit the system develops a third-order phase transition which separates a region with a single fixed point and a region with large the number of fixed points.

To analyse the large-$N$ behaviour of our system, we will use Theorem 2.1 together with techniques from Random Matrix Theory. First, we note that the expectation on the right-hand side of (7) only depends on the product matrix

where each $\mathbf{J}_{d}$ is a (rectangular) random matrix with i.i.d. centred Gaussian entries. In fact, due to invariance of the determinant under similarity transformations, only the eigenvalues of the product matrix $\mathbf{X}_{D}$ matters. We note that the product matrix $\mathbf{X}_{D}$ is a real matrix, thus its eigenvalues are either real or complex conjugate pairs. Let us assume that $\mathbf{X}_{D}$ has exactly $n$ real eigenvalues denoted $\lambda_{1},\ldots,\lambda_{n}$ and $m=(N-n)/2$ complex conjugate pairs denoted $z_{1},z_{1}^{*},\ldots,z_{m},z_{m}^{*}$, then we make the trivial (but nonetheless important) observation that

We note that $n$ must have the same parity as $N$ (i.e. $n\equiv N\mod 2$), since the complex eigenvalues are paired. The form (13) is important, since statistical properties of the eigenvalues of products of independent real Gaussian matrices have been studied rather extensively over the last five years ; see for an overview of recent progress on products of random matrices. A remarkable result is that the eigenvalues of the product matrix $\mathbf{X}_{D}$ belongs to a certain class of Pfaffian point processes . We state this result more precisely through the following proposition.

Three further comments regarding Proposition 3.1 are in order before we proceed. First, by setting $\nu_{d}\geq 0$ we have implicitly assumed that $N=N_{0}=N_{D}$ is the smallest matrix dimension. Going back to our multi-layered system (2), we see that this assumption can be made without loss of generality, since it is merely a question about on which manifold we are looking for fixed points. On the random matrix side, this is equivalent to the cyclic invariance of the absolute determinant,

and we note that our problem does not depend on the variances $\sigma_{1},\ldots,\sigma_{D}$ individually, but only on their product. Consequently, it is useful to introduce the geometric mean

as the important parameter of our problem.

Third, we emphasise that we intentionally refer to (14) as the “joint density function” but not the “joint probability density function”, since it is normalised to $p_{N,n}^{(D)}$ rather than unity. Due to this fact, it also useful to introduce the partial expectation with respect to joint density function (14) defined as

Here, we have used (13), (18), and (19). We recall that terms on the right-hand side for which $n$ and $N$ have different parity is equal to zero by definition.

Important quantities for our purposes are the (mean) spectral density for real and complex eigenvalues given by

It worth noting that the first sum which appear in the spectral density of the real eigenvalues (22) starts at one rather than zero, because matrices with no real eigenvalues (obviously) do not contribute to the spectral density of real eigenvalues.

An essential property the spectral densities (22) is that under proper rescaling they concentrate their mass on regions with compact support. This is often referred to as the global (or macroscopic) scaling regime and we state the known result as proposition.

Using notation as above, with $\nu_{1},\ldots,\nu_{D-1}$ and $D$ kept fixed, we have

which may be considered a corollary of Proposition 3.2. In other words, in the large-$N$ limit the fraction of real eigenvalues tends to zero, such that the spectrum is completely dominated by the complex eigenvalues. The limit (26) can be obtained using techniques from free probability , while the limit (25) is more challenging, since the real spectrum is subdominant. The spectral form (25) was first conjectured in while a proof (based on explicit formulae derived in ) was provided in . We refer the aforementioned papers for the details of the proof.

We now return to our original question, namely to determine the mean number of fixed points for large $N$. The leading-order behaviour is captured by the so-called ‘complexity’ defined as

It is evident that by using the complexity (28) we presuppose that the number of fixed points grows exponentially with $N$. We will verify below that such exponential growth is indeed present beyond a certain threshold. From Theorem 2.1, we have

where $\mathbf{X}_{D}$ denotes the product matrix (12). As emphasised earlier, the right-hand side of (29) depends only on the eigenvalues of the product matrix. Furthermore, for large $N$ the contribution from the complex eigenvalues dominate, which allows us to write

A change of variables, $z\mapsto\hat{z}=z/N^{D/2}$, yields

which by comparison with Proposition 3.2 tells us that the natural scale of our problem is set by

In the limit $N\to\infty$ the complexity is therefore given by the integral

Moreover, the global density (26) is rotational invariant in the complex plane, thus we can perform the above integral using the identity

Thus, changing to polar coordinates in (33) and performing the angular part of the integral using (34), we obtain a simple expression for the complexity

We note that similar method was used in to find the complexity for a certain type of neural networks. In fact, the matrix expectation considered in corresponds to our single layer case ($D=1$).

The complexity (28) in our problem plays a role similar to that of the free energy in equilibrium statistical mechanics. We note that the complexity (35) is twice continuously differentiable, but its third derivative $C^{\prime\prime\prime}(\hat{\sigma})$ is discontinuous at $\hat{\sigma}=1$. Thus, by analogy to conventions from statistical mechanics we shall say that our system develops a third order phase transition at the critical value $\hat{\sigma}_{c}=1$ in the large-$N$ limit. Links between the spectral edge behaviour in Random Matrices and third order phase transitions in physical systems have recently received considerably attention in the literature, see for review.

Returning to the definition (28), we see that the complexity gives the leading order asymptotic behaviour of the mean number of fixed points assuming exponential growth with increasing $N$. We have seen above that the exponential growth assumption is indeed justified beyond the threshold $\hat{\sigma}_{c}=1$. In other words, we have

for $\hat{\sigma}>1$. Clearly, this is a rather crude approximation since the error term (albeit sub-exponential) is not forbidden to grow with $N$. Even worse, the complexity tells us nothing about the number of fixed points below the threshold ($\hat{\sigma}<0$) except that mean growth must be sub-exponential. In order to go beyond the complexity, we need to perform a more detailed analysis. The main property that we will use for this analysis is a relation between the mean value of the absolute determinant (21) and the mean spectral density of the real eigenvalues (22). This relation is so important that we will state it as a separate lemma.

Let $\lambda$ be a real constant. Using the joint density function (14) in the definition of the partial expectation (20), we get

The Vandermonde determinant (15) can be factorised as

Now, using the structure of the joint density function (14) once again, we see that

If we set $\lambda=1/\overline{\sigma}^{D}$, multiply both sides by $\overline{\sigma}^{ND}$, and sum over $n$, then we can recognise the right- and left-hand side as (21) and (22), respectively. This completes the proof. ∎

Lemma 3.3 is rather surprising: a priori, the left-hand side of the identity (37) represent a problem which depends on both real and complex eigenvalues of the product matrix, but it turns out that this can be reduced to a problem involving only real eigenvalues at the expense of increasing the matrix dimension by one. The lone matrix case ($D=1$) of Lemma 3.3 was shown in with a proof based on Householder reflection. The Householder reflection method has the benefit that we do not need to know the full structure joint density for the eigenvalues, but it also less suitable for generalisations.

Together with Proposition 3.2, Lemma 3.3 gives us an intuition for the origin of the phase transition at $\hat{\sigma}=1$ that found earlier. The right-hand side of equation (37) asks us to evaluate the mean spectral density for the real eigenvalues at $1/\overline{\sigma}^{D}$, but in the large-$N$ limit the mean spectral density concentrate on an interval with compact support. Dependent on our choice of $\sigma$, we may be in either a high or a low density region, which will result in a large or a small number of fixed points, respectively. The non-analyticity of the limiting mean spectral density at the edge of its support is the origin of the phase transition. We will study this more carefully below.

We already know that the appropriate scaling is set by $\hat{\sigma}=\overline{\sigma}/N^{1/2}$, thus using this scaling in (37) from Lemma 3.3, we get

The right-hand side in (45) is a so-called Meijer $G$-function, see e.g. [27, §16]. We note that

but refer to the literature for a more detailed description of Meijer $G$-functions and their properties.

For present purposes, we are interested in large-$N$ behaviour. An approximation for the ratio of normalisation constants (44) can be found using a Poincaré-type expansion for the Gamma functions [27, §5.11] which gives

An asymptotic expansion for the Meijer $G$-function is also known [31, §5.9]. To leading order, we have

By comparison with (46), we see that the leading term in the expansion (48) is exact for $D=1$.

Now, inserting the approximations (47) and (48) back into (43), we get

As alluded to earlier, the approximation of the spectral density depends on whether we are in the high density region ($\hat{\sigma}>1$) or in the low density region ($\hat{\sigma}<1$). In the high density region we have

We note that the leading $N$ behaviour is in agreement with our result for the complexity (35) obtained using a different method. It is also worth mentioning that while the logarithmic leading $N$ behaviour (i.e. the complexity) is independent of $\nu_{1},\ldots,\nu_{D-1}$, the sub-leading terms are not.

Evaluation of the mean spectral density in the low density region is trickier. Here, the real global spectral density defined by the limit (25) is zero. However, this does not imply that the finite-$N$ density is zero but rather that this region is dominated by rare events. We expect to have a ‘large deviation principle’ for the form

since the spectrum of the product matrix concentrating near the origin in this scenario. This implies that we have

for $\hat{\sigma}\ll 1$. In fact, it is expected that this approximation holds up to the threshold $\hat{\sigma}_{c}=1$, where the density develops a discontinuity. While this is difficult to prove, numerics (see Figure 2) leaves little doubt about its validity. We will verify it analytically for the single-layer ($D=1$) case only.

The finite-$N$ real spectral density is given by

with $\nu_{0}=0$. In the single layer case ($D=1$), the weight function is a Gaussian (46) and the sum in (55) can be expressed in terms of an incomplete gamma function. This allows integration over $x$, which yields

where $\Gamma(N,x)=\int_{x}^{\infty}dt\,e^{-t}t^{N-1}$ and $\gamma(N,x)=\int_{0}^{x}dt\,e^{-t}t^{N-1}$ are incomplete gamma functions. We are interested in $\lambda=N^{1/2}/\hat{\sigma}$ for which a saddle approximation gives

It is seen that the first term on the right-hand side is dominant for $\hat{\sigma}>1$, while the last term is dominant for $\hat{\sigma}<1$. Thus, we have

consistent with both (50) and (54). The behaviour near the critical value $\hat{\sigma}=1$ is given by

which is the so-called local edge regime .

Summary and outlook

In this paper, we studied the mean number of fixed points for a special class of multi-layered random dynamical systems. The class of systems we have studied may be considered as a null model for general multi-layered systems. Each layer was represented by a zero-mean Gaussian random map chosen (statistically) independent from the other layers, and was furthermore chosen to be homogeneous as well as domain- and codomain-isotropic. Our main result about such multi-layered systems was twofold.

First, we showed that asking for the mean number of fixed point of the aforementioned multi-layered random dynamical system is equivalent to an otherwise separate question within the framework of Random Matrix Theory (see Theorem 2.1). More precisely, to find the mean number of fixed points we can calculate the mean of the absolute value of the characteristic polynomial of a product of independent Gaussian matrices (also known as real Ginibre matrices ). This result was found using a general framework which have been build around the so-called Kac–Rice formula, see e.g. . This result is important for two main reasons: (i) it shows that the mean number of fixed points is a universal quantity in the sense specified in Section 2 and (ii) the random matrix problem is much easier to study both numerically and analytically.

Our second main result is an asymptotic expression for the mean number of fixed points in the high-dimensional limit. We found that the mean number of fixed point for our multi-layered system (2) with dimension $N$ and depth $D$ behave as

for $N\to\infty$. In other words, the large-$N$ limit of our multi-layered system (2) as defined in Section 2 has two phases: a phase with single fixed point for $\hat{\sigma}<1$ and a phase where the number of fixed points grows exponentially with $N$ for $\hat{\sigma}>1$. This type of transition from a ‘trivial’ landscape to ‘complex’ landscape has been observed in large number of models over the recent years, see e.g. . It has been suggested to refer to such transitions as topological trivialisation . It is intriguing that such topological trivialisation appears to be a relatively generic feature of high-dimensional non-linear systems.

We note that the single-layered systems ($D=1$) studied in this paper is closely related the (continuous-time) systems studied in . In fact, Fyodorov and Khoruzhenko looked for equilibrium points in a continuous-time system $d\mathbf{x}(t)/dt=-\mathbf{x}(t)+\mathbf{f}(\mathbf{x}(t))$ which is equivalent to looking for fixed points in the discrete-time system $\mathbf{x}(t+1)=\mathbf{f}(\mathbf{x}(t))$. Thus, the model considered in this paper is in direct correspondence to the Fyodorov-Khoruzhenko model, albeit the random maps in this paper is chosen in different way than in . It is therefore not surprising that the our result for the mean number of fixed points (60) with $D=1$ is identical to the result in . However, our model generalises the model by Fyodorov and Khoruzhenko by considering multi-layered systems ($D>1$). Furthermore, it reasonable to expect that it will be easier to study quantities such as periodic orbits within the discrete-time setting (as considered in this paper) compared to within the continuous-time setting (as considered in ). This is important since periodic orbits plays a crucial role in our understanding of the dynamical properties of complex systems.

We would like thank Yan Fyodorov for sharing a draft of . The work is part of a research program supported by the Australian Research Council (ARC) through the ARC Centre of Excellence for Mathematical and Statistical frontiers (ACEMS). PJF also acknowledge partial support from ARC grant DP170102028.

Appendix A Proof of main theorem

To proof Theorem 2.1 we will first show the single-layer case ($D=1$) and then the multi-layer case ($D\geq 2$). The reason for dividing the proof into two parts is that the structure of the single-layer case and the multi-layer cases differ slightly. On the other hand, the conceptual idea behind the proof is the same for both situations. Thus, by first understanding the single-layered case we hopefully make the generalisation to multi-layered cases more transparent.

where the expectation on the right-hand side is with respect to joint distribution of the random vector field $\mathbf{f}$ and the random Jacobian field $\nabla\mathbf{f}=(\partial f_{i}/\partial x_{j})_{ij}$ (the latter is an $N\times N$ matrix-valued field). In (61), we have used

to denote the Dirac delta of vector-valued argument and $\mathbf{I}_{N}$ to denote the $N\times N$ identity matrix.

We recall that the map $\mathbf{f}$ is assumed to have correlation function

where we have used the constraints (6). The second equality in (66) implies that the fields $\mathbf{f}$ and $\nabla\mathbf{f}$ are uncorrelated if evaluated at a common point $\mathbf{x}$, and since the fields are furthermore Gaussian this implies stochastic independence. Due to this independence, we may rewrite (61) as

The first expectation in (67) is trivial since the field $\mathbf{f}$ is a centred Gaussian, and we have

In order to evaluate the second expectation (67), we first note that due to homogeneity this expectation is in fact independent of the location $\mathbf{x}$. Thus, the latter expectation depends on a single Gaussian matrix-valued random variable $\mathbf{J}=\nabla\mathbf{f}(\mathbf{0})=(J_{ij})_{ij}$ with

where have used the notation $\sigma:=\sqrt{-\kappa^{\prime}(0)}>0$. We can recognise $\mathbf{J}=(J_{ij})_{ij}$ as a random matrix with i.i.d. centred Gaussian entries with variance $\sigma^{2}$, i.e. a matrix from the so-called real Ginibre ensemble. Thus, we may write the second expectation in (67) as

Finally, using the evaluations (68) and (70) in (67) and performing the integration over $\mathbf{x}$ proves the single-layer version of the Theorem 2.1.

Multi-layer case:

For $D\geq 2$, we must look for solutions to the system

with notation as in Section 2. In order to do so, we introduce

Thus, similar to the single-layer case, we can apply the Kac–Rice formalism, which gives an expression for the mean number of fixed points

We recall that the layer $\mathbf{f}^{(1)},\ldots,\mathbf{f}^{(D)}$ is chosen such that each layer is independent of the others and

Thus, similar to the single-layer case, we have

and consequently that the fields $\mathbf{F}$ and $\nabla\mathbf{F}$ are independent if evaluated at a common point $\mathbf{X}$. We can therefore write (75) as

Now, switching back to our original notation, we have

where each $\nabla\mathbf{f}^{(d)}=(\partial f_{i}/\partial x_{j})_{ij}$ denotes an $N_{d-1}\times N_{d}$ matrix-valued Gaussian field.

The first expectation in (79) is straightforward to evaluate since the layers are independent, and we have

In order to simplify the determinant which appear within the second expectation in (79), we will employ the following general determinant identity:

which holds for any matrices $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ of dimensions $n\times n$, $n\times m$, and $m\times n$, respectively. Successive use of this identity yields

Now, to evaluate the second expectation in (79), we note that the expectation is independent of the location $\mathbf{X}=(\mathbf{x}^{(0)},\ldots,\mathbf{x}^{(D-1)})$ due to homogeneity of the field $\mathbf{F}$. So, in complete analogue to the single-layer case, we can introduce matrices $\mathbf{J}_{1}=\nabla\mathbf{f}^{(1)}(\mathbf{0}),\ldots,\mathbf{J}_{D}=\nabla\mathbf{f}^{(D)}(\mathbf{0})$ and write

It follows from (77) that each matrix $\mathbf{J}_{d}=(J^{(d)}_{ij})_{ij}$ has covariance matrix

where have used the notation $\sigma_{d}:=\sqrt{-\kappa_{d}^{\prime}(0)}>0$. Thus, we recognise $\mathbf{J}_{d}$ as an $N_{d-1}\times N_{d}$ random matrix with i.i.d. centred Gaussian entries with variance $\sigma_{d}^{2}$. Moreover, since the fields $\mathbf{f}^{(1)},\ldots,\mathbf{f}^{(D)}$ are stochastically independent so are the matrices $\mathbf{J}_{1},\ldots,\mathbf{J}_{D}$. In other words, we have

The Theorem follows by inserting (A) and (83) into (79) and performing the integrals over $\mathbf{x}^{(0)},\ldots,\mathbf{x}^{(D-1)}$.

References