1 Introduction↩︎

Stochastic dynamical systems are ubiquitous as models for natural data. They are perfectly suited for application to time-series data, and therefore also a good candidate to characterize systems that reach a steady state in the limit. If a system is governed by some stochastic differential equation (SDE) and the same system is observed under different interventions, ideally one would learn the underlying parameters governing the dynamics, and guarantee accurate prediction under new interventions.

However, in many natural settings, data is modeled as following an SDE even if one does not have access to explicit trajectories. Studies of ecological systems focus on the long-term survival of multiple species modeled by the quasi-stationary state of SDEs with environmental factors as perturbations [1]. The application of flow cytometry to protein signaling networks under perturbation [2] is destructive and yields protein quantification at one time point, modeled using the stationary distributions of linear SDEs in [3].

One highly motivating application is single-cell genomic sequencing with high-throughput CRISPR perturbations. Biologists are often interested in inferring the gene regulatory network (GRN) that characterizes the dynamics of gene expression, informing which genes should be targeted for treatment [4]. But the destructive nature of sequencing makes it impossible to observe the trajectory of a single cell at multiple time-points, and in general it is difficult to obtain any time-series genomics data due to the high expense. Therefore, practitioners often only collect data at the end of an experiment, i.e., from the stationary distribution of the system.

Understanding the dynamics is essential for extrapolating to unseen settings, but noise and latent confounding makes it non-trivial to determine the true dynamics. Causal disentanglement aims to learn causal factors in spite of these confounders, mainly focusing on directed acyclic graph (DAG) based methods. To demonstrate these methods are well-founded, there is considerable effort devoted to understanding which models have identifiability guarantees [5]. However, these models suffer from inherent weakness, in particular 1) being unable to represent cycles or 2) approximate continuous-time dynamical models.

There has been renewed interest in modeling with stochastic differential equations (SDE) directly [6]. In the genomic context, there is precedent for this type of modeling to represent the so- called “Waddington landscape” [7], the hypothetical energy surface of cells. Furthermore, SDEs are commonly used for simulating transcriptomic datasets from a given gene regulatory network [8], [9].

As demonstrated in the context of diffusion models, SDEs are fully expressive in practice and can accurately generate observational data [10]. But to identify the true underlying SDE requires assumptions on the model. Foundational theoretical works on identifying dynamical systems typically learn from many trajectories or even the infinitesimal generator of the dynamics [11], leaving open the harder setting of observing only the stationary distribution.

Our interest in this work is to verify which parametric assumptions are necessary for dynamical systems to have identifiability guarantees without trajectory data. Namely:

1.0.0.1 Contributions

In this work, we offer the first analysis of identifiability of interventional stochastic differential equations, with data restricted to the stationary measure. Specifically:

• We characterize tight bounds on the number of interventions necessary for identifiability of linear SDEs with shift interventions.

• We extend this analysis to nonlinear SDEs in the small noise regime, showing that identifiability is possible even without knowing the activation function of the true model.

• We apply this insight to synthetic data and semi-synthetic genomic data to confirm the efficacy of learned activations in causal SDEs, which improve expressiveness without sacrificing a simple structure and enable the inference of gene regulatory networks.

2 Setup↩︎

2.1 Notation↩︎

We will write the elementary basis vectors as \(\{e_i\}_i\). We will consider \(\sigma: \mathbb{R}^n \rightarrow \mathbb{R}^n\) as any elementwise function, i.e. \(\sigma_i(x) = \sigma_i(x_i)\). This includes elementwise activations as they are applied in multilayer perceptrons (MLPs), but we also allow for elementwise functions where each component acts differently. Writing \(\sigma'\) will, unless otherwise described, denote the map \(\mathbb{R}^n \rightarrow \mathbb{R}^n\) that applies elementwise the derivative of each component function, i.e. \(\sigma'(x) = [\sigma_1'(x_1), \dots, \sigma_n'(x_n)]\). We will use \(Jf\) to denote the Jacobian of a vector valued function \(f\), and \(\Delta f\) to denote the Laplacian of \(f\). We will use the notation \(diag\) to map vectors to diagonal matrices, or to map matrices to exclusively their diagonal part.

We use \(s_i(A)\) to denote the \(i\)th singular value of matrix \(A\). We let \(P_A\) denote the orthogonal projection onto the image of matrix \(A\), and \(P_A^\perp = I - P_A\) the orthogonal projection onto its complement. We let \(A^\dagger\) denote the pseudoinverse of \(A\), and note that if \(A\) has linearly independent columns it is a left inverse such that \(A^\dagger A = I\), likewise for rows and right inverse. We let \(\|\cdot\|\) denote the spectral norm and \(\|\cdot\|_F\) the Frobenius norm. We also introduce a minimal signed permutation distance \(d_P(X, Y) = \min_{\Lambda, \Pi} \|X - Y\Pi \Lambda\|\) where \(\Lambda\) is minimized over diagonal matrices with entries in \(\pm 1\) and \(\Pi\) is minimized over permutation matrices.

