Critical scaling in hidden state inference for linear Langevin dynamics

Let us consider a generic network where only the dynamics of a subnetwork of nodes is observed while the others are hidden and form what we call the “bulk”. The indices $i,j=1,...,N^{\rm b}$ and the superscript b are used for the hidden or bulk variables; similarly the indices $a,b=1,...,N^{\rm s}$ and the superscript s for the observed or subnetwork nodes of the network. Assuming linear couplings $\{J_{ij}\}$, $\{K_{ia}\}$ between hidden and observed variables $x_{i}$, $x_{a}$, their dynamical evolution is described by

$J_{ij}$ (respectively $J_{ab}$) gives the hidden-to-hidden (respectively observed-to-observed) interactions, while the coupling between observed and hidden variables is contained in $K_{ia}$ and $K_{aj}$. $\lambda$ is a self-interaction term acting as a decay constant and providing the basic timescale of the dynamics. The dynamical noises $\xi_{i}$, $\xi_{a}$ are Gaussian white noises with zero mean and diagonal covariances $\Sigma_{i}$, $\Sigma_{a}$

The application of the extended Plefka expansion to this problem with hidden and observed nodes [5] allows us to obtain a closed system of integral equations describing the second moments of the Gaussian posterior distribution over the hidden dynamics, i.e. conditioned on the observed trajectories. These second moments are denoted by $C_{i}(t,t^{\prime})$, the hidden posterior variance, $R_{i}(t,t^{\prime})$, the hidden posterior response, $B_{i}(t,t^{\prime})$ and $B_{a}(t,t^{\prime})$, the posterior correlations of auxiliary variables introduced to represent the hidden and observed dynamics respectively à la Martin–Siggia–Rose–Janssen–De Dominicis [12, 13, 14]. They are diagonal in the site index but functions of two times by construction of the Extended Plefka Expansion [6], which gives an effectively non-interacting approximation of the dynamics where couplings among trajectories are replaced by a memory term (related to $R_{i}(t,t^{\prime})$) and a time-correlated noise (whose covariance is related to $B_{i}(t,t^{\prime})$).

Our main interest is in the posterior statistics: the equal-time variance of inference errors gives a mean-square prediction error, while the time correlations between the inference errors define a posterior relaxation time. To simplify the equations derived in [5], we consider long times, where a stationary regime is reached: all two-time functions become time translation invariant (TTI) and the Laplace-transformed equations yield a system of four coupled equations for $\tilde{C}_{i}(z)$, $\tilde{R}_{i}(z)$, $\tilde{B}_{i}(z)$, $\tilde{B}_{a}(z)$. Let us briefly recall it here, as it will be the starting point for the analysis

We expect the extended Plefka approach to give exact values for posterior means and variances in the case of mean-field type couplings (i.e. weak and long-ranged), in the thermodynamic limit of an infinitely large system. More precisely we define the thermodynamic limit as the one of an infinitely large bulk and subnetwork, $N^{\rm b},N^{\rm s}\rightarrow\infty$ at constant ratio $\alpha=N^{\rm s}/N^{\rm b}$. For the mean-field couplings we assume, as in [6], that $\{J_{ij}\}$ is a real matrix belonging to the Girko ensemble [15], i.e. its elements are independently and randomly distributed Gaussian variables with zero mean and variance satisfying

The parameter $\eta\in[-1,1]$ controls the degree to which the matrix $\{J_{ij}\}$ is symmetric. In particular, the dynamics is non-equilibrium – it does not satisfy detailed balance – whenever $\eta<1$. Similarly, $\{K_{ai}\}$ is taken as a random matrix with uncorrelated zero mean Gaussian entries of variance

