Gaussian approximation of dynamic cavity equations for linearly-coupled stochastic dynamics

Major comments

We thank the reviewer for suggesting additional references. We have carefully considered them and decided to include them in the revised manuscript.

Minor comments

1) In the derivation, Sec. 2.1, I would be more explicit on the precise reason why linear interactions allow for graphical model construction presented in the paper, maybe contrasting this with the more complicated situation that would occur for non-linear interactions (this is briefly mentioned in the conclusions, bit I feel that this point should be stressed more in the first part of the paper).

In the case of linear couplings, the disentanglement of the locally loopy structure of the factor graph into a locally tree-like structure is straightforward. For nonlinear interactions—such as terms of the form $h\left(\sum_{j=1}^N a_{ij} J_{ij} x_j(t)\right)$, with $h(x)$ a generic nonlinear function—this decoupling remains, in principle, possible but requires the introduction of auxiliary fields. This significantly complicates the factor graph representation and the derivation of the dynamical cavity equations.

2) In the abstract and introduction, the authors state that the formalism developed in this work can be applied to systems subject to global constraints, and in Sec. 3, they discuss the example of a spherical constraint. Could the authors comment on whether they could treat the case in which constraints are imposed uniformly on the continuous variables (for example, their positivity, as needed in applications to ecosystems dynamics), which are perhaps not straightforward to enforce using Lagrange multipliers?

In the manuscript, we focused on constraints that are either global but analytically tractable—such as the spherical constraint, which can be enforced via a Lagrange multiplier—or local constraints acting on individual nodes, as in observation-constrained dynamics. Uniform constraints such as positivity are more challenging to enforce within this framework. However, for models relevant to ecosystem dynamics, we believe that explicit enforcement of positivity is not strictly necessary. The presence of an absorbing state at zero ensures that variables initialized in the positive regime either remain positive or are absorbed at zero, which effectively constrains their dynamics without requiring hard constraints.

3) In the conclusion, the authors claim that while extending their equations to include nonlinearities in interactions—such as those found in neural networks and learning models—is challenging and intricate, this is not the case for multi-body nonlinearities, such as those in p-spin models with p≥3. Could they provide some intuition for why this is the case and how general is this statement?

In the case of multi-body nonlinearities, such as those arising in $p$-spin models with $p \geq 3$, the structure of the dynamical generating functional remains amenable to a factor graph representation. Specifically, one can write a factorized form of the dynamical partition function similar to Eq. (5), where the interaction term involves $p$ variables—for instance, $\hat{x}_i$, $x_{j_1},\dots,x_{j_{p-1}}$. This naturally leads to a bipartite graphical model in which factor nodes encode the $p$-body interactions. The resulting graph remains locally tree-like in sparse topologies, and a procedure analogous to that used in our work—based on the cavity method and a small-coupling expansion—can be applied. One thus obtains a closed set of equations for the local means, correlations, and responses.

Comments

1) In section 3.4, the authors assume that the cavity quantities do not depend on the edjes. This is not trivial at all due to the disorder, specially for small $K$. The authors should clarify on this.

This remark is certainly valid for a single instance, because of the random realization of couplings, and it could be relevant for ensemble averages. In order to obtain ensemble average results for sparse random regular graphs, that could be compared with classical results for fully-connected ones, we assumed that the disordered-averaged cavity response functions and cavity correlation functions are identical on every edge. This is reasonable because the 2-spin model does not show replica symmetry breaking—it is a disguised ferromagnet—and the aging behavior has a purely dynamical nature. In the manuscript we now state clearly this is an hypothesis. We checked through Monte Carlo simulations that the disorder-averaged correlation functions of different nodes are comparable, which suggests that our hypothesis is correct in the thermodynamic limit.

2) I suggest the authors to introduce a short discussion in the introduction making explicit the difference between their approach and the ones in [44] and [48].