We will use \(\lesssim\) to denote inequality up to constant factor (treating some problem parameters as constants where specified), and similarly reserve the notation \(\tilde{C}\) for any absolute constant greater than zero. Finally iid indicates independent and identically distributed.

2.2 Stochastic Differential Equations↩︎

We consider SDEs of the following form, where \(v\) is a vector field and \(B_t\) is standard Brownian motion: \[\begin{align} dX_t = v(X_t)dt + \sqrt{\epsilon}dB_t. \end{align}\] We will only consider autonomous systems, i.e. where the drift and noise terms have no dependence on time \(t\). We will enforce the weak conditions on the drift and noise to guarantee a unique stationary distribution [12], with a density \(p\) that satisfies the Fokker-Planck equation, \[\begin{align} 0 = -\nabla \cdot (pv) + \frac{\epsilon}{2} \Delta p. \end{align}\]

2.3 Linear SDEs↩︎

We need some classical facts about linear SDEs, which are better understood than their nonlinear cousins, since Fokker-Planck can be solved explicitly.

2.4 Interventional SDEs↩︎

We focus on the setting where we only observe the SDE through the induced stationary density under \(k\) different shift interventions. Specifically, there are vectors \(\{c_i\}_{i=1}^k\) with each \(c_i \in \mathbb{R}^n\), and we observe the stationary distribution of the SDE, \[\label{eq:sde} dX_t = (v(X_t) + c_i)dt + \sqrt{\epsilon}dB_t.\tag{1}\]

We denote the concatenated intervention column vectors by the matrix \(C \in \mathbb{R}^{n \times k}\). In the causal disentanglement literature, shift interventions typically give fewer guarantees [14], [15], and even in the case of linear SDEs, the drift is trivially not identifiable without knowledge of the intervention vectors. Therefore, we assume knowledge of the interventions \(C\).

3 Related Work↩︎

3.1 Causal Representation↩︎

In terms of modeling causality [16], the most popular underlying model is the structural causal model (SCM), which characterizes the conditional distribution of a random variable under arbitrary intervention. Learning an SCM typically requires very strong assumptions such as sparsity [17], interventional data [5], parametric assumptions [18], among many other results. Sparsity is a very common theme in these models, though it may also be expressed in an assumption that the number of latent variables is small, i.e. a low-rank constraint [19].

3.2 Causal Disentanglement with Interventional Data↩︎

The bulk of the literature on causal disentanglement focuses on SCMs with an underlying DAG. Relevant to our work are results that assume access to multiple interventional environments, either acting directly on the observed variables [20] or identifying a latent model under some distributional assumptions [5], [14], [15].

Some of these works have addressed the crucial limitation of acyclicity [21]–[23], but without necessarily incorporating dynamics. Many works require hard interventions, where an intervened variable is a function of exogenous noise, although some can handle soft interventions [24].

3.3 Dynamical System Methods↩︎

Previous work has considered modeling perturbations with dynamical systems [6]. One can prove identifiability of SDE parameters from the generator [11] or trajectories [25], [26]. Other methods work in our harder setting where observed data is drawn from the stationary distribution under an SDE, and match a learned SDE under numerous interventions. These works focus on linear drift and intervention-dependent parameters [27] or non-linear drift with shift interventions [28]. Notably, neither paper gives theory to confirm if these models are identifiable.

There are also methods specific to a particular scientific domain. In genomics, some methods act on pseudotime, an inferred notion of time from cell states when very few “real” timepoints are available [29], [30]. Although harder to obtain due the destructive nature of RNA sequencing, genuine temporal data (even with very few timepoints) can also be modeled with the intent of extracting a GRN [31].

The closest work to ours proves an identifiability bound for linear SDEs under a strong sparsity assumption [32]. However, their result doesn’t consider any interventional data, and the exact pattern of sparsity in the drift matrix is assumed to be known a priori, which is rarely the case in applied problems of interest. Another closely related work is [25], which can simultaneously infer the drift and diffusion. However, this work applies only in an easier setting that assumes linearity, doesn’t study interventions, and assumes multiple temporal marginals rather than just the stationary distribution. Similarly [33] proves identifiability with interventions but assumes access to the distribution over trajectories.

4 Main Results↩︎

We consider two main parameterizations of the drift as linear or a two-layer neural network (an MLP) to verify when the model may be uniquely identified up to appropriate invariances. Prior work has considered sparsity in the linear setting only [32], but we focus on parameterizations that project to a low-dimensional space. In other words, although the ambient dimension \(n\) may be large, we assume the hidden dimension \(r\) is much smaller, at least \(n>2r\), and ideally the number of interventions to uniquely identify the parameters should scale with \(r\). This focus is driven by the empirical observation that high-dimensional dynamics are often driven by a low-dimensional subset, for example in gene modules [34].