We have introduced amplitude parameters $j$ for $\{J_{ij}\}$ (hidden-to-hidden) and $k$ for $\{K_{ai}\}$ (hidden-to-observed) here. The scaling of both types of interaction parameters with $1/\sqrt{N^{\rm b}}$ ensures that, when the size of the hidden part of the system increases, the typical contribution it makes to the time evolution of each hidden and observed variable stays of the same order. The high connectivity, where all nodes interact with all other ones, implies that in the thermodynamic limit all nodes become equivalent. Local response and correlation functions therefore become identical to their averages over nodes, defined as

As a consequence, all site indices can be dropped and the correlation and response functions can be replaced by their mean values, $\tilde{C}_{i}(z)\equiv\tilde{C}^{\rm bb|s}(z)$, $\tilde{R}_{i}(z)\equiv\tilde{R}^{\rm bb|s}(z)$, $\tilde{B}_{i}(z)\equiv\tilde{B}^{\rm bb|s}(z)$. The superscripts bb$|$s here emphasize that we are considering moments conditioned on subnetwork values though for brevity we drop them in the rest of the calculation. Similarly, for the correlations of auxiliary variables related to observations we can set $\tilde{B}_{a}(z)\equiv\tilde{B}^{\rm ss}(z)$. Of primary interest is then $\tilde{C}(z)$, the Laplace transformed posterior (co-)variance function of prediction errors, which as a site-average can be thought of as a macroscopic measure of prediction performance. This should become self-averaging in the thermodynamic limit. To see this, consider e.g. the sum $\sum_{j}J_{ij}J_{ji}\tilde{R}_{j}(z)$. Replacing $\tilde{R}_{j}$ by $\tilde{R}$, the prefactor is a sum of $N^{\rm b}$ terms so converges to its average $N^{\rm b}\eta j^{2}/N^{\rm b}=\eta j^{2}$ in the thermodynamic limit. So, as in [6], we can replace

Making this and similar substitutions in the system (2.3), and choosing scalar noise covariances $\Sigma_{i}=\sigma_{\rm b}^{2}$ and $\Sigma_{a}=\sigma_{\rm s}^{2}$ as in [9], one finds

Here $\tilde{B}^{\rm ss}(z)=-1/(\sigma_{\rm s}^{2}+k^{2}\tilde{C}(z))$ has already been substituted into equations (2.10a) and (2.10c).

We comment briefly on the relation of the above results to the Fluctuation Dissipation Theorem (FDT), which in terms of Laplace transforms reads

This can be compared with the expression for $\tilde{R}(z)$ that follows from (2.10b) (taken for $z$ and $-z$)

The r.h.s. of the FDT is then

Comparing with (2.11), the FDT is satisfied for symmetric couplings ($\eta=1$) as expected, while there are progressively stronger deviations from FDT as $\eta$ decreases towards $-1$.

As a consistency check and to support our claim of exactness in the thermodynamic limit, we briefly compare our results with expressions for $\tilde{C}(z)$ and $\tilde{R}(z)$ that can be worked out by alternative means.

In general, from (2.10b) we can get an expression for $\tilde{R}(z)$ in terms of $\tilde{C}(z)$: this is (2.12). Substituting into (2.10a), $\tilde{B}(z)$ can also be worked out as a function of $\tilde{C}(z)$. Using these expressions for $\tilde{R}(z)$ and $\tilde{B}(z)$ in equation (2.10c), one finds a closed algebraic equation for the posterior variance $\tilde{C}(z)$

This is the same expression as obtained by calculations in [9] using an explicit average over the quenched disorder variables $J_{ij}$ and $K_{ia}$. Particular cases are also further validated by random matrix theory, as follows.

We would like first to understand how $\tilde{C}(\text{i}\omega)$ depends on the parameters $\lambda$, $j$, $k$, $\sigma_{\rm s}$, $\sigma_{\rm b}$, $\eta$ and $\alpha$. The last two of these are already dimensionless. By extracting the appropriate dimensional scales, we can reduce the other five parameters to only two dimensionless combinations.