In short, our method is developed as a small-coupling expansion at the single-disorder-instance level. This distinguishes it from [44], which performs a large-connectivity (i.e., small-coupling) expansion on the disorder-averaged problem. Our approach retains sample-specific information and is thus more general, allowing it to be employed as an algorithm capable of capturing finite-size effects and node-level fluctuations. In addition, while the DMFT framework proposed in [48] is designed for directed sparse graphs with nonlinear interactions, our method can deal with both directed and undirected graphs with linear interactions. Furthermore, our formulation provides a closed set of equations for the average quantities (means, response functions, and correlations), whereas [48] relies on a population dynamics algorithm to compute disorder-averaged observables.

Comment on heterogeneous networks

We thank the referee for this important suggestion. In the revised manuscript, we extend our analysis to heterogeneous interaction networks—Erdős-Rényi (ER) and Random Configuration Models (RCMs) graphs. We consider ferromagnetic couplings in all cases and also include bimodal disordered interactions specifically for the ER case.

This extension is significant, as only in the regular random graph with ferromagnetic couplings the interaction matrix has an analytical closed-form spectral density. In contrast, the spectral properties of the heterogeneous cases are not analytically tractable. We show that the Gaussian Expansion Cavity Method (GECaM) accurately reproduces equilibrium correlations and captures finite-size effects, confirming its applicability to a broad class of networks.

Major comments

1) I might have misunderstood figure 1, but it seems that the interactions between variables at different times do not correspond to the graphical structure shown. I don't see why the graph in figure 1 includes an interaction between $\hat{x}^n_i$ and $x^{n−1}i$ at the previous time step, since this term is absent from eq. (5). The authors also mention on page 5 that the local tree-like structure of the graph is a consequence of the linear coupling between variables. I think it'd be important to elaborate further on this point and better explain why linear interactions allow one to disentangle the loopy structure.

We thank the referee for pointing this out. Regarding the figure, the referee is correct; there was an inconsistency in the graphical structure, which we have now corrected in the revised version. As for the second part, we agree that our original wording may have been unclear. What we intended to convey is that, in the case of linear couplings, the disentanglment of the locally-loopy structure of the factor graph along a locally tree-like structure is straightforward. For nonlinear interactions—such as when the dynamics involves a term like $h\left(\sum_{j=1}^N a_{ij} J_{ij} x_j(t)\right)$ with $h(x)$ a generic nonlinear function—this decoupling is still in principle possible but requires the introduction of auxiliary fields. This significantly complicates the factor graph representation and the derivation of the dynamical cavity equations.

2) Regarding eq. (17), I have the following issue. If one expands eq. (17) up to second order in $\alpha$, one should recover the expansion of eq. (6), but this does not seem to be the case. The reason is that the $O(\alpha^2)$ contribution coming from the term with $\mu_{k\setminus i}^n$ in eq. (17) does not appear to cancel out with the $O(\alpha^2)$ contribution in eq. (6). Please clarify this point.

The key aspect here is that the correlation function $C_{i\setminus j}^{n,n'}$ in eq. (17) is defined as connected, i.e.,

3) On page 8 the authors mention that $B_{i\setminus j}(t,t')$ is zero because the dynamics is causal. In the generating functional formalism, the moments of conjugate variables are zero because the generating functional is normalized. It would therefore be important to elaborate further on why causality implies that these quantities vanish in the present context.

We agree that stating causality alone implies vanishing conjugate moments may be misleading. In the exact generating functional formalism, these moments vanish due to normalization. However, since cavity messages are derived from an approximation, we cannot rigorously prove their normalization. Instead, we show self-consistently that conjugate moments vanish under the assumption of causal dynamics. This argument would not hold for non-causal dynamics—e.g., observation-constrained dynamics as in Ref. [27]—where constraints introduce feedback from future times and the conjugate moments are indeed non-zero.

4) In the long time limit, eq. (50) has a stable solution if all eigenvalues of $J−\lambda$ have negative real parts. Besides that, the spectral density of, for instance, a sparse and symmetric random matrix $J−\lambda$ typically has unbounded support for $N \to \infty$, which makes the linear dynamics described by eq. (50) unstable. It is thus important to make explicit the conditions on $J−\lambda$ that render the linear dynamics stable

We now clarify in the manuscript that stability requires all eigenvalues of $J - \lambda$ to have negative real parts, and that the analysis is restricted to parameter regimes where this condition holds.