4.1 Linear case↩︎

We start with the linear case and consider the parameterization \(v(x) = (AB - D)x\) where \(A \in \mathbb{R}^{n \times r}\), \(B \in \mathbb{R}^{r \times n}\), \(D \in \mathbb{R}^{n \times n}\). Clearly \(AB\) is a redundant parameterization of a rank \(r\) matrix, but we use this notation to contrast the non-linear setting in Section 4.2. To ensure the drift is Hurwitz and the SDE has a stationary distribution, we assume \(\|AB\| \leq \gamma < 1\) and \(D \succeq I\). Intuitively, \(AB\) drive the dynamics while \(D\) is a decay term to prevent unbounded dynamics. This parameterization is similar to the one used in [27], which also uses linear SDEs but considers hard interventions with having entirely new rows in the drift matrix, rather than shift interventions.

It is impossible to get a good deterministic guarantee if the dynamics and interventions are chosen adversarially, as seen in the following proposition:

The proof is provided in Appendix ¿sec:sec:linear95deterministic?. Intuitively, one can choose interventions that don’t affect the system dynamics at all. But this situation is pathological, and in the linear case under some weak distributional assumptions we can show identifiability.

We take care to explain why these assumptions are not particularly restrictive. The low-rank constraint is already enforced by the drift function so \(A\) and \(B\) are almost surely full-rank when drawn from any density. The spectral norm bound is necessary to guarantee the SDE has a stationary distribution. And the rotational invariance assumption encodes an uninformative prior on which low-rank subspace governs the dynamics. Two simple choices that satisfy these assumptions are sampling \(A\) and \(B^T\) from either the uniform measure on the Stiefel manifold (rectangular matrices with orthonormal columns) or with iid Gaussian entries, and then scaling down to ensure a spectral norm strictly smaller than one, with \(C\) sampled from any density.

Additionally, for theoretical tractability, we presume knowledge of the decay term \(D\). This is a stronger assumption, but plausible in some applied contexts. For example, in single-cell genomics the decay rate of genes can be estimated with external experiments (e.g., from BRIC-seq [35]).

Altogether, we can now present a nearly tight identifiability result in the linear case.

See Appendix ¿sec:sec:linear95tight? for the proof. Naively, one might count \(2nr\) unknown parameters in the entries of \(A\) and \(B\), and assume the \(n^2\) entries of the stationary covariance \(\omega\) would be enough to identify them, but this isn’t the case. One needs exactly enough interventions to account for the hidden rank.

4.1.0.1 Identifying the decay.

If diagonal decay term \(D\) is not inferred in advance, learning it simultaneously is comparable to the setting of robust PCA [36] or recovery of a diagonal plus a positive semidefinite low rank term [37]. However, unlike the usual matrix completion setting where entries of the matrix are revealed uniformly at random, we have access to correlated low-rank measurements of the form \(e_ic_j^T\) for \(\{e_i\}_{i=1}^n\) the elementary basis. This setting is well studied with sub-Gaussian measurement vectors [38], but the tools do not readily apply to our setting with non-random measurements and the constraints of the Lyapunov equation. Nevertheless, we observe in Section 5 that identifiability with an unknown diagonal decay matrix \(D\) is empirically achieved, while still subject to the lower bound established in Theorem 3.

4.2 Nonlinear case↩︎

The nonlinear case is more challenging, because there is no longer a nearly closed form characterization of the stationary distribution. To apply tools from above, we consider when a linearized SDE can approximately capture the true dynamics, by restricting to contractive drift and small noise. We consider the vector field \(v(x) = A\sigma(Bx) - x\), with the constraints that \(\|A\|, \|B\| \leq 1\).

The stronger assumptions on \(A\) and \(B^T\) are necessary for quantiative bounds. Furthermore, their columns have unit norm, which is necessary to prevent scaling issues that would lead to non-identifiability where one could move constant factors from \(A\) or \(B\) into the definition of \(\sigma\).

The strongest constraint here is condition 1, upper bounding the first derivative of the activation, as it implies the noiseless dynamics are globally contractive, but it guarantees a stationary distribution exists for any intervention vector \(c\). Furthermore, these conditions still allow for a wide class of possible activations \(\sigma\), mainly constrained to be increasing with bounded first and second derivative. We observe the improved expressiveness of the stationary distribution of nonlinear SDEs in Figure 1. Note constraint 4 and 5 rule out functions that are linear (or locally linear). This enables a stronger possible identifiability guarantee than only recovering \(A\) and \(B\) up to rotation as in the linear case.

4.2.0.1 Noiseless Setting.