From (2.1) one sees that $j$, $k$, $\lambda$ have dimensions $t^{-1}$, while the dimension of $\sigma_{\rm s/b}^{2}$ is $x^{2}t^{-1}$. We can build from these the dimensionless parameters $\gamma=j/\lambda$ and $p=\lambda/\sigma$. Here $\sigma=\sigma_{\rm b}k/\sigma_{\rm s}$, which has dimension $t^{-1}$ and contains the observation “intensity” $k$ as well as the ratio between the dynamical noises $\sigma_{\rm b}/\sigma_{\rm s}$. The latter is a third dimensionless parameter but as it only enters one prefactor we will not need to keep it separately.

Extracting appropriate dimensional amplitudes for all four two-point functions, we write them as

Here $\Omega=\omega/\sigma$ is a dimensionless frequency; similarly $\mathcal{C}_{\alpha,\gamma,\eta,p}(\Omega)$, $\mathcal{R}_{\alpha,\gamma,\eta,p}(\Omega)$, $\mathcal{B}_{\alpha,\gamma,\eta,p}(\Omega)$ and $\mathcal{B}^{\rm ss}_{\alpha,\gamma,\eta,p}(\Omega)$ are dimensionless and depend on the dimensionless parameters $\alpha$, $\gamma$, $\eta$ and $p$: for the sake of brevity, we do not write the subscripts indicating this dependence in the following. Let us briefly comment on (3.1a). One sees that $\tilde{C}(z)$ is directly proportional to $\sigma_{\rm s}^{2}$ and inversely proportional to $k^{2}$: the weaker the hidden-to-observed coupling and the stronger the dynamical noise acting on the observed variables, the less information one can extract from the subnetwork trajectories and the more uncertain the predictions for the behaviour of the bulk.

To summarize, we switch from eight original parameters $\{\alpha,\eta,\lambda,j,k,\sigma_{\rm s},\sigma_{\rm b},\omega\}$ to a set of five dimensionless parameters $\{\alpha,\eta,\gamma,p,\Omega\}$. For the dimensionless second moments (3.1), the system (2.10) becomes

where we have dropped the frequency argument and introduced the shorthand $\mathcal{R}_{-}=\mathcal{R}(-\Omega)$. The solution for $\mathcal{B}^{\rm ss}$, already taken into account by substitution in (3.2), is given by $-1/(1+\mathcal{C})$.

We next analyze in the parameter space $\alpha,\gamma,p$ (at fixed $\eta$) the singularities of $\mathcal{C}(0)$, the (dimensionless) zero frequency posterior covariance. $\mathcal{C}(0)$ is a convenient quantity to look at in the first instance as it helps one detect where qualitative changes of behaviour occur. We will see that these then show up also in the relaxation time and the inference error itself, i.e. the equal time posterior variance $C(t-t)=C(0)$. Their behaviour will be addressed respectively in sections 3.4, 3.5 and 3.8, 3.9, while in sections 3.6 and 3.7 we clarify further the interpretation of our results for the relaxation times by looking at the shape of the correlation functions.

Independently of $\eta$, we find two critical regions that are shown graphically in figure 1:

1. 1.

   $\forall p$, $\alpha=0$ and $\gamma>\gamma_{c}$

2. 2.

   $\forall\gamma$, $p=0$ and $0<\alpha<1$

The first case gives back the dynamics without observations ($\alpha=0$), for which $\gamma<\gamma_{c}=1/(1+\eta)$ is the condition of stability beyond which trajectories typically diverge in time (see [6]). Interestingly, as soon as $\alpha>0$, the constraints from observations make the solution stable irrespective of whether $\gamma$ is smaller or bigger than the critical value. For $\gamma>\gamma_{c}$ the observed trajectories would then be divergent, and so would the predicted hidden trajectories, while the error (posterior variance) of the predictions would remain bounded. It is difficult to conceive of situations where divergent mean trajectories would make sense, however, so we only consider the range $\gamma\leq\gamma_{c}$ in our analysis. In what follows, we will plot all the curves that refer to this region of the parameter space in blue.