5) The authors mention on page 14 that the relation between random matrix theory and the Gaussian solution for the cavity equations is double-sided. Equations (45) and (48) for the response function are exactly the same as the cavity equations for the diagonal elements of the resolvent of sparse symmetric matrices (in this sense, it'd be instructive to cite a few relevant papers where these equations have been originally derived, e.g., works of Perez-Castillo, Rogers, Metz, Biroli, etc). On the other hand, the spectral density of $J−\lambda$ determines the correlation and response of the linear dynamical system, but this is only valid for graphical models whose spectral density has a bounded support (see point 4). Is it not correct to think that this restriction should also manifest itself in the Gaussian solution of the dynamical cavity equations?

Yes, this is correct. If the system is unstable, the TTI regime is never reached, and the Laplace-transformed GECaM equations are not defined. We clarify this point in the revised manuscript.

6) In section 3.1, the authors solve the equations for the Laplace transform of the response and the correlation functions in the TTI regime, in the case of random regular graphs. Technically speaking, this problem is equivalent to solving the resolvent equations for RRG's, so this section does not present new material. It would be interesting to explore further the consequences of the solution. For instance, how do the behaviours of the two-point functions compare with those in the fully-connected case?

We used our method to compute the equilibrium correlation functions for random regular graphs with ferromagnetic couplings and compared them with the fully-connected (FC) case. We observed that sparseness leads to a slower relaxation of the system. In particular, for small enough connectivity values, the decay of the correlation function deviates from a purely exponential form, exhibiting a fast initial decay followed by a slower relaxation tail. This illustrates how network topology—specifically, the finite connectivity—can qualitatively affect the dynamical behavior even in absence of disorder.

7) On page 24, it is not clear why one can assume that the local quantities are independent of the site index. Although the graph is regular, the coupling strengths are random variables, and local quantities should fluctuate from site to site. Please, clarify the status of the homogeneity assumption in this case.

This remark is certainly valid for a single instance, because of the random realization of couplings, and it could be relevant for ensemble averages. In order to obtain ensemble average results for sparse random regular graphs, that could be compared with classical results for fully-connected ones, we assumed that the disordered-averaged cavity response functions and cavity correlation functions are identical on every edge. This is reasonable because the 2-spin model does not show replica symmetry breaking—it is a disguised ferromagnet—and the aging behavior has a purely dynamical nature. In the manuscript we now state clearly this is an hypothesis. We checked through Monte Carlo simulations that the disorder-averaged correlation functions of different nodes are comparable, which suggests that our hypothesis is correct in the thermodynamic limit.

8) Equations (108-111) are equivalent to those for a ferromagnetic model with $J_{ij}=J$. The solution of eqs. (108-111) exhibits aging, which shows that in the sparse regime there is no need for randomness in the coupling strengths to observe aging effects. I believe this is an interesting conclusion from eqs. (108-111). It would be also important to clarify the meaning of the fully-connected limit in this context, because the model is ferromagnetic for finite $K$, while we recover the Cugliandolo-Kurchan equations for $p=2$ spherical model with Gaussian couplings in the limit $K\to\infty$.

The equations for a ferromagnet would be indeed the same as Eqs. (108-11). The catch is that while for the bimodal disordered case the coupling $J$ needs to scale as $1/\sqrt{K}$, for the ferromagnetic case it needs to scale as $1/K$. Hence, it is correct that the ferromagnetic model also shows the same sort of aging in the sparse regime (in fact, also the disordered case is just a disguised ferromagnet). In the large connectivity limit, instead, the ferromagnetic model does not show any aging because the terms of order $J^2$ vanish (because of the scaling $J\sim 1/K$), while the disordered model still preserves aging behavior as known from literature. We also verified numerically that this behavior is observed for the ferromagnetic model.

Minor comments

a) After eq. (7), what do the authors mean by the term "quasi-probability distribution"?

By quasi-probability distribution, we refer to a mathematical object similar to a probability distribution that relaxes some of Kolmogorov’s axioms. In our case, normalization holds, but non-negativity is not guaranteed—there may be regions where the function takes negative values.

b)–h)

We thank the referee for the suggestions and have accordingly corrected the manuscript.