One might be tempted to simply consider the noiseless case, where the SDE reduces to an ordinary differential equation (ODE), and rather than recovering the stationary distribution one recovers instead the unique global stable point. However, the noiseless case cannot take advantage of the low-rank structure, and requires \(\Omega(n)\) interventions a.s. for identifiability.

The proof is in Appendix ¿sec:sec:ode95failure?. Intuitively, the equilibrium points of the noiseless ODE approximate the mean of the SDE for small noise (note this is not true for larger noise, for example see [39]). This suggests that in order to make use of the low-rank assumption on the dynamics, one must at least take advantage of second order moments of the SDE even in the small noise limit.

4.2.0.2 Moments in the zero-noise limit.

As the noise converges to zero, the stationary distribution converges to a dirac centered on the global stable point, and there are no higher order moments. However, the rescaled second-order moments have a non-trivial limit as noise goes to zero:

The proof is given in Appendix ¿sec:sec:mean95cov?. This is one instance of perturbation theory for SDEs [40], [41]. Equipped with this fact, one can consider access to the stationary distribution of an SDE, specifically the first and second-order moments, and inspect what happens now as the noise goes to zero. Due to issues with scaling and permutation of \(A\) and \(B\) being unidentifiable from the model, we measure recovery using the signed permutation distance \(d_P\).

These polynomial dependencies treat all other problem parameters, i.e. \(\gamma, \tau, M, \nu\) as constants. Without prior knowledge of the activation, we can still guarantee identifiability but with a limiting bound instead:

We give a brief sketch of the proof, to capture the novel elements and how the identifiability of \(A\) and \(B\) appears here but not in the linear case. The zero-noise limit approximately linearizes each intervened SDE. Unlike the exactly linear case, the stationary covariance is not identical across interventions. By Theorem 5, the covariance under the \(i\)th intervention approximately matches the SDE with drift \(A(J\sigma(Bx_i^*))B\) where \(x_i^*\) is the unique zero of the intervened drift \(v(x) + c_i\).

Sufficient variability in these fixed points (guaranteed by our assumptions on \(\sigma\)) enables the identification of \(A\) and \(B\). Suppose one had access to the drift matrices: because the Jacobian term in the drift is diagonal, a linear combination of drifts across \(r\) interventions can yield a rank-one matrix, which must correspond to the outer product of a column of \(A\) and a row of \(B\). But we don’t have access to the drift matrices directly, only the covariances they induce, and so the proof relies on finding a different set of matrices. The robust version of this argument is based on the tensor decomposition analysis in [42].

We hypothesize the constraint that the SDE drift have a single absorbing point is unnecessary. If the ODE given by \(\frac{dx}{dt} = v(x) + c_i\) had multiple stable equilibria, then in the limit of small noise, the SDE stationary distribution should approach a Gaussian mixture centered around each equilibrium point. The proof only requires certain linear independence conditions of vectors induced by each fixed point, so one intervention with multiple equilibria would offer essentially the same information as additional unimodal interventions. This is consistent with the observations in [28] that nonlinear SDE models generalize well with a limited number of interventions.

5 Experiments↩︎

To validate the theory given above, we consider experiments demonstrating the recovery of SDE parameters and generalization to unseen interventions. To show the broad applicability of the theory, we consider multiple loss functions: a loss acting directly on the parameters in the linear case [27], a kernelized Stein discrepancy [43], and rollout loss for neural SDEs [44].

5.1 Linear SDE Recovery↩︎

Because we are fitting linear SDEs where the stationary distribution is Gaussian, we can choose a very simple loss that matches the population mean and covariance of each intervented distribution, to verify the claim of identifiability (Theorem 3). This loss is akin to the one used in the Bicycle method introduced in [27] that fits interventional linear SDEs. Assuming the noise scale \(\epsilon\) is known and \(L\) denotes the true drift, we fit the parameters \(\hat{A}\), \(\hat{B}\), and \(\hat{D}\). Letting our estimate for the drift be denoted by \(\hat{L} := \hat{A}\hat{B} - \hat{D}\), the loss function matches the mean and covariance of each intervention: \[\label{eq:lin95loss} \mathcal{L}_{lin}(\hat{A}, \hat{B}, \hat{D}) = \|\hat{L}^{-1}C - L^{-1}C\|_F + \|\hat{L} \omega + \omega \hat{L}^T + \epsilon I\|_F.\tag{2}\]

We evaluate the recovery of \(AB\) under different numbers of interventions for two different ambient dimensions \(n\) and two different true ranks \(r\) in Figure 2. For each sampled SDE, we use the best train error from 100 independent initializations to deal with the nonconvexity of the objective. The only exception is the oversampled \(k=r\log(n)\) where we need only 5 initializations, as the landscape is empirically easier to learn. Further details are given in Section 9.