The second limit, $p\to 0$, corresponds, for fixed ratio between noises $\sigma_{\rm s}/\sigma_{\rm b}$, to $k\gg\lambda$: we call this scenario an “underconstrained” hidden system. In general for large $k$ the posterior variance decreases as $1/k^{2}$, as used in the scaling (3.1a). But for $\alpha<1$ there are directions in the space of hidden trajectories that are not constrained at all by subnetwork observations, and their variance will scale as $1/\lambda^{2}$ instead. These directions give a large contributions to the dimensionless $\mathcal{C}(\Omega)$ that diverges in the limit $k\gg\lambda$. In general one has a similar effect when $k/\sigma_{\rm s}\gg\lambda/\sigma_{\rm b}$, where the noise in the dynamics acts to effectively reduce the relevant interaction or decay constant. This behaviour is broadly analogous to what happens in learning of linear functions from static data [10, 11]: there the prediction error will also diverge when $\alpha$, which is then defined as the ratio of number of examples to number of spatial dimensions, is less than unity and no regularization is applied. For the curves belonging to this second region we choose the colour red.

Close to the two critical regions in the space of $\alpha$, $\gamma$ and $p$, $\mathcal{C}(0)$ is expected to exhibit a power-law dependence on the parameters specifying the distance away from the critical point. One can study this behaviour by using standard scaling techniques developed for the study of critical phase transitions [17]. The aim of this analysis, which we describe in A and summarize in the next section, is to find the master curves and associated exponents that describe the approach to the critical point(s). From the master curves for the power spectrum we can then derive (B) the corresponding scaling behaviour of the relaxation times and the inference error, which is discussed in sections 3.4–3.9 below.

The master curves of the power spectrum describe its shape in the critical regions and for small frequencies $\Omega$ of the order of a parameter-dependent characteristic frequency $\Omega^{*}$ (see A). These master curves are useful because they highlight the two main trends (with $\eta$ and with $\alpha$) that will manifest themselves in the relaxation time and in the inference error. In the first critical region where there are few observations, the conditional dynamics is dominated by hidden-hidden interactions. Their degree of symmetry $\eta$ therefore plays a key role and in fact determines a non-trivial crossover in the critical behaviour (see figure 9). As developed in A.1, this crossover from equilibrium to generic non-equilibrium dynamics can be studied by including $\epsilon=1-\eta$ as a small parameter in the scaling analysis. By contrast, in the second critical region the scaling functions do not depend on $\eta$ (see figure 8, right). The intuition here is that the critical behaviour is dominated by whether a direction in the hidden trajectory space is constrained by observations or not, i.e. by hidden-to-observed rather than hidden-to-hidden interactions.

Regarding the trend with $\alpha$, the main feature is the decrease in $\mathcal{C}(0)$ with increasing $\alpha$ meaning that with many observations the predictions for the hidden dynamics will become arbitrarily precise, in line with intuition. This can be seen in figure 2, where one observes also that at $\Omega\sim O(1)$ all spectra collapse into a Lorentzian tail. This indicates an exponential decay of the correlations between prediction errors in the temporal domain. As the amplitude of the tail is largely $\alpha$-independent, the typical time of this exponential decay decreases with $\alpha$: with many observations, errors in the prediction of the hidden states become progressively less correlated with each other.

This already shows that from the power spectrum $\mathcal{C}(\Omega)$ one can extract useful information on relaxation times and the same is true for the inference error, as we analyze more systematically in the next sections. Importantly, we find nontrivial power-law dependences on $\alpha$ that arise from the dynamical nature of our problem and have no analogue in static learning  [10, 11].

We look at the relaxation time, which is a measure of time correlations in the errors of inferred hidden values. We study in particular how it depends on the number of observations and the interaction parameters.

The relaxation time can be defined in a mean-squared sense as

Given this relation, results for power spectra can be directly used to obtain the master curves for relaxation times, see B.1 and B.2. To construct these master curves, let us introduce the dimensionless version of the relaxation timescale,

We look separately at the two critical regions, starting in this subsection with the first one, i.e. in the vicinity of the “no observations” critical point. We summarize here the main findings, leaving details for B.1. The distance from the singularity is controlled by $\alpha$ itself and $\delta\gamma=\gamma_{c}-\gamma$. As mentioned in the overview in section 3.3, the behaviour in this region depends crucially on the symmetry parameter $\eta$ and so we break down the results accordingly.

We next discuss the behaviour of the relaxation time in the second critical region, $\alpha\to 1$ (i.e. $N^{\rm s}=N^{\rm b}$) and $p\to 0$ (i.e. $k\gg\lambda$ at fixed $\sigma_{\rm s}/\sigma_{\rm b}$). The distance from the critical point is therefore represented by $p$ itself and by $\delta\alpha=\alpha-1$. As is shown in B.2, when both are small the results are independent of the degree of symmetry $\eta$ of the interactions among the hidden variables. For the dimensionless relaxation time we find for positive $\delta\alpha$ (see B.2)

and this is consistent with numerical data for $\alpha>1$, see figure 3 (right). The figure also shows that as $\alpha$ decreases below unity the relaxation time reaches a plateau, whose height can be worked out as $1/p\sqrt{1+\gamma^{2}}$. The plateau and the decay in eq. (3.7) are obtained as the two limiting behaviours of an $\eta$-independent master curve that applies when $\delta\alpha$ is of order $p$.

To understand the behaviour across the entire range $0<\alpha<1$ one can work out a separate master curve for fixed $\alpha$ and $p\to 0$ (see Eq. (B.12) in B), which is plotted in figure 4 (left) for $\eta=-1$ along with numerical results. This curve does have a nontrivial dependence both on $\eta$ and on $\alpha$; for $\alpha\to 0$ it approaches a constant, whose value can be checked to be consistent with the $\alpha=0$ result (Eq. (A.1)) as it must.

A comparison of figures 3 (right) and 4 (left) shows that the range around $\alpha=1$ where relaxation times are independent of $\eta$ depends not just on $p$ but also on $\gamma$. For small $\gamma$ as in figure 3 (right) we see good agreement with the $\eta$-independent master curves down to at least $\alpha\approx 0.5$. For $\gamma=1$, the value we take to plot curves for $\eta<0$, figure 4 (left) shows that already at $\alpha\sim 0.9$ there is visible variation across different values of $\eta$.

Figure 4 (left) contains a further insight: at fixed $\gamma$ and moderate $\alpha$, varying $\eta$ can shift the system away from the second critical region ($p\to 0$) towards the first one ($\delta\gamma\to 0$). Indeed, if $\gamma=1$, $\delta\gamma=\gamma_{c}-\gamma$ becomes small when $\eta$ is close to zero because $\gamma_{c}=1/(1+\eta)$. Accordingly for the $\eta=-0.0001$ curve in the figure, the relaxation time for small $\alpha$ is captured better by the master curve for the first critical region (dark blue line with crosses) than the one for the second region. At fixed $\eta$ and for appropriate combinations of small $\delta\gamma$ and $p$ one can even see both critical regions as $\alpha$ is increased. Figure 4 (right) shows a case in point.

As a general trend in the relaxation times one sees that they decrease significantly with increasing $\alpha$ (see e.g. figure 4, right): as the values of the hidden variables become constrained increasingly strongly by those of the observed ones, the remaining uncertainty in the prediction becomes local in time. There are surprising counter-examples, however. Figure 4 (left) shows that for strongly anti-symmetric interactions, e.g. $\eta=-1$, the relaxation time $\tau$ increases with $\alpha$ up until $\alpha=1$. This unusual trend is linked to a change of shape in the decay of the time-dependent correlation function, which is strongly non-exponential for $\alpha<1$ but becomes close to exponential beyond (see the inset of figure 4, left).