We confirm our theory: when the decay is known, \(r-2\) interventions fails to grant identifiability with poor performance and extremely high variance, even with many independent initializations, while \(r\) interventions clearly succeed. When the decay isn’t known, we use a modest oversampling of \(r\log(n)\) interventions. This amount is informed by the literature on matrix completion, where \(\Omega(rn\log(n))\) revealed matrix entries were proven to be necessary for a rank \(r\) matrix recovery problem [45]. We scale down by a factor of \(n\) because in our setting each intervention yields the mean vector of the perturbed linear SDE in \(\mathbb{R}^n\), hence \(n\) measurements.

5.2 Nonlinear SDE Recovery↩︎

We repeat this verification of identifiability in the setting of Theorem 6, using a very small value of \(\epsilon=10^{-5}\) to approximately linearize. For some randomly chosen \(A\) and \(B\), we define the true drift \(v(x) = A\sigma(Bx) - x\) and use each intervened drift \(v + c_i\) to calculate the corresponding sample mean \(m_\epsilon^i\) and sample covariance \(\Sigma_\epsilon^i\).

We define a parameterized drift of the form \(\hat{v}(x) = \hat{A}\sigma(\hat{B}x) - x\) where \(\sigma(x) = 0.7x + \tanh(x)\) as this satisfies the Assumption 2. We train with the loss function

\[\begin{align} \mathcal{L}_{nonlin}(\hat{A}, \hat{B}) = \sqrt{\sum_{i=1}^k \left\|\hat{v}(m_\epsilon^i) + c_i\right\|^2} \\ + \sqrt{\sum_{i=1}^k \left\|J\hat{v}(m_\epsilon^i)(\Sigma_\epsilon^i / \epsilon) + (\Sigma_\epsilon^i / \epsilon) J\hat{v}(m_\epsilon^i)^T + I\right\|^2} \end{align}\]

The motivation for this loss comes from Theorem 5. The first term enforces the approximate fixed point property of the mean, and the second term enforces the approximate Lyapunov equation of the covariance. Thus we can train \(\hat{v}\) to match the first and second moments of the true stationary distribution of \(v\), which will be approximately Gaussian for sufficiently small \(\epsilon\).

To empirically evaluate identifiability, we calculate the normalized error in recovering \(A\) and \(B\) from the parameters \(\hat{A}\) and \(\hat{B}\), chosen from the run with minimum training loss across 100 independent runs. Further details are given in the Appendix. We observe that \(r^2\) interventions enable consistently good recovery of the underlying drift parameters in the known activation setting.

5.3 Synthetic Nonlinear SDE Generalization↩︎

For an evaluation of Theorem 7, we apply the insight that identifiability doesn’t require knowing \(\sigma\) and consider learnable MLPs. Learnable activations [46] have been proposed before [47] but to our knowledge not previously applied to SDE parameterizations.

To assess generalizability, we consider a loss function proposed by [28]: the kernel deviation from stationarity (KDS), or equivalently the kernelized Stein discrepancy of the stationary distribution [43]. When the parametric drift is defined as \(v_\theta(x) = A\sigma(Bx) - Dx\) and the target stationary distribution is \(\mu\), the KDS is, \[\label{eq:kds95loss} \mathcal{L}_{KDS}(\theta) = E_{x\sim\mu}E_{x'\sim\mu} \mathcal{A}_x \mathcal{A}_{x'} k(x,x').\tag{3}\] where \(k\) is a given kernel and \(\mathcal{A}_x\) is the SDE generator acting on the variable \(x\) (for more details see [28]). We use samples from the true stationary distribution \(\mu\) to approximate these expectations, further training details are given in the Appendix. This loss is nevertheless unsuitable when using noise small enough to reach the limit proven in Theorem 7, and therefore we assess generalizability according to performance on unseen interventions.

We follow [28], [48] and evaluate with the mean squared error (MSE) between the true and predicted distribution. Further details are given in Section 9.3. The results are given in Table 2. We observe a clear benefit for low noise SDEs, which gets smaller as the noise gets bigger and further away from our proven settings. This is consistent with our theory, and also intuitive: as the noise gets larger, the stationary distribution gets smoother and has less dependence on the intervention vectors. Nevertheless, in the low noise regime we see that learnable activations improve prediction on test interventions without overfitting.

5.4 Simulated GRN SDE Generalization↩︎

We consider an applied setting for learnable activations, and how they impact the recovery of GRNs. In the genomics context, \(n\) is the number of expressed genes in the thousands, and \(r\) corresponds to a much smaller number of latent gene modules/pathways. Because there are few fully known regulatory networks we rely on semi-synthetic data to produce data with similar characteristics to true perturbed transcriptomic samples and a ground-truth GRN.

We use the PerturbODE model from [31], similar to our setup and modified to handle a loss for parameterized SDEs [44]. We consider the inferred GRN as the matrix multiplication of the weight matrices in the MLP that parameterizes the drift, corresponding to a first-order Taylor approximation of the dynamics. GRN recovery can be measured as a classification task of individual edges. As a baseline, we also consider the Bicycle interventional linear SDE method [27]. To simulate a known GRN, we use the boolean ODE based-simulator BEELINE [8]. BEELINE uses the provided GRN to parametrize an SDE in the space of mRNA expression and protein expression, with physically plausible regulation dynamics. More precise details are given in the Appendix in Section 9.4.

GRN recovery is still extremely difficult (random guessing is \(\approx 0.0586\) and maximum AUPRC for dynamical methods \(\approx 0.07\)), but we observe improvement using nonlinear SDEs over linear SDEs, and substantial improvement using learnable activation in the SDE parameterization (Figure 3).

6 Discussion↩︎

We confirm identifiability directly in the linear case and non-linear, low-noise limit, and see improved generalization using learned activations for small noise. It may be beneficial to increase the capacity of the SDE beyond two-layer neural networks, if the extracted drift were still identifiable.

We note some limitations of the given theory. The results demonstrate identifiability in a setup with infinitesimal small noise and restrictions on the possible drift functions. Generalizing to larger noise regimes may be possible with higher order moment information about the stationary distribution, although understanding the stationary distribution outside of linearization is very challenging. Broader parameterizations of the drift or interventions are interesting extensions for future work.

7 Conclusion↩︎

In this work, we’ve given the first provable results regarding identifiability of interventional SDEs, in both the linear and nonlinear case. Although such models are currently less popular and less studied then comparable SCM-based modeling, the ability to obtain provable guarantees for dynamical systems without any temporal data suggests they may become more fruitful for causal inference in the future. As longitudinal genomics datasets often have a small number of time-points, future work may consider how to incorporate sparse trajectory information with the stationary distribution for better identifiability guarantees or more effective architectures for recovering regulatory networks.

8 Proofs of Results↩︎

Proof of Proposition 2↩︎

Consider block matrices such that, \[\begin{align} A &= \left[\begin{array}{ c | c } \frac{1}{\sqrt{2}} I_r & 0 \\ \hline 0 & 0 \end{array}\right] \\ B &= \left[\begin{array}{ c | c } \frac{1}{\sqrt{2}} I_r & 0 \\ \hline 0 & 0 \end{array}\right] \end{align}\] and set \(D = I\) and \(\epsilon = 1\). If we also set \[\begin{align} \omega &= \left[\begin{array}{ c | c } I_r & 0 \\ \hline 0 & \frac{1}{2} I_{n-r} \end{array}\right] \end{align}\]

then it is straightforward to confirm that the drift matrix \(L = AB - I\) satisfies the Lyapunov equation

\[\begin{align} L\omega + \omega L^T + I = 0 \end{align}\]

Now, consider any skew-symmetric matrix \(Q\) that is only supported in the top left \(r \times r\) block. Then \(AB + Q \omega^{-1}\) will only be supported in the top left block and therefore still have rank \(r\). Furthermore, \(Q\) obeys the trivial Lyapunov equation \(Q\omega + \omega Q^T = 0\), so we have that the drift matrix \(\hat{L} := AB + Q \omega^{-1} - I\) also satisfies the Lyapunov equation \(\hat{L}\omega + \omega \hat{L}^T + I = 0\). Choosing the norm of \(Q\) small enough guarantees that \(\hat{L}\) is still Hurwitz, while choosing \(Q \neq 0\) guarantees \(L \neq \hat{L}\).

In the worst case, every intervention is of the form \(e_i\) for \(i > r\). The block structure and Woodbury matrix identity imply that \((AB + Q \omega^{-1} - I)^{-1}e_i = -e_i\) for any \(Q\) chosen as above, and therefore \(-\hat{L}^{-1}e_i = -L^{-1}e_i\).

We conclude that for each intervention \(e_i\) with \(i>r\), by Theorem 1 the drift matrices \(L\) and \(\hat{L}\) both induce identical stationary distributions. Thus, even with \(n-r\) interventions, the drift is not identifiable.

Proof of Theorem 3↩︎

We show separately the upper and lower bounds for number of interventions needed for almost sure identifiability.

Proof of Proposition 4↩︎

In the zero-noise setting, the SDE reduces to an ODE, and in the limit as time goes to infinity, the stationary distribution is simply a point mass on the unique fixed point of the drift. Therefore, one only gets access to the stationary points \(x_i^*\) for each vector field \(v + c_i\): other higher moments are not defined.

Suppose we have \(n-r-1\) interventions. Since \(A\) has \(r\) columns, there must be a vector \(u\) in \(\mathbb{R}^n\) orthogonal to the columns of \(A\) and \(C\).

Let \(\hat{B} = B + bu^T\) for any non-trivial \(b \in \mathbb{R}^r\). We claim \(\hat{B}\) induces the same data as \(B\). First, note that because \(u\) is orthogonal to every \(c_i\), and \(x_i^* - c_i\) is in the image of \(A\), it follows \(u\) is also orthogonal to each \(x_i^*\).

Thus, one can confirm that \(\hat{B}x_i^* = Bx_i^*\), so the fixed points don’t change. But \(\hat{B} \neq B\), and they are clearly not equivalent up to permutation or scaling invariance.

Proof of Theorem 5↩︎

We will require two lemmas.

Using these results, we can proceed to the main perturbation theory result:

Proof of Theorem 6 and Theorem 7↩︎

We begin with a series of lemmas. We first need a result about the equilibrium points of the drift, ruling out any degeneracies. As before, we let \(x_i^*\) denote the unique zero of \(v(x) + c_i = A\sigma(Bx) -x + c_i\). Additionally, in the sequel we will consider \(\lesssim\) to absorb all constants \(\gamma, \tau, M, \nu\).

The following lemma is essentially an application of the main result in [42], rephrased in our setting:

Our last lemma requires some exposition. Consider \(A\) and \(B\) in the support as defined in Assumption 2, and diagonal \(D \in \mathbb{R}^{r \times r}\) in the support of \(\sigma'(\mathbb{R}^r)\). Let \(\omega\) be the solution to the Lyapunov equation \((ADB - I)\omega + \omega(ADB-I)^T + I = 0\). Define \(\Gamma = A^T \omega^{-1} A\).

Given diagonal matrices with free parameters \(\alpha, \beta \in \mathbb{R}^{r \times r}\), define \(X = \alpha (DBA - diag(DBA)) + \beta\).

We introduce the notation \(\sigma_{[i]}'\) to denote the Jacobian of \(\sigma\) evaluated at \(Bx_i^*\), to distinguish it from an element of \(\sigma\). Because \(\sigma\) acts elementwise, \(\sigma_{[i]}'\) is necessarily a diagonal matrix.

Finally, for ease of notation in the perturbation analysis, given matrices \(M_1, M_2\) we will write \(M_1 = M_2 + O(\epsilon)\) to mean \(\|M_1 - M_2\| = O(poly(n,r,k)\epsilon)\). We use similar notation for vectors in the \(l_2\) norm.

9 Experimental Details↩︎

9.1 Linear SDE Recovery↩︎

For the true simulated linear SDE, \(A\) and \(B\) are sampled with iid Gaussian entries, and then normalized to have spectral norm equal to \(0.9\). The true decay matrix \(D\) has each diagonal entry sampled iid from the uniform distribution on the interval \([1,2]\). The model is trained with the loss given in Equation 2 , with the Adam optimizer [53] with an initial learning rate of \(0.005\) and \(3000\) iterations. Additionally, due to the non-convexity of the objective, on each training instance we run 100 independent initializations and pick the one with the smallest training error. This excludes the oversampled setting with \(k=r * \log(n)\), where we get good performance with only 5 independent initializations.

The plotted mean and standard deviations in Figure 2 are based on 5 separate instantiations of the true model and fitting a new model on the stationary distribution. Experiments were run on CPU. All experiments (including those below) were run on a Linux system.

9.2 Nonlinear SDE Recovery↩︎

The true simulated \(SDE\) has \(A\) and \(B^T\) sampled uniformly from the Stiefel manifold, and rescaled by \(\gamma = 0.995\), with \(\epsilon=10^{-5}\) and \(\sigma(x) = 0.7x + 0.3\tanh(x)\). We use ambient dimension \(n=8\) and true low-rank dimension \(r=2\). To get training data, we sample 1000000 samples from each perturbed SDE, initialized at the stationary point with thinning of 300 and burnin of 100. We then calculate the mean and covariance for use in the loss.

In training, the learned network is trained for 10000 iterations and learning rate \(0.002\) also in Adam, using the loss \(\mathcal{L}_{nonlin}\) introduced in Section 5.2. The interventions are iid Gaussians with standard deviation of \(0.5\). We do 100 runs, choose the drift with smallest training loss, and measure the normalized Frobenius error against the true matrices \(A\) and \(B\) (after unit scaling and minimized over all column / row permutations, as the identifiability is only up to scaling and permutation), then report the mean and standard deviation across 10 independent instances.

9.3 Synthetic Nonlinear SDE Generalization↩︎

To parameterize learnable activations, given an input \(x \in \mathbb{R}^n\), we learn a function that acts elementwise on \(x\) without necessarily being the same function on each element, i.e. \(\sigma(x) = \left[\sigma_1(x_1), \dots, \sigma_n(x_n) \right]\). We choose to parameterize each \(\sigma_i\) as two-layer MLP using the actual sigmoid \(\tilde{\sigma}\) as the activation function. As a warm start, we also add \(0.1\tilde{\sigma}(x)\).

We consider a fixed true SDE, with hidden dimension \(r = 3\), where the rectangular matrices \(A\) and \(B^T\) have one down the diagonal and zero elsewhere, \(D = I\), \(\gamma = 0.98\), and the elementwise activation function is given by \[\begin{align} \sigma_1(x) &= 3\cos(3(x-0.5))\\ \sigma_2(x) &= 2\sin(2(x+1.5)) - 1\\ \sigma_3(x) &= \sigma_1(x) \end{align}\] The hidden dimension of the activation MLPs is fixed at 20.

Each intervention is sampled as a Gaussian vector with mean zero and variance \(0.1\) on each entry, and from each intervened SDE we collect 5000 samples. To do this, we use the Euler-Maruyama scheme to solve the SDE, and use MCMC to draw approximately independent samples from the stationary distribution, using \(dt = 0.01\), a thinning factor of \(300\) and the first 500 samples treated as burn-in. We use an radial basis function kernel to parameterize the KDS loss.

We train with 10 interventions, and evaluate on 20 held-out interventions. Training is done with AdamW [54], with initial learning rate \(0.003\) and 50000 iterations. The hyperparameter ranges considered are given in Table 3. All experiments were run on CPU.

9.4 Simulated GRN SDE Generalization↩︎

We first give a summary of the simulator. BoolODE is scRNA-seq data simulator that samples mRNA readouts via a stochastic differential equation (SDE). The underlying cyclic GRN is encoded in the SDE via a Hill function based approximation to a boolean logical circuit. \(x_i\), \(p_i\), and \(R[i]\) each represents the mRNA concentration for gene \(i\), the concentration of the protein encoded by gene \(i\), and the set of proteins that regulate gene \(i\). Further, \(r\), \(l_p\), \(l_x\), and \(s\) are scalar coefficients representing translation rate, RNA decay rate, protein decay rate, and noise standard deviation respectively. The functions \(f_i(\cdot)\) encodes the gene regulatory network. The default SDE is given by, \[\begin{align} \frac{dx_i}{dt} &= f_i(p_{R[i]}) - l_x x_i + s \sqrt{x_i} dB_t \\ \frac{dp_i}{dt} &= r x_i - l_p p_i + s \sqrt{p_i} dB'_t \end{align}\] where \(B\) and \(B'\) are independent Brownian motions.

We simulate overexpression (e.g., CRISPR-a) experiments by inducing perfect intervention coupled with increased transcription for a set of intervened genes \(I_k\). For each intervened gene \(j \in I_k\) we set \(f_j(\cdot)=0\) and add a positive shift (set to 20 in the experiments) to \(dx_j/dt\). We consider \(K\) such intervention regimes, plus the observational setting (\(k=0\)) with no intervention, i.e, \(I_0=\emptyset\), for a total of \(K+1\) regimes. With this setup, we obtain a family of distributions \(\{\rho_k\}_{k=0}^K\) in gene expression space.

We adapt the neural ODE architecture proposed in [31]. The base SDE model under \(I_k\) parametrizes the change of gene expression via the drift function as below. \(A\) and \(B\) are the coefficient matrices encoding the modular graph. \(\alpha\) is a scaling vector controlling the rate of non-linear activation in the module, where the activation \(\sigma\) is the logistic sigmoid. Further, \(\beta\) is a bias term that shifts the activation threshold of the modules. The GRN is extracted as \(A \, \text{diag}(\alpha \circ \mathbf{1}_N) B\). The intervention is modeled as a combination of standard basis vectors, \(\sum_{j \in I_k} e_{j}\), specifying the overexpression. The model also learns a single diffusion coefficient scalar. Lastly, \(M\) is a masking matrix blocking the signals from the intervened genes’ regulators. Altogether, the drift is given by, \[\begin{align} v_k(x) = M_{k} A \sigma( \alpha \circ (B x - \beta) ) + \sum_{j \in I_k} e_{j} - Dx. \end{align}\] We alternatively consider the architecture using a learnable MLP to model the nonlinear activation of modular signals, \[\begin{align} v_k(x) = M_{I_k} A \sigma_* (B x - \beta) + \sum_{j \in I_k} e_{j} - Dx \end{align}\] where \(\sigma_*\) denotes the MLP as described above. The model loss is the Sinkhorn divergence [55] applied to the true perturbed distribution and the samples drawn from the learned SDE.

Each interventional distribution \(\rho_k\) is obtained by taking the observational data \(\rho_0\) as initial distribution and simulating the SDE via the Euler-Maruyama method.

We fit the hyperparameters by testing the ODE architecture in [31] on the same simulated datasets, and then doing zero-shot transfer of the hyperparameters to the SDE models. For Bicycle, we sweep the hyperparameters in Table 4 and use the package defaults for the remaining values. Experiments are run on an nvidia tesla V100 gpu.

References↩︎