Having studied overall relaxation timescales, we next look more comprehensively at the time correlation functions $C(t-t^{\prime})$, which can be obtained numerically by inverse Fourier transform of the power spectrum.

For the first critical region we refer to figure 5 (left). Here the master curves are the result at $\alpha=0$ from [6]: for $\eta=0$ and $\eta=1$ these are given respectively by a pure exponential and by an exponentially weighted integral of a modified Bessel function $I_{1}$. In the latter case the correlation function has a $1/\sqrt{t-t^{\prime}}$ regime indicated in the figure, before crossing over to $(t-t^{\prime})^{-3/2}$ times an exponential cutoff (which kicks in around $t-t^{\prime}=250$ in the case shown in the figure). For intermediate symmetries (e.g. $\eta=0.85$ in figure 5 left) one has a crossover between this behaviour and the pure exponential for $\eta=0$. This crossover occurs at a time $\sim 1/\epsilon^{2}$ (here $\epsilon=1-\eta$ as before), which is consistent with the frequency space crossover at frequency $\sim\epsilon^{2}$, see e.g. figure 9 (right).

For $\eta=-1$, in the $\alpha=0$ case $C(t-t^{\prime})$ is a Bessel function of the first kind $J_{1}$ [6]. This has oscillations, and we find that this behaviour persists as $\alpha$ is increased. This is shown in the inset of figure 5 (right) where the oscillations exhibit a longer decay time for $\alpha=0.5$ consistently with the findings on relaxation times of figure 4 (left). Note that small oscillations are also present in the curve $\alpha=0.5$ of figure 4 (left, inset) on longer times than shown in the graph.

We show an example of $C(t-t^{\prime})$ in this second critical region in figure 5 (right). The master curve for the power spectrum in this region (see A.2) for $\delta\alpha\to 0$ tends to $1/\sqrt{1+(\Omega/p\sqrt{\gamma^{2}+1})^{2}}$, independently of $\eta$. Its inverse Fourier transform gives a modified Bessel function $K_{0}(\Omega/p\sqrt{\gamma^{2}+1})$ in the time domain, which is plotted in the figure alongside the numerical results (dashed red line with circles).

We turn finally to the behaviour of the inference error for the prediction of hidden unit trajectories. This is given by the equal time posterior correlator

where $\mathcal{C}_{0}$ is a dimensionless equal time posterior variance. We see that the size of the error is generically proportional to the noise acting on the dynamics of hidden and observed variables, and inversely proportional to the hidden-to-observed interaction strength. We summarize in the following subsections our results for $C(0)$, with a focus on the $\alpha$-dependence, and leave details to B.3 and B.4. We begin with the first critical region.

In the second critical region, the dominant terms of the integral (3.8) can be shown to scale near $\alpha=1$ as $1-\alpha$ for $\alpha<1$ and $\mbox{const.}-\ln(\alpha-1)$ for $\alpha>1$ (see B.4). These terms are plotted (red line) in figure 7 where the linear scale inset clearly shows the linear dependence on $1-\alpha$. A separate $p\to 0$ master curve for fixed $\alpha$ below 1 can also be derived (by integrating the solution of (A.36)) and matches smoothly to the $1-\alpha$ behaviour around $\alpha=1$.

As a common trend across the two critical regions we have the intuitively reasonable result that the inference error decreases when the number of observed variables gets bigger. Figure 6 (right) summarizes this and shows that, as in the case of the relaxation times, the behaviour with increasing $\alpha$ can cross over from the first critical region to the second provided that both $\delta\gamma$ and $p$ are small. It is worth re-emphasizing the non-trivial power law dependences of the inference error on $\alpha$ that occur in our dynamical setting and are quite distinct from the simpler behaviour in static learning scenarios [10, 11].