Science, Faculty of
Mathematics, Department of
DSpace
UBCV
Rahmani, Hossein
Master of Science - MSc
University of British Columbia
Stochastic differential equations (SDEs) are a fundamental tool for modelling dynamic processes, including gene regulatory networks (GRNs), contaminant transport, financial markets, and image generation. However, learning the underlying SDE from data is a challenging task, especially if individual trajectories are not observable. Motivated by research in single-cell datasets, we present the first comprehensive approach for jointly identifying the drift and diffusion of an SDE from its temporal marginals. First, a detailed explanation of the trajectory inference problem is given, and its connection to entropic optimal transport formulated as a convex variational problem is established. Assuming linear drift and additive diffusion, it is then proved that these parameters are identifiable from marginals if and only if the initial distribution lacks any generalized rotational symmetries. Further, it is shown that the causal graph of any SDE with additive diffusion can be recovered from the SDE parameters, thereby linking the problem to causal discovery. To complement this theory, entropy-regularized optimal transport is adapted to handle anisotropic diffusion, and APPEX (Alternating Projection Parameter Estimation from X0) is introduced as an iterative algorithm designed to estimate the drift, diffusion, and trajectories solely from temporal marginals. It is demonstrated that APPEX iteratively decreases KullbackâLeibler divergence to the true solution, and its effectiveness is illustrated on simulated data from linear additive noise SDEs.
Identifying drift, diffusion, and causalstructure from temporal snapshotsbyHossein RahmaniB.Sc., Sharif University of Technology, 2023A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)October 2025© Hossein Rahmani 2025The following individuals certify that they have read, and recommend tothe Faculty of Graduate and Postdoctoral Studies for acceptance, the thesisentitled:Identifying drift, diffusion, and causal structure from temporalsnapshotssubmitted by Hossein Rahmani in partial fulfillment of the requirementsfor the degree of Master of Science in Mathematics.Examining Committee:Elina Robeva, Associate Professor, Department of Mathematics, UBCSupervisorGeoffrey Schiebinger, Associate professor, Department of Mathematics, UBCSupervisory Committee MemberiiAbstractStochastic differential equations (SDEs) are a fundamental tool for modellingdynamic processes, including gene regulatory networks (GRNs), contami-nant transport, financial markets, and image generation. However, learningthe underlying SDE from data is a challenging task, especially if individualtrajectories are not observable. Motivated by research in single-cell datasets,we present the first comprehensive approach for jointly identifying the driftand diffusion of an SDE from its temporal marginals. First, a detailed ex-planation of the trajectory inference problem is given, and its connectionto entropic optimal transport formulated as a convex variational problem isestablished. Assuming linear drift and additive diffusion, it is then provedthat these parameters are identifiable from marginals if and only if the ini-tial distribution lacks any generalized rotational symmetries. Further, itis shown that the causal graph of any SDE with additive diffusion can berecovered from the SDE parameters, thereby linking the problem to causaldiscovery. To complement this theory, entropy-regularized optimal trans-port is adapted to handle anisotropic diffusion, and APPEX (AlternatingProjection Parameter Estimation from X0) is introduced as an iterative al-gorithm designed to estimate the drift, diffusion, and trajectories solely fromtemporal marginals. It is demonstrated that APPEX iteratively decreasesKullbackâLeibler divergence to the true solution, and its effectiveness is il-lustrated on simulated data from linear additive noise SDEs.iiiLay SummaryMany important systems in nature and society cannot be fully observed, andwe must rely on partially observed data to understand them. Examples in-clude how genes regulate each other in a single cell or how financial marketsevolve. This research develops mathematical tools to uncover these hiddendynamics using only snapshots of the system taken at different times. Theresults show how the underlying forces that drive change, the role of ran-domness, and the network of causal relationships between variables can beidentified. A new method, called APPEX, is introduced to estimate thesehidden processes step by step, using ideas from optimal transport. Experi-ments on simulated data confirm that the method can reliably reconstructboth the evolution and the structure of complex systems, making it appli-cable to fields such as biology, physics, and finance.ivPrefaceThe content of this thesis is based on the joint work of the author with hissupervisor, Elina Robeva, and collaborators Vincent Guan, Joseph Janssen,Andrew Warren, Stephen Zhang, and Geoffrey Schiebinger on âIdentifyingDrift, Diffusion, and Causal Structure from Temporal Snapshots.â A versionof this work has been submitted for publication and is available at arXiv.Chapter 1 was written by the author to provide background and a detailedexplanation of the trajectory inference problem and its connection to en-tropic optimal transport. Chapters 2â6 and the Appendix are based on thesubmitted paper, with only minimal changes to the introduction, the de-scription of the parameter estimation method (APPEX), and the abstract.The contributions of the authors to this work are as follows:â¢ Project conception: The project idea originated with Elina Robeva.The iterative refinement strategy for trajectory inference and parame-ter estimation was developed by the author and Geoffrey Schiebinger.â¢ Theory: Identifiability results in Chapter 3 were proved primar-ily by Vincent Guan, Joseph Janssen, and Elina Robeva, with ad-ditional consultation from Andrew Warren, Stephen Zhang, and Ge-offrey Schiebinger.â¢ Algorithm design: The APPEX algorithm and the argument foriterative KullbackâLeibler divergence decrease were developed jointlyby all authors, with significant contributions from the author to thealgorithmâs design.â¢ Experiments: Numerical experiments were carried out in Python bythe author and Vincent Guan, and in R by Joseph Janssen.â¢ Writing: All authors contributed to writing the manuscript [1]. Theauthor also wrote Chapter 1 (background) and adapted parts of Chap-ters 2â6 and the Appendix for the thesis.vPrefaceThe use of generative AI in this thesis was limited to using ChatGPT-5 andGoogle Gemini to verify the grammatical correctness and clarity of sentencesand paragraphs. These tools were not used for producing mathematicalcontent, proofs, algorithms, or experimental results.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Stochastic Differential Equations (SDEs) . . . . . . . . . . . 11.1.1 Inference from temporal marginals . . . . . . . . . . . 21.2 Optimal Transport (OT) . . . . . . . . . . . . . . . . . . . . 21.2.1 Discrete and Continuous Formulations . . . . . . . . 31.2.2 Entropic Regularization and Relative Entropy . . . . 41.3 Trajectory Inference via Schro¨dinger Bridges . . . . . . . . . 51.4 Identifiability of Linear additive noise SDE parameters . . . 62 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2 Our contributions . . . . . . . . . . . . . . . . . . . . . . . . 93 System identification from temporal marginals . . . . . . . 113.1 SDE parameter identifiability . . . . . . . . . . . . . . . . . . 123.2 Causal graph identification from SDE parameters . . . . . . 18viiTable of Contents4 Our parameter estimation method . . . . . . . . . . . . . . . 214.1 Trajectory inference via anisotropic entropy-regularized opti-mal transport . . . . . . . . . . . . . . . . . . . . . . . . . . 244.2 Parameter estimation via MLE . . . . . . . . . . . . . . . . . 254.3 The APPEX algorithm . . . . . . . . . . . . . . . . . . . . . 275 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 295.1 Revisiting non-identifiable SDEs . . . . . . . . . . . . . . . . 305.2 Higher dimensional randomly generated linear additive noiseSDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 305.3 Causal discovery . . . . . . . . . . . . . . . . . . . . . . . . . 336 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 376.1 Limitations and future work . . . . . . . . . . . . . . . . . . 38Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41AppendicesA Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51A.1 Non-identifiable linear additive noise SDEs . . . . . . . . . . 51A.2 Properties of the Linear Additive Noise SDE . . . . . . . . . 53A.3 Weak formulation of the Fokker-Planck equation . . . . . . . 54A.4 Properties of auto-rotationally invariant r.v.s . . . . . . . . . 55A.5 Additional proofs for Theorem 3.3 . . . . . . . . . . . . . . . 57A.6 SDEs and causality . . . . . . . . . . . . . . . . . . . . . . . 60A.7 Entropy-regularized Optimal Transport . . . . . . . . . . . . 63A.8 Maximum Likelihood Estimation . . . . . . . . . . . . . . . . 67A.9 Additional experiments . . . . . . . . . . . . . . . . . . . . . 72A.9.1 Consistency experiments . . . . . . . . . . . . . . . . 72viiiList of Tables5.1 Mean absolute error (MAE) of estimated drift and diffusionfor dimensions 3-10. . . . . . . . . . . . . . . . . . . . . . . . 325.2 Correlation between estimated and true drift and diffusionfor dimensions 3-10. . . . . . . . . . . . . . . . . . . . . . . . 325.3 Average Structural Hamming Distance (lower is better) withvarying dimensions and random edge probabilities p. . . . . 355.4 Average Structural Hamming Distance (HD) for simple edges,and latent confounders (lower is better) with varying dimen-sions and random edge probability p = 0.25. . . . . . . . . . 35ixList of Figures2.1 Outline of our theoretical (a) and algorithmic (b) contribu-tions. In (b), our algorithm Alternating Projection ParameterEstimation from X0 alternates between Anisotropic Entropy-regularized Optimal Transport (AEOT) and Maximum Like-lihood Estimation (MLE). . . . . . . . . . . . . . . . . . . . . 93.1 Four examples of auto-rotationally invariant distributions areillustrated. By Theorem 3.3, these initial distributions leadto non-identifiability of linear additive noise SDEs from tem-poral marginals. . . . . . . . . . . . . . . . . . . . . . . . . . 134.1 Visualization of our parameter estimation algorithm APPEX.Given observed temporal marginals from an underlying SDE(left), APPEX alternates between trajectory inference (mid-dle) and MLE parameter estimation (right) in order to findthe SDE parameters that best represent the temporal marginalobservations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 225.1 The mean absolute error for estimates of A and H usingAPPEX is shown per iteration for all three pairs of SDEsfrom Section A.1 . . . . . . . . . . . . . . . . . . . . . . . . . 315.2 The correlation between the estimated and true SDE param-eters is plotted per iteration across 10 random linear additivenoise SDEs for dimensions 3 and 10 . . . . . . . . . . . . . . 335.3 The true and estimated causal graphs by WOT and APPEXfor four random SDEs of dimensions are illustrated. Rededges i â j represent negative edge weights, such that X(i)negatively regulates X(j). Green edges i â j represent pos-itive edge weights, such that X(i) positively regulates X(j).Dashed edges represent the effect of a latent confounder Ui,jon observed variables i, j through correlated diffusion. . . . . 36xList of FiguresA.1 Marginals are plotted at various times for the SDE pairsfrom Section 3 of the Appendix, following various non-auto-rotationally invariant initial distributions. . . . . . . . . . . . 53A.2 The three types of auto-rotational invariance is illustrated.Each type leads to non-identifiability from marginals, wheninitialized at the distribution. . . . . . . . . . . . . . . . . . . 58A.3 The mean squared error of the classical MLE drift estimator(Classical) converges at a slower rate compared to the MLEdrift estimator (4.11) derived for multiple trajectories (Tra-jectory). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72A.4 Average mean squared error (MSE) of A and GGâ¤ as thenumber of observations per marginal increases across 50 ran-dom 2d systems. Note this is a log-log plot, and the dottedline in the APPEX plot shows that the convergence follows apower law relationship with power â0.125. . . . . . . . . . . . 72A.5 Average mean squared error (MSE) of A and GGâ¤ as thenumber of observations per marginal increases across 70 ran-dom 12d systems. Note this is a log-log plot, and the dottedline in the APPEX plot shows that the convergence follows apower law relationship with power -1. . . . . . . . . . . . . . 73xiList of Symbolsb(Xt) Drift function of a general SDEA Drift matrix of a linear SDEG Diffusion matrix of an additive noise SDEH = GGâ¤ Observational diffusion matrixpt â¼ Xt Temporal marginal distribution of the SDE at time tp(x, t) Probability density of the process at time t, assuming existenceXi = Xti Distribution of SDE at ith observation timeâXi Xi+1 âXipi â¼ Xti Temporal marginal distribution of the SDE at measured time tiN Number of observed temporal marginalspËi Empirical temporal marginal of the SDE at measured time tix(j)tijth sample from ith empirical marginalMi Number of samples for the ith empirical temporal marginal pËiÎ£ Covariance of an initial distribution X0ϵ Scalar parameter for regularizing entropy-regularized optimal transport(EOT)Σ(θ) (Possibly) anisotropic matrix parameter for regularizing anisotropicentropy regularized optimal transport (AEOT)xiiAcknowledgementsThere are many people to whom I am thankful. First and foremost, I amdeeply grateful to my supervisor, Professor Elina Robeva, for her supportand guidance throughout my masterâs studies. From her I learned not onlycausal inference but also how to approach research. Her encouragementand commitment were essential for my growth and made this program anexciting experience.I am also grateful to Professor Geoffrey Schiebinger for his valuable guidanceand for the many insightful discussions we had, particularly on topics inoptimal transport. His perspectives and expertise greatly enriched my workand helped me refine the direction of my research.I would like to thank my collaborators and coauthors Vincent Guan, JosephJanssen, Andrew Warren, and Stephen Zhang for their generous collabora-tion and the many stimulating discussions we shared. Working alongsidethem has been an invaluable learning experience, and I have greatly bene-fited from their ideas, insights, and examples.Finally, I thank my family for their constant support, encouragement andbelief in me. Their love and care made this journey possible.xiiiChapter 1BackgroundThis chapter reviews the mathematical foundations for statistical inferencein the âmarginal onlyâ observational setting, commonly referred to as thetrajectory inference problem and gives literature review on identificationthe parameters of the underlying process that generated these trajectories.Section 1.1 introduces stochastic differential equations (SDEs). Section1.2 reviews the theory of optimal transport and its entropic regulariza-tion, further presents how the entropic optimal transport problem can beseen as minimizing relative entropy. This connection naturally leads to theSchro¨dinger bridge formulation, which is presented in Section 1.3 as a coreproblem to be solved for trajectory inference. Finally, Section 1.4 shifts fo-cus from the trajectories to the system itself, discussing the critical challengeof identifiability of SDE parameters from marginal data alone and reviewingthe assumptions required by current methods.1.1 Stochastic Differential Equations (SDEs)The mathematical framework of SDEs, or drift-diffusion processes, origi-nates from statistical mechanics. The drift models the velocity of particlesflowing in a carrier fluid, while diffusion is driven by Brownian motion [2],and governs their dispersion due to random collisions [3]. Formally, a d-dimensional time-homogeneous SDE is driven by an m-dimensional Brown-ian motion Wt, and is parameterized by a drift function b(x) : Rd â Rd anda diffusion function Ï(x) : Rd â RdÃm, such thatdXt = b(Xt) dt+ Ï(Xt) dWt, X0 â¼ p0, (1.1)where b(x) and Ï(x) are Lipschitz to ensure the solutionâs existence anduniqueness [4].In this section, we consider SDEs with additive noise.dXt = b(Xt) dt+G dWt, X0 â¼ p0. (1.2)11.2. Optimal Transport (OT)We also refer to H = GGâ¤ â RdÃd as the (observational) diffusion, sinceit is only possible to observe H [5, 6]. Time-homogeneous linear additivenoise SDEs are of particular interest [5, 7, 8], as they generalize the popularOrnstein-Uhlenbeck process by removing the assumption of stationarity. Atime-homogeneous linear additive noise SDE obeys the formdXt = AXt dt+G dWt, X0 â¼ p0, A â RdÃd. (1.3)We overview important properties of this SDE in Section A.2 of the ap-pendix, including its closed form solution and transition density.1.1.1 Inference from temporal marginalsThe marginals-only observational setting dates back to at least Aristotle,who dissected bird eggs at various developmental stages, in order to inferthe processes and mediators (e.g., heat) behind embryogenesis [9â11]. If theunderlying system is modeled by an SDE, then the evolution of its temporalmarginals pt are defined by its Fokker-Planck equation. For example, astraightforward computation shows that the marginals of a linear additivenoise SDE evolve according to the following Fokker-Planck equation:ââtp(x, t) = ââ · ((Ax)p(x, t)) + 12â · (Hâp(x, t)), p(x, 0) = p0(x). (1.4)Contemporary research in statistical inference from temporal marginals islargely based on single-cell biology. Many works focus on trajectory infer-ence for cell fates [12â15], while others infer the causal graph, e.g. GRN[16â19], or perform parameter estimation [20, 21]. Recent works jointly esti-mate a subset of these quantities by leveraging their intrinsic relationships.For instance, [22, 23] are primarily interested in trajectory inference, and it-eratively estimate drift and trajectories. Meanwhile, [8] iteratively estimatesdrift and trajectories, while additionally applying permanent interventionson the drift dynamics for network inference.1.2 Optimal Transport (OT)Optimal Transport (OT) was first proposed by Gaspard Monge in the 1780sto efficiently relocate piles of earth [24]. Today, OT is used in diverse areasof science and engineering, from fluid mechanics [25] to climate science [26]to pure math [27]. In the following section, we provide a brief introductionto the theory of optimal transport. We refer the reader to [28] for a detailedintroduction of optimal transport and its computational aspects.21.2. Optimal Transport (OT)1.2.1 Discrete and Continuous FormulationsThe theory of optimal transport originated with the Monge problem,which seeks the most efficient way to move a distribution of mass froma set of source locations {x1, . . . , xn} in Rd with corresponding masses{a1, . . . , an} to a set of targets {y1, . . . , ym} that are to receive masses{b1, . . . , bm}. The cost of moving a unit of mass from xi to yj is given by acost function c(xi, yj). In its discrete form, the goal is to find a deterministicmapping, or pushforward map S, that transports mass while minimizing thetotal transportation cost [29]:minSnâi=1c(xi, S(xi))over all S : {x1, x2, ..., xn} â {y1, y2, ..., ym} such thatâj â {1, 2, ...,m}, bj =âi:S(xi)=yjaiHowever, a solution to the Monge problem may not exist if, for example, asingle source must be split and sent to multiple destinations [24].To overcome this, Leonid Kantorovich developed a relaxation of the problemthat allows for the splitting of mass [30]. The Kantorovich problemseeks a transport plan (or coupling) P , which can be thought of as a jointprobability distribution, that minimizes the total cost. In the discrete case,this is a linear programming problem:minPâÎ (a,b)nâi=1mâj=1Ci,jPi,jwhere C â RnÃm is the cost matrix, with Ci,j denoting the cost of trans-porting one unit of mass from xi to yj , and P â RnÃm+ is the transport plan,where Pi,j represents the amount of mass moved from xi to yj . The feasibleset isÎ (a, b) := {P â RnÃm+ : P1m = a and Pâ¤1n = b}.These concepts extend to the continuous setting for probability measuresµ and ν on metric spaces. The continuous Kantorovich problem isformulated as:minÎ³âÎ (µ,Î½)â«XÃYc(x, y)dγ(x, y)31.2. Optimal Transport (OT)where Î (µ, Î½)Î (µ, ν) := {γ â P (X Ã Y) : γ|X = µ and γ|Y = ν}.Unlike the Monge problem, the Kantorovich formulation guarantees the ex-istence of a solution under general conditions [31](see Theorem 4.2).1.2.2 Entropic Regularization and Relative EntropyA crucial modification to classical optimal transport is the addition of en-tropic regularization, which provides key benefits for theory and compu-tation [32].â¢ Algorithmic: The problem becomes strongly convex and can besolved much faster than the standard linear program version.â¢ Modeling: It enables probabilistic modeling of biologically importantphenomena, such as âpartially fatedâ cells, by producing soft transportplans instead of deterministic matchings.â¢ Theoretical: It provides a rigorous foundation for trajectory infer-ence by creating a deep connection to Brownian motion.The ϵ-entropy regularized Optimal Transport problem for continuous mea-sures with densities f and h is formulated as finding a transport planγ â Î (f, h) that solves:minÎ³âÎ (f,h)â«RdÃRdc(x, y)dγ(x, y) + ϵH(γ)where H(γ) is the entropy of the plan γ and ϵ > 0 is the regularizationstrength. Note that this definition of entropy, H(γ) =â«dÎ³ log dγ, is thenegative of the standard Shannon entropy.This regularized problem is deeply connected to the concept of relativeentropy, or Kullback-Leibler (KL) divergence. The relative entropy of ameasure µ with respect to another measure ν (where µ is absolutely contin-uous with respect to ν, denoted µâª ν) is defined as:H(µ|ν) =â«logdµdÎ½(y)dµ(y)The entropy-regularized Optimal Transport objective can be shown to beequivalent to minimizing the relative entropy of the transport plan γ with41.3. Trajectory Inference via Schro¨dinger Bridgesrespect to a carefully chosen reference measure. Let us define a joint densityKϵ (Gibbs kernel) based on the cost function c(x, y):Kϵ(x, y) = f(x) 1cÏµexp(âc(x, y)ϵ)where cϵ is a normalizing constant. The relative entropy of a valid transportplan γ with respect to Kϵ can be expanded as follows:H(γ|Kϵ) =â«dÎ³(x, y) log(dγ(x, y)dKϵ(x, y))= H(Î³)ââ«dÎ³(x, y) log(Kϵ(x, y))= H(Î³)ââ«dÎ³(x, y) log(f(x)) +1Ïµâ«dÎ³(x, y)c(x, y)â log(cÏµ)=1Ïµ(â«dÎ³(x, y)c(x, y) + ÏµH(Î³))âH(f)â log(cϵ)Since f and cϵ are given, the termsH(f) and log(cϵ) are constant with respectto the optimization variable γ, minimizing the term in the parenthesis isequivalent to minimizing the relative entropy H(γ|Kϵ). Thus, the entropy-regularized OT problem is mathematically equivalent to solving:minÎ³âÎ (f,h)H(γ|Kϵ)This equivalence provides a powerful theoretical link between optimiza-tion, probability, and information theory, rooted in the classical Schro¨dingerbridge problem, and it underpins many modern computational approachesto optimal transport.1.3 Trajectory Inference via Schro¨dinger BridgesThe Trajectory Inference problem asks whether the full path-measure of adynamic process can be recovered given only its marginal distributions {Pt}at various time points.We often model these dynamics with a stochastic differential equation (SDE):dXt = b(Xt, t) dt+ ϵ dBtwhere Xt is the state, b is a drift field, and ϵdBt represents Brownian motion.In general, the marginals {Pt} do not uniquely determine the law P ontrajectories. However, a unique solution exists in the important case wherethe drift is conservative (b = ââÎ¦), as formalized in the following theorem.51.4. Identifiability of Linear additive noise SDE parametersTheorem 1.1. Let P be the law of the SDEdXt = ââÎ¦(Xt, t)dt+ ϵdBtwith initial condition H(P0|W ϵ0) < â. If R is any other law on paths withthe same marginals (Rt = Pt for all t), thenH(P |W ϵ) â¤ H(R|W ϵ)with equality if and only if P = R [13].This theorem provides a variational principle known as the Schro¨dingerBridge problem[15, 21]:For a pair of marginals, given the transition law K of a reference Wienerprocess, an initial distribution µ, and a final distribution ν, the Schro¨dingerBridge solution is the law of the stochastic process, which satisfies bothmarginal constraints, while minimizing relative entropy to K:Ïâ = ProjKLÎ (µ,ν)(K) = argminÏâÎ (µ,Î½)DKL(Ïâ¥K). (1.5)This reframes trajectory inference as a convex optimization problem, pro-viding a rigorous foundation for reconstructing dynamics from static data.In practice, this problem is solved by reformulating it as an entropic optimaltransport problem between the discrete-time marginals. This resulting con-vex optimization can then be solved efficiently using methods such as theSinkhorn algorithm [29, 32]. More technical details are available in A.7.1.4 Identifiability of Linear additive noise SDEparametersThe problem of non-identifiability occurs if there exists an alternative drift-diffusion pair, e.g. (AË, HË) Ì¸= (A,H) for linear-additive noise SDEs, whichshares the same marginals. Equivalently, the processes satisfy the sameFokker-Planck equation across all observed times, given initial distributionp0. Non-identifiability can occur due to many factors. For example, distinctSDEs can share the same stationary distribution, and would hence be non-identifiable from marginals if p0 is initialized at equilibrium [13, 21]. Also,if the initial distribution is standard normal p0 â¼ N (0, Id), then it wouldbe impossible to distinguish between a pure Brownian motion, and an SDE61.4. Identifiability of Linear additive noise SDE parameterswith an additional rotational vector field [21]. We overview three classicalexamples of non-identifiability in Section A.1 of the Appendix.Because of non-identifiability, previous works require assumptions about theunderlying SDE to perform inference. The standard assumptions are thatthe drift is an irrotational vector field âÏ [13, 15, 22, 33, 34] and that thediffusion is both known and given by isotropic Brownian motion ÏdWt[8, 13,15, 22, 23, 33, 34]. It is also common to leverage additional perturbationaldata from interventions [8, 19], and to impose regularization [8, 35], to betterlearn the system. However, these constraints impose unrealistic conditionson the data. An irrotational drift field cannot model negative feedbackloops or repressilator dynamics [33], which are common in biological data.Similarly, isotropic diffusion cannot model variable noise scales or correlatednoise structures, which commonly arise in real data [19, 36, 37]. While mostworks use a known diffusion estimate to infer drift, [38] infers diffusion, butrequires a drift estimate. In both cases, the accuracy of the prior estimate iscrucial, because misspecified diffusion typically leads to poor drift estimationand vice-versa [39]. This is problematic, since knowing drift or diffusiona priori is unrealistic in practice. Finally, we note that concurrent workestimates drift and diffusion, but does not consider the uniqueness of learnedparameters, due to non-identifiability [34, 40].To advance the field of statistical inference from temporal marginals, wework to remove previous assumptions by determining conditions for iden-tifiability and designing a method that leverages them. [5] showed thatidentifiability of a linear additive noise SDE from ground truth trajectorieswith fixed X0 â Rd is equivalent to a non-degenerate rank condition basedjointly on X0, the linear drift A and the observational diffusion GGâ¤. Incontrast, our work studies identifiability of linear additive noise SDEs in themore general âmarginals-onlyâ setting. We also demonstrate that identifi-ability is guaranteed across all linear additive noise SDEs solely from theproperties of the initial distribution. To the best of our knowledge, this hasnot been shown even in the case of trajectory-based observations.7Chapter 2Introduction2.1 MotivationThis work presents the first comprehensive approach for jointly identify-ing the drift and diffusion of a stochastic differential equation (SDE) fromobserved temporal marginals. SDEs provide a powerful framework for mod-eling complex systems influenced by both deterministic forces and randomfluctuations, captured by drift and diffusion terms, respectively [41]. Infer-ring these parameters is crucial for understanding the underlying dynamicsof a system.While parameter estimation has been studied extensively using trajectorydata, either from one long trajectory [41â48], or multiple short trajectories[5, 49â52], it is often impossible to observe individual trajectories whentracking large populations. For example, in single-cell RNA sequencing(scRNA-seq) datasets, the destructive nature of sequencing technologies pre-vents specific cells from being tracked over time [12, 13, 53â55]. The recentexplosion of scRNA-seq data has renewed interest in trajectory inferenceto model developmental processes, often conceptualized by Waddingtonâsepigenetic landscape, where SDEs can describe how cells navigate a high-dimensional gene expression space toward different fates [12, 56]. Similarly,in hydrogeology, sensors cannot track distinct contaminant particles, pro-viding only data on overall plume migration [57â64]. In such cases, observa-tional data are limited to samples from the populationâs temporal marginalsat various measurement times, a setting we refer to as marginals-only.However, the marginals-only setting presents a fundamental challenge whichis non-identifiability. Many stochastic processes are non-identifiable fromtheir marginals alone, since different underlying dynamics can produce iden-tical population-level distributions over time [33, 65]. To circumvent thisissue, previous works often simplify the problem by assuming that eitherthe drift or diffusion is already known [8, 13, 15, 22, 23, 38]. Furthermore, it82.2. Our contributionsis common to impose restrictive modeling assumptions, such as irrotationaldrift and isotropic diffusion, which may not hold in many real-world systems[8, 13, 15, 22, 23, 33]. This thesis addresses these limitations by establishingconditions for identifiability and providing a practical algorithm for joint in-ference without such constraints. This thesis is a minimally modified versionof what has been prepared for publication and is available on Arxiv.2.2 Our contributionsWe prove that for a large class of initial distributions, the drift, diffusion, andunderlying causal graph of an SDE with linear drift and additive diffusioncan be identified from its temporal marginals. We also introduce a practicalalgorithm for estimating each of these quantities solely from the marginalsof a general additive diffusion process, and demonstrate the effectiveness ofour method on simulated data. We illustrate our theoretical contributionsin Figure 2.1a and our algorithmic contributions in Figure 2.1b. A detailedoverview of our contributions is also given below.Temporal marginalsobserved froma non-auto-rotationallyinvariant X0Ground truthSDE parametersCausal graph GTheorem 3.3Lemma 3.5(a) Theoretical contributions: iden-tifying drift, diffusion, and causalstructure from temporal marginalsEstimate ofSDE parametersEstimate oflaw on pathsAEOTMLE(b) Algorithmic contributions:APPEX algorithmFigure 2.1: Outline of our theoretical (a) and algorithmic (b) contributions.In (b), our algorithm Alternating Projection Parameter Estimation fromX0 alternates between Anisotropic Entropy-regularized Optimal Transport(AEOT) and Maximum Likelihood Estimation (MLE).1. We provide a theoretical foundation for system identification from tem-92.2. Our contributionsporal marginals in Section 3.(a) In Theorem 3.3, we fully characterize the conditions for guar-anteeing identifiability of linear drift and additive diffusion fromobserved marginals. Identifiability is guaranteed if and only if theinitial distribution X0 is not auto-rotationally invariant (Defini-tion 3.2).(b) In Lemma 3.5, we connect parameter estimation to causal struc-ture learning via dynamic structural causal models [66]. We showthat the causal graph of any additive noise model can be recov-ered from the SDE parameters.2. In Section 4, we introduce our algorithm: Alternating Projection Pa-rameter Estimation from X0 (APPEX). It is the first method designedto estimate drift and diffusion from temporal marginals of additivenoise SDEs, without prior knowledge.(a) We show that with each iteration, APPEXâs estimates approachthe true solution. APPEX alternates between a trajectory in-ference step, which is a Schro¨dinger bridge problem, and a pa-rameter estimation step, which is solved via maximum likelihoodestimation.(b) To solve the Schro¨dinger bridge problem for trajectory inference,we adapt entropic optimal transport for parameter dependenttransport costs, including arbitrary anisotropic diffusion.(c) We test APPEXâs efficacy across a wide range of experiments inSection 5. We first use APPEX to empirically demonstrate thatclassical non-identifiable linear additive noise SDEs are resolvedwhen X0 is initialized appropriately according to Theorem 3.3.We then demonstrate APPEXâs effectiveness over a large-scale ex-periment on higher dimensional SDEs with randomly generatedparameter sets. Our results demonstrate that APPEX identi-fies arbitrary linear additive noise SDEs with significantly higheraccuracy than the widely-used Waddington-OT (WOT) method[12]. Our final experiments demonstrate that APPEX can beused to identify the systemâs causal graph, even in the presenceof correlated diffusion from latent confounders.10Chapter 3System identification fromtemporal marginalsWe present our main theorem in Section 3.1. We prove that generalizedrotational symmetries in the initial distribution X0 provide necessary andsufficient conditions for identifying the drift and diffusion parameters oflinear additive noise SDEs from its time marginals. We then show that theunderlying causal graph of a general additive noise SDE can be recoveredfrom its drift and diffusion parameters in Section 3.2.We first define the desired rotational properties. Recall that classical rota-tions in Rd for d â¥ 2 are given by matrix exponentials eAt, such that t is therotation angle and A is skew-symmetric, i.e. A+Aâ¤ = 0. We may similarlydefine rotations over anisotropic geometries.Definition 3.1. Let d â¥ 2. We define a Σ-generalized rotation in Rd as amatrix exponential eAt, such that t is the rotation angle and A â RdÃd isskew-symmetric with respect to Σ âª° 0, i.e. AΣ+ Î£Aâ¤ = 0.If Σ â» 0, then Σ-generalized rotations are equivalent to classical rotationsin a transformed space, such that lengths and angles are modified by scalingand shearing(see Lemma A.4 in appendix). Intuitively, ellipsoids, defined byEΣ = {x â Rd : xâ¤Î£â1x = 1}, are invariant under Σ-generalized rotations,since eAtEΣ = EΣ. We may similarly extend this definition for randomvariables.Definition 3.2. Let X be a d-dimensional r.v. with covariance Σ. X isauto-rotationally invariantâ¢ if there exists a nontrivial Σ-generalized rotation eAt with A Ì¸= 0 sucheAtX â¼ X ât â¥ 0, orâ¢ if there exists an orthogonal transformation P , such that at least onecomponent of Pâ¤X is mean 0 univariate Gaussian, i.e. (Pâ¤X)1 â¼113.1. SDE parameter identifiabilityN (0, Ï2).Conversely, X is not auto-rotationally invariant, if eAtX â¼ X ât â¥ 0 ad-mits only the solution A = 0, and none of its cross sections are univariateGaussian.Auto-rotational invariance is reminiscent of Maxwellâs theorem, which statesthat if a random vector X has independent components, and is invariant toall classical rotations, thenX â¼ N (0, Ï2Id) [67]. Indeed, all mean 0 isotropicGaussians are auto-rotationally invariant, since they are invariant under eAt,for any skew-symmetric A. To extend our definition to one dimensionalsubspaces, where rotations are undefined, we define mean 0 Gaussians asthe only univariate auto-rotationally invariant distributions by convention.However, auto-rotational invariance generalizes Maxwellâs characterizationin a few important ways. First, auto-rotational invariance generalizes rota-tions for anisotropic geometries defined by the covariance Σ, e.g. N(0,Σ) isauto-rotationally invariant (see Lemma A.5 in the appendix). Second, auto-rotational invariance only requires that there exists one matrix A such thateAtX â¼ X, rather than invariance to all rotations, e.g. a stack of GaussiansX = Unif(â1, 1) Ã N(0,Σ) would be auto-rotationally invariant. Further-more, auto-rotational invariance does not assume independent components,and therefore permits any rotationally symmetric distribution (after shear-ing and scaling Rd accordingly), e.g. uniform distributions over ellipsoidsand cylindrical distributions are auto-rotationally invariant. Finally, thematrix A is not required to be related to a rotation, even after shearing andscaling the space. Indeed, all degenerate r.v,âs are auto-rotationally invari-ant, since eAtX â¼ X would hold ifX is in the nullspace of A with probability1. It turns out that if X is non-degenerate and satisfies eAtX â¼ X for someA Ì¸= 0, thenX has a density whose level sets are ellipsoids (see Lemma A.6 inthe appendix). Thus, almost all random variables are not auto-rotationallyinvariant. Figure 3.1 illustrates examples of auto-rotationally invariant ran-dom variables, and Figure A.2 in the appendix shows that there are onlythree types of auto-rotational invariance: rank-degenerate random variables,random variables that are invariant to a non-degenerate Σ-generalized rota-tion, and random variables containing a Gaussian cross-section.3.1 SDE parameter identifiabilityWe are now ready for our main identifiability theorem, which characterizesidentifiability of linear additive noise SDEs based on the auto-rotational123.1. SDE parameter identifiability1.000.750.500.250.000.250.500.751.00x1.000.750.500.250.000.250.500.751.00y1.000.750.500.250.000.250.500.751.00z(a) X0 is spherical and invariant underall classical rotations.0.40.20.00.20.4x0.30.20.10.00.10.20.3y0.40.20.00.20.4z(b) X0 is cylindrical (after shearing andscaling transformations), and is invari-ant under rotation in the xy plane aftertransformation.4 3 2 1 0 1 2 3 4x1.000.750.500.250.000.250.500.751.00y(c) X0 is a Gaussian pancake, which ismean 0 univariate Gaussian in x.0.50 0.25 0.00 0.25 0.50x0.500.250.000.250.50y(d) X0 is degenerate and its support liesin the nullspace of A =[1 â10 0].Figure 3.1: Four examples of auto-rotationally invariant distributions areillustrated. By Theorem 3.3, these initial distributions lead to non-identifiability of linear additive noise SDEs from temporal marginals.133.1. SDE parameter identifiabilityinvariance of the initial distribution X0.Theorem 3.3. Let p0 be a probability distribution on Rd with finite mo-ments, and letÎ = {(A,H) : A â RdÃd, G â RdÃm, H = GGâ¤}denote the set of drift-diffusion pairs under a linear additive noise model(1.3). Then, given the initialization X0 â¼ p0, any parameter set (A,H) â Îis identifiable from its marginals {pt : t â¥ 0}A,H , in the sense thatif {pt : t â¥ 0}A,H = {pt : t â¥ 0}AË,HË then (A,H) = (AË, HË),if and only if p0 is not auto-rotationally invariant.In other words, identifiability of (A,H) from marginals is guaranteed if p0is not auto-rotationally invariant. However, if p0 is auto-rotationally in-variant, then there exist multiple drift-diffusion pairs (A,H) â Î that wouldobey the same marginals from p0.Proof. We first show that if X0 is auto-rotationally invariant, then there ex-ist multiple systems with equivalent marginals for each time, but distinct pa-rameter sets (A,H = GGâ¤) Ì¸= (AË, HË = GËGËâ¤). Since X0 is auto-rotationallyinvariant, then there exists a matrix A Ì¸= 0 such that eAtX0 â¼ X0, or thereexists a cross-section that is univariate Gaussian, i.e. (Pâ¤X0)1 â¼ N (0, Ï2)such that Pâ¤ = Pâ1.If there is A Ì¸= 0 such that eAtX0 â¼ X0, then it would be impossible todifferentiate between 0 drift and Ax, provided that A is skew-symmetricwith respect to the diffusion H = GGâ¤. Indeed, given that eAtX0 â¼ X0 andAH +HAâ¤ = 0, the SDEsdXt = GdWt X0 â¼ p0 (3.1)dXËt = AYt dt+GdWt Y0 â¼ p0 (3.2)would obey the same marginals. Intuitively, A induces a Σ-generalized ro-tation that is non-identifiable from X0, where Σ = Cov(X0). We may checkthat Xt â¼ Yt for all t â¥ 0 by comparing their solutions (A.7):Xt = X0 +â« t0GdWs = X0 +GWt â¼ X0 +N (0, Ht)Yt = eAtX0 +â« t0eA(tâs)GdWs â¼ X0 +N (0, Ht),143.1. SDE parameter identifiabilitywhere the fact thatâ« t0 eA(tâs)GdWs â¼ N (0, Ht) follows from Lemma A.9.Then, since the Ito integral is independent of X0, it follows that Xt â¼ Yt forall t â¥ 0.Alternatively, suppose that (Pâ¤X0)1 â¼ N (0, Ï2). Recall thatN (0, Ï2) is thestationary distribution of any one-dimensional OU process with drift:diffusivityratio â1 : 2Ï2 [13, 40]. Thus, the SDEsd(Pâ¤Xt) =[â12 00 A2,2]Pâ¤Xtdt+[Ï 00 G2,2]dWt X0 â¼ p0 (3.3)d(Pâ¤Yt) =[âÎ±22 00 A2,2]Pâ¤Ytdt+[Î±Ï 00 G2,2]dWt Y0 â¼ p0. (3.4)would obey the same marginals Pâ¤Xt â¼ Pâ¤Yt for all t â¥ 0. Both SDEs aredistributionally identical in the first component, since they would stay atthe stationary distribution N (0, Ï2). The SDEs also trivially obey the samelaw in the complementary subspace. It then follows that Xt â¼ Yt wouldhold ât â¥ 0 given the SDEsdXt = P[â12 00 A2,2]Pâ¤Xtdt+ P[Ï 00 G2,2]dWt X0 â¼ p0 (3.5)dYt = P[âÎ±22 00 A2,2]Pâ¤Ytdt+ P[Î±Ï 00 G2,2]dWt Y0 â¼ p0. (3.6)We now prove that identifiability is guaranteed if X0 is not auto-rotationallyinvariant. Suppose that Xt â¼ Yt for all t â¥ 0, such that the first system hasparameters (A,H = GGâ¤) and the second system has parameters (AË, HË =GËGËâ¤). Recall that linear additive noise SDEs admit closed form solutions:Xt = eAtX0 +â« t0eA(tâs)GdWsYt = eAËtX0 +â« t0eAË(tâs)GËdWs,such that eAtX0 â¥â« t0 eA(tâs)GdWs due to the independent increments ofBrownian motion. Furthermore,â« t0eA(tâs)GdWs â¼ N (0,Σt) such that Σt =â« t0eA(tâs)HeAâ¤(tâs)ds andâ« t0eAË(tâs)GËdWs â¼ N (0, Î£Ët) such that Î£Ët =â« t0eAË(tâs)HËeAËâ¤(tâs)ds.153.1. SDE parameter identifiabilitySinceXt â¼ Yt, then the cumulant generating functionsKXt(v) = log(E[evâ¤Xt ]) =KeAtX0(v)â 12vâ¤Î£tv and KYt(v) = log(E[evâ¤Yt ]) = KeAËtX0(v)â 12vâ¤Î£Ëtv areequal for all t â¥ 0. This gives us the equationKX0(eAâ¤tv)âKX0(eAËâ¤tv) =12vâ¤(Î£t â Î£Ët)v ât â¥ 0. (3.7)We may rewrite the cumulant generating functions in terms of their powerseries expansions [68]KX0(eAâ¤tv) =âân=01n!Tn(eAâ¤tv, . . . , eAâ¤tv) andKX0(eAËâ¤tv) =âân=01n!Tn(eAËâ¤tv, . . . , eAËâ¤tv),where the tensor Tn is the symmetric tensor of order n, which defines thenth cumulant of X0. To simplify the analysis, we perform the change ofvariables w(t) = eAâ¤tv, v = eâAâ¤tw(t). Let A¯ = AËâ¤ â Aâ¤, and supposefor contradiction that A¯ Ì¸= 0. Applying this change of variables to (3.7)yieldsâân=01n!(Tn(w(t), . . . , w(t))â Tn(eA¯tw(t), . . . , eA¯tw(t)))=12w(t)â¤(â« t0eâAsHeâAâ¤sdsâ eA¯â¤tâ« t0eâAËsHËeâAËâ¤sds eA¯t)w(t).Since the RHS is purely quadratic in w(t), it follows that for all n â¥ 3, wehaveTn(w(t), . . . , w(t)) = Tn(eA¯tw(t), . . . , eA¯tw(t)) (3.8)In other words, for every choice of n-repeated arguments (w, . . . , w), thevalue of Tn remains the same if we multiply each of the n components byeA¯t. Now, consider the Jordan form A¯ = P (D +N)Pâ1, such thateA¯t = PeDteNtPâ1Then we may apply (3.8) on the repeated argument (Pw, . . . , Pw)Tn(Pw, . . . , Pw) = Tn(PeDteNtw, . . . , PeDteNtw)163.1. SDE parameter identifiabilityHowever, we note that Sn(w, . . . , w) := Tn(Pw, . . . , Pw) represents the nthcumulant of the r.v. Pâ¤X0. Hence, ât â¥ 0, we have equality of the sym-metric tensorsSn(w, . . . , w) = Sn(eDteNtw, . . . , eDteNtw). (3.9)Now, we first suppose that the diagonal matrix D contains at least onenonzero eigenvalue λ1. Then, applying (3.9) on the argument (Ëi, . . . , iË),where iË is the first canonical unit vector, yieldsSn(Ëi, . . . , iË) = Sn(eNtiË, . . . , eNtiË)enÎ»1t= Sn(Ëi, . . . , iË)enÎ»1twhere we have used multilinearity, and the fact that for any upper triangularJordan block, N · iË = 0. For this equality to hold ât â¥ 0, we conclude thatSn(Ëi, . . . , iË) = 0, which implies that the nth cumulant of (Pâ¤X0)1 is zeroân â¥ 3. However, a classical result states that only Gaussian distributionshave finitely many nonzero cumulants [69]. Therefore, the first componentof Pâ¤X0 is a univariate mean-0 Gaussian if Σt Ì¸= Î£Ët. This contradicts thefact that X0 is not auto-rotationally invariant. We also note that if Σt = Î£Ët,then (3.7) implies that eA¯tX0 â¼ X0, since equality in all cumulants impliesequality in distribution. This also contradicts the non-auto-rotational in-variance of X0, unless A¯ = 0. Hence, we conclude that A¯ = 0 =â A = AË.This trivially yields eAtX0 â¼ eAËtX0, and forces the RHS of (3.7) to be 0,which in turn implies H = HË, from the definition of Σt and Î£Ët.A similar argument holds in the case where all eigenvalues of D are zero,such that eDteNt = eNt. Since N is an upper triangular Jordan block, wehave for some k â N smaller than d,eNt = I +Nt+(Nt)22+ · · Â·+ (Nt)kk!=kâj=0(Nt)jj!,In particular, all entries along the main diagonal of eNt are 1, i.e. (eNt)i,i =1, and all entries along the super diagonal are t, i.e. (eNt)i,i+1 = t. We maypick the vector v =[t 1 0 . . . 0]â¤and easily check thateNtv =kâj=0(Nt)jj!v =ï£®ï£¯ï£¯ï£¯ï£¯ï£¯ï£°t+ t10...0ï£¹ï£ºï£ºï£ºï£ºï£ºï£» = v + tËi173.2. Causal graph identification from SDE parametersTherefore, applying (3.9) and multilinearity yieldsSn(v, . . . , v) = Sn(v + tËi, . . . , v + tËi)=nâk=0(nk)tkSn(v, . . . , vï¸¸ ︷︷ ï¸¸nâk, iË, . . . , iËï¸¸ ︷︷ ︸k).Note that the k = 0 term yields Sn(v, . . . , v). We may rearrange and furtherapply multilinearity to yield:0 =nâk=1(nk)tkSn(v, . . . , vï¸¸ ︷︷ ï¸¸nâk, iË, . . . , iËï¸¸ ︷︷ ︸k)To satisfy this equation for all t, we must have Sn(jË, . . . , jËï¸¸ ︷︷ ï¸¸nâk, iË, . . . , iËï¸¸ ︷︷ ︸k) = 0for all k = 1, . . . , n. The case k = n implies that Sn(Ëi, . . . , iË) = 0. Asargued in the previous case D Ì¸= 0, this implies that (Pâ¤X0)1 â¼ N (0, Ï2),or Σt = Î£Ët. This yields identifiability A = AË and H = HË provided that X0is not auto-rotationally invariant.3.2 Causal graph identification from SDEparametersGiven that the drift and diffusion are identified, we can gain important in-sights into the systemâs causal structure. It has already been shown thatknowing drift and diffusion provides the systemâs post-intervention distribu-tions [5, 70]. We prove under basic conditions that the causal graph G of anadditive noise SDE can be recovered from the drift b(Xt) and observationaldiffusion H = GGâ¤.Consider a d-dimensional additive noise SDE (1.2). By isolating the jthvariable and integrating, we obtain:X(j)t = X(j)0 +â« t0bj(Xα(j)s ) ds+mâk=1Gj,kW(k)t , (3.10)where α(j) â [d] is the set of variables that are used as input in bj(Xt).Then, (3.10) is an example of a dynamic structural causal model [66]. In183.2. Causal graph identification from SDE parametersparticular, we can define a causal graph G over vertices V = [d] and edgesE = {i â j : j â [d], i â α(j)} [66][Definition 2]. Intuitively, given j â V ,we include the edge i â j if X(i)t influences the evolution of X(j)t throughthe drift. Under an additive noise model, the diffusion Ï(Xt) = G is space-independent, and does not contribute edges between vertices. However, ifG is anisotropic, then the diffusion induces correlations between variablesthrough shared driving noise [37]. In this paper, we interpret shared drivingnoise as evidence of a latent confounder.Definition 3.4. Let X evolve according to a d-dimensional additive noiseSDE: dXt = b(Xt) dt + G dWt. We define a multi-directed edge eË =(i1, . . . ip) â G if variables X(i1), . . . , X(ip) depend on a common noise sourceW(k)t , k â [m].This definition follows standard convention in the static setting, where la-tent confounders introduce dependencies between noise terms in a functionalcausal model. In this framework, multidirected edges represent sets of vari-ables whose noise terms are dependent due to a shared latent confounder[71]. In the dynamic setting, shared driving noise has similarly been usedto to model latent confounders (see [72] and [66][Example 4]). However,we note that this is not a comprehensive characterization of latent con-founders in the dynamic setting, since confounding effects cannot always beattributed to diffusion. To address this, [36] uses blunt edges to specificallydenote correlated noise between two variables.Under this causal framework, we show that the drift and diffusion identifyall directed edges in the underlying causal graph, as well as all pairwisemultidirected edges, i.e. bidirected edges.Lemma 3.5. Let Xt evolve according to a d-dimensional additive noiseSDE: dXt = b(Xt)dt+GdWt, and let G = (V = [d], E, EË) be its causal graph,over the variables X(1), · · · , X(d), with directed edge set E and multidirectededge set EË.a. There exists a directed edge e = i â j in G if and only if bj(Xt)depends on X(i)t .b. There exists a multidirected edge eË containing (i, j) in G, if and onlyif Hi,j = (GGâ¤)i,j Ì¸= 0.Thus, G is identified by b(Xt) and H = GGâ¤ if no more than two variablesshare a common noise source.193.2. Causal graph identification from SDE parametersProof. The first claim follows immediately from the definition of the dynamicstructural causal model (3.10).To prove the second claim, it is easy to see from (3.10) that X(i1), . . . , X(ip)depends on a common noise source if and only if Gi1,k, . . . Gip,k Ì¸= 0 forsome column k â [m] in the diffusion G. Indeed, the existence of such acolumn implies that X(i1), . . . , X(ip) each depend on W(k)t , and the otherdirection is given by the fact that Wt is a multivariate Gaussian process,with independent components. If two variables share a noise source, thiscondition is equivalent to Hi,j = Gi · Gj Ì¸= 0. We therefore conclude thatthere is a multidirected edge containing (i, j) in G if and only if Hi,j Ì¸= 0.Lemma 3.5 enables arbitrary causal structures to be learned from the drift,as well as pairwise shared noise structures from the diffusion. CombiningTheorem 3.3 and Lemma 3.5 shows that we can identify the full causalstructure of a linear additive noise SDE from its marginals if and only if X0is non-auto-rotationally invariant. This result is particularly important forGRN inference, since our mathematical framework encompasses rotationaldrift, which can model feedback and repressilator dynamics [33].20Chapter 4Our parameter estimationmethodIn this section, we introduce our method, Alternating Projection ParameterEstimation from X0 (APPEX), which estimates a processâ drift, diffusion,and causal graph from observed temporal marginals. Figure 4.1 illustratesthe general strategy of the APPEX algorithm. We focus on linear additivenoise SDEs, since Theorem 3.3 guarantees identifiability given that X0 is notauto-rotationally invariant. However, our algorithm can be used to estimateany additive noise SDE in the form of Equation 1.2. In this section, wegive an intuitive explanation at the population level, using the processâ truetemporal marginals pi â¼ Xti . In Chapter 5, we perform experiments usingempirical marginals pËi, with finite samples per marginal.Let p0, . . . , pNâ1 denote a set of observed temporal marginals. Let (A(0), H(0))denote initial guesses for the drift and diffusion. The idea is to use an al-ternating optimization algorithm to obtain increasingly better estimates forthe drift and diffusion matrices. Formally, we consider the spacesD = {q : q0 = p0, . . . , qNâ1 = pNâ1} (4.1)A ={p : â (A,H) â RdÃd Ã Sdâª°0 s.tp(x, t+ât | y, t) = N(eAâty,â« ât0eAsH eAâ¤s ds)}(4.2)D is the set of laws on paths, which coincide with the N temporal marginalsat their respective times. A is the set of laws given by linear additivenoise SDEs. To find the underlying linear additive noise SDE in the in-tersection A â© D, our method alternates between information projections,argminqâDDKL(qâ¥p) and moment projections, argminpâADKL(qâ¥p). Theinformation projection step corresponds to trajectory inference, whose goal21Chapter 4. Our parameter estimation methodFigure 4.1: Visualization of our parameter estimation algorithm APPEX.Given observed temporal marginals from an underlying SDE (left), APPEXalternates between trajectory inference (middle) and MLE parameter esti-mation (right) in order to find the SDE parameters that best represent thetemporal marginal observations.is to find the law on paths q â D, which best aligns with the reference mea-sure p. The moment projection then corresponds to maximum likelihoodparameter estimation, since it amounts to determining the law on pathsp â A over linear additive noise SDEs, which best aligns with the law of theestimated trajectories q. Indeed, in the infinite-data limit, MLE is equiv-alent to minimizing relative entropy, i.e. KL divergence [73, 74]. At thepopulation level, an iteration entails the updates:q(k) = argminqâDDKL(qâ¥pA(kâ1),H(kâ1)) (4.3)A(k) = argminAâRdÃdDKL(q(k)â¥pA,H(kâ1)) (4.4)H(k) = argminHâSdâª°0DKL(q(k)â¥pA(k),H). (4.5)At iteration k, we start by using the SDE parameters estimated from theprevious iteration to define a reference SDE pA(kâ1),H(kâ1) â A. Using thisreference SDE and the set of observed temporal marginals, we perform tra-jectory inference (4.3), to estimate the path distribution q(k) â D that min-imizes relative entropy to the reference SDE, while satisfying the marginalconstraints (4.1) (see the middle plot of Figure 4.1). We then perform atwo-step maximum likelihood estimation, in order to estimate the drift anddiffusion from the estimated trajectories. This in turn updates the referenceSDE pA(k),H(k) â A for the next iteration. To perform APPEX on a general22Chapter 4. Our parameter estimation methodadditive noise SDE (1.2), the spaces A and the drift maximum likelihoodparameter estimation step (4.4) would be modified accordingly. We show inCorollary A.18 of the appendix that there is a closed form diffusion MLEfor general additive noise SDEs, once the drift has been estimated.Previous works [8, 22, 23] implement similar iterative schemes, which al-ternate between trajectory inference and MLE parameter estimation. Incontrast to our method, each of these works assumes that the diffusion ofthe process is known, and hence invariant across iterations. One reason forthis assumption is to ensure finite KL divergence between the estimatedtrajectories and estimated SDEs for all iterations [22]. Given continuouslyobserved marginals of a d-dimensional process, the KL divergence betweentwo laws on paths q, p is taken over the path space â¦ = C([0, T ],Rd), suchthatDKL(qâ¥p) =â«â¦log(dqdp(Ï))dq(Ï).Thus, DKL(qâ¥p) will only be finite between two laws on paths if we candefine the Radon-Nikodym derivative dqdp over path space â¦. By Girsanovâstheorem, this is only ensured when both processes share the same diffu-sion [22, 23]. However, if we only consider measurements from a finitenumber of observed marginals, the KL divergence of the discretized pro-cesses qN , pN over the path space projected on the N measurement times,ΣD = C({ti}Nâ1i=0 ,Rd), can be decomposed as follows, see [75]:DKL(qNâ¥pN ) =Nâ2âi=0DKL(qi,i+1â¥pi,i+1)âNâ2âi=1DKL(qiâ¥pi), (4.6)where pi,i+1, qi,i+1 are the joint probability measures of q, p, restricted totimes ti and ti+1, and similarly qi, pi are the marginals of q, p at time ti. Wenote that DKL(qiâ¥pi) â¤ DKL(qi,i+1â¥pi,i+1) holds by the data-processing in-equality. Hence, DKL(qNâ¥pN ) will be finite as long as DKL(qi,i+1â¥pi,i+1) =â«RdÃRd log(dqi,i+1dpi,i+1(x))qi,i+1(x) dx <â for each i = 0, . . . , Nâ1. In particu-lar, if p and q are the laws of two different drift-diffusion SDEs with distinctnon-degenerate diffusions, then the KL divergence DKL(qNâ¥pN ) over thediscretized path space ΣD will be finite. Since we consider a setting with afinite number of temporal marginal observations in this work, we may con-sider laws on paths discretized over these observations. This allows us toboth consider diffusions which are not known in advance, and to improveour diffusion estimates with respect to KL divergence after each iteration.234.1. Trajectory inference via anisotropic entropy-regularized optimal transport4.1 Trajectory inference via anisotropicentropy-regularized optimal transportTo begin each iteration of APPEX, we infer trajectories given the set oftemporal marginals {pi}Nâ1i=0 and a reference SDE pA(kâ1),Hkâ1) â A, bysolving an associated Schro¨dinger Bridge problem. Under certain conditions,the Schro¨dinger Bridge problem (1.5) can be solved using the variationalentropic optimal transport (EOT) formulationÏâ = argminÏâÎ (µ,Î½)â«c(x, y) dÏ(x, y) + Ïµ2DKL(Ïâ¥µâ ν). (4.7)Intuitively, Ïâ transports the probability measure µ to ν while minimizinga total cost, based on the function c(x, y) and the entropic regularizationϵ2 > 0. Objectives (1.5) and (4.7) are equivalent [8, 29], under the referencemeasure dK(x, y) = K(x, y)dµ(x)dν(y) with density K(x, y) â e c(x,y)ϵ2 . Theconvention in the literature is to set c(x, y) = â¥yâxâ¥22 , to define an isotropicGaussian law K [8, 13, 15, 76]. The EOT problem (4.7) can then be solvedvia Sinkhornâs algorithm.However, the standard EOT problem (4.7) only considers a scalar regu-larization parameter ϵ2, which means that it can only directly model ref-erence SDEs with isotropic diffusion (see Remark A.17 in the appendix).To generalize trajectory inference for SDEs with anisotropic diffusion, weformalize âanisotropic entropy-regularized optimal transportâ (AEOT), byconsidering transition kernels K with custom mean and covariance, e.g.parameterizing K with A and H = GGâ¤ rather than the scalar ϵ2 > 0. In-deed, we can consider more general Gaussian transition kernels, Kθ(x, y) =exp( (yâµ(x,Î¸))â¤Î£(Î¸)â1(yâµ(x,Î¸))2 ). This allows us to model a transition y|x â¼N (µ(x, θ),Σ(θ)), whose covariance Σ(θ) is possibly anisotropic. Under thisformulation, the cost is given by the inner productcθ(x, y) =(y â µ(x, Î¸))â¤Î£(Î¸)â1(y â µ(x, θ))2.Trajectory inference would amount to solving the AEOT problem:Ïâ = argminÏâÎ (µ,Î½)â«(y â µ(x, Î¸))â¤Î£(Î¸)â1(y â µ(x, Î¸))2dÏ(x, y) +DKL(Ïâ¥µâ ν),(4.8)244.2. Parameter estimation via MLEwhere the entropic regularization is captured in the matrix Σ(θ) and functionµ(x, θ). The AEOT problem over a pair of marginals (4.8) can be solved viaSinkhornâs algorithm with inputs µ, ν,Kθ. Proposition A.16 in the Appendixthen shows that the solution to the multi-marginal trajectory inference step(4.3) is the joint distribution given by the Markov concatenation of couplings(see Definition A.15 in the Appendix)Ïâ = Ï0,1 â¦ . . . â¦ ÏNâ2,Nâ1,where Ïi,i+1 is the AEOT solution (4.8) with marginals µ = pi, ν = pi+1 anda transition kernel determined by the reference SDE.Remark 4.1 (Application to empirical marginals). Proposition A.16 is for-malized for exact temporal marginals pi. However, in practice, we observeempirical measures pËi, which converge in distribution to pi, in the limit ofinfinite data. It is for these empirical measures pËi which we actually computethe estimated couplings ÏËi,i+1 in practice. This corresponds to numericallyestimating the AEOT solutionÏËi,i+1 = arg minÏâÎ (pËi,pËi+1)â«(y â µ(x, Î¸))â¤Î£(Î¸)â1(y â µ(x, Î¸))2dÏ(x, y)+DKL(Ï||pËi â pËi+1).(4.9)By the main theorem in [77], it holds that ÏËi,i+1 converges in distributionalmost surely to the minimizer of (4.8) with µ = pi and ν = pi+1 as thenumber of samples goes to infinity for each time i and i+1. Combining thiswith the previous proposition, we see that estimating (??) for each pair oftimes i and i + 1 is asymptotically equivalent to the KL minimization stepin (4.3).We also note that (??) is practically solved using Sinkhornâs algorithm withempirical marginals pËi, pËi+1 and the discretized transition kernel Kθ of thereference SDE. For example, if we have M samples per empirical marginal,then pËi â¼ Unif(x(j)ti : j = 1, . . .M) and pËi+1 â¼ Unif(x(j)ti+1: j = 1, . . .M)would both be discrete uniform distributions over their samples, and Kθ âRMÃM would be a square matrix, such that entry Kθj,k is obtained by apply-ing the given transition kernel on the data points x(j)tiand x(k)ti+1.4.2 Parameter estimation via MLETo optimize objectives (4.4) and (4.5) for each iteration of APPEX, we re-quire maximum likelihood estimators for the SDE parameters, given multiple254.2. Parameter estimation via MLEobserved trajectories from [0, T ]. In the context of iteration k of APPEX,these are the trajectories sampled from the law on paths q(k) obtained fromthe trajectory inference step. We derive closed-form maximum likelihoodestimators for the linear additive noise SDE from Equation 1.3, given mul-tiple trajectories in Proposition 4.2. Closed form drift MLE solutions aregenerally unavailable for additive noise SDEs, however we derive a closedform diffusion estimate (which depends on the drift) for general additivenoise SDEs in Corollary A.18 of the appendix.Proposition 4.2 (MLE estimators for drift and diffusion of SDE (1.3) frommultiple trajectories). Given M trajectories over N equally spaced times:{X(j)iât : i â 0, ..., N â 1, j â 0, ...,M â 1} sampled from the linear addi-tive noise SDE (1.3), the maximum likelihood solution for linear drift isapproximated byAË =1âtï£«ï£­Nâ2âi=0Mâ1âj=0(âX(j)i )X(j)â¤iï£¶ï£¸ï£«ï£­Nâ2âi=0Mâ1âj=0X(j)i X(j)â¤iï£¶ï£¸â1 (4.10)Mââââ 1ât(Nâ2âi=0Epi,i+1 [(âXi)Xâ¤i ])(Nâ1âi=0Epi [XiXâ¤i ])â1, (4.11)where pi,i+1 is the joint measure over Xti and Xti+1. Similarly, the maximumlikelihood solution for diffusion is approximated byHË =1MTNâ2âi=0Mâ1âj=0(âX(j)i âAX(j)i ât)(âX(j)i âAX(j)i ât)â¤)(4.12)Mââââ 1TNâ2âi=0Epi,i+1 [(âXi âAXiât) (âXi âAXiât)â¤] (4.13)Proof. See Section A.8 of the Appendix. We note that estimators AË andHË were derived using the discretized transition kernel, Xi+1|Xi â¼ N (Xi +AXiât,Hât). We derive the maximum likelihood estimators with the exacttransition kernel Xi+1|Xi â¼ N (eAâtXi, Hât) in the one dimensional casein Section A.8.We note that the MLE drift estimator for A does not depend on the diffusionH, but the MLE estimator for H depends on A. We therefore estimate driftfirst in each iteration of APPEX.264.3. The APPEX algorithm4.3 The APPEX algorithmWe summarize APPEXâs implementation in Algorithm 1. In each iteration,trajectory inference (4.3) is performed by solving the Anisotropic EntropyRegularized Optimal Transport (AEOT) problem (4.8), and the drift anddiffusion estimates are updated via MLE. By default, we use Sinkhornâs algo-rithm to solve AEOT and we use closed form MLE solutions when available.Algorithm 1 Parameter estimation for an additive noise SDE from tem-poral marginals with APPEX1: Input: Observed marginals pËi, i = 0, . . . , N â 1, number of iterationsK, ât2: Result: Estimated drift function bË and additive noise HË3: Initialize: bËâ 0, H â Ï2Id, k â 04: while k < K do5: for i = 1, . . . , N do6: ÏËi,i+1 â AEOT(bË, HË, pËiâ1, pËi,ât)7: end for8: Sample-Trajectoriesâ ÏËNâ1,N â¦ . . . â¦ ÏË1,2(p0)9: bËâ MLEfit(Sample-Trajectories)10: HË â MLEfit(Sample-Trajectories, bË)11: k â k + 112: end while13: G â Estimate-Causal-Graph(bË, HË, ϵ)Remark 4.3 (Application to causal discovery). As proven in Lemma 3.5and visualized in Figure 2.1, the causal graph can be derived from the SDEparameters. For nonlinear drift b(Xt), we include the edge i â j if thejth component of the drift, bj, is a function of X(i). For linear drift onemay simply examine Aj,i. For multi-edges, one can check Hi,j to determinewhether there is an unobserved confounder causing X(i) and X(j). We there-fore include a causal discovery step in line 13 of Algorithm 1. We acknowl-edge that many causal algorithms and independence tests are applicable, andleave the implementation as general as possible.Because the trajectory inference and parameter estimation subproceduresare optimal with respect to minimizing KL divergence, we can show thatAPPEXâs estimates reduce relative entropy to the true solution.Lemma 4.4. Suppose that A and H are the true drift and diffusion parame-274.3. The APPEX algorithmters of a linear additive noise SDE with temporal marginals p0, ..., pNâ1. Letq(k) â D be the law of the estimated trajectories at iteration k of APPEX,and pAk,Hk â A be the law of the law of the estimated linear additive noiseSDE, such that A0 â RdÃd and H0 â Sdâª°0. Then, the relative entropy be-tween the estimated law of the trajectories and the law of the estimated linearadditive noise SDE is decreasing with each iteration [21]:DKL(q(k+1)â¥pAk+1,Hk+1) â¤ DKL(q(k)â¥pAk,Hk) âk â¥ 0.Furthermore, if p0 is not auto-rotationally invariant, then as the number ofobserved marginals N â â, the relative entropy between the estimated lawof the trajectories and the law of the estimated SDE is uniquely minimizedby the underlying SDE pA,H :infqâD,pâADKL(qâ¥p) = 0 ââ q = p = pA,HProof. Proposition A.16 shows that q(k+1) â D minimizes relative entropyto pAk,Hk . Similarly, any MLE solution pA(k+1),H(k+1) minimizes relativeentropy to q(k+1). In particular, Proposition 4.2 approximates the MLEparameters Ak+1 â RdÃd and H(k+1) â Sdâª°0. By construction of the iterativescheme, we have that By construction of the iterative scheme, we have thatDKL(q(k+1)â¥pA(k+1),H(k+1)) â¤ DKL(q(k+1)â¥pA(k+1),H(k))â¤ DKL(q(k+1)â¥pA(k),H(k))â¤ DKL(q(k)â¥pA(k),H(k)), (4.14)which proves the first claim of the lemma. Then, by Theorem 3.3, if X0 isnot auto-rotationally invariant, it follows that pA,H is the law of the uniquelinear additive noise SDE, which obeys the marginal constraints p0, . . . pNâ1as N ââ. Equivalently, we have A â©D = {pA,H} andinfqâD,pâADKL(qâ¥p) = 0 ââ q = p = pA,H ,Since DKL(q(k)â¥pA(k),H(k)) is non-increasing for successive iterations, thisimplies that the drift and diffusion estimates A(k), H(k) progressively fit thedata with each iteration. In this sense, our method approaches the truesolution, which is unique, if the data is appropriately initialized and hassufficient observations.28Chapter 5ExperimentsIn a wide range of experiments on simulated data from linear additive noiseSDEs, we find that APPEX accurately identifies the underlying SDE fromits temporal marginals and that it outperforms the popular Waddington-OT method [12]. In Section 5.1, we assess the validity of Theorem 3.3 byshowing that APPEX can learn the drift and diffusion parameters of clas-sical non-identifiable SDEs, if X0 is not auto-rotationally invariant. Wethen demonstrate APPEXâs effectiveness on randomly generated higher-dimensional SDEs in Section 5.2. Finally, in Section 5.3, we show thatAPPEX can use its parameter estimates to recover the underlying causalgraph, thus supporting Proposition 3.5. An empirical experiment for con-sistency is also provided in Appendix A.9.1, which shows that as the numberof samples per time marginal increases, APPEXâs drift and diffusion esti-mates converge to the true parameters.For each experiment, we simulate data for linear additive noise SDEs us-ing Euler-Maruyama discretization with dtEM = 0.01, such that we gen-erate M = 500 trajectories, each observed across 100 time steps. Eachd-dimensional trajectory is initialized with p0 â¼ Unif{xi}di=1 where {xi}di=1are random linearly independent vectors with each entry having a magni-tude uniformly sampled between 2 and 10. It follows from Proposition A.8and Theorem 3.3, that this initial distribution ensures that the drift anddiffusion are theoretically identifiable. To model the âmarginals-onlyâ set-ting, we subsample at the rate dt = 0.05 to produce N = 20 marginals withM = 500 observations per time.To perform parameter estimation, we use 30 iterations of our APPEX al-gorithm, such that the first reference SDE is an isotropic Brownian motiondXt = ÏdWt, even when the true process has anisotropic noise. To modelthe realistic setting where diffusion is not precisely known, we assume thatthe initial guess of the diffusionâs trace is within an order of magnitude ofthe ground truth trace tr(H). Specifically, we initialize A(0) = 0, H(0) =295.1. Revisiting non-identifiable SDEsÏ2Id s.t. Ï2 â¼ tr(H)10Unif(â1,1). Because the first reference SDE is ÏdWt,we note that the first iteration of APPEX is equivalent to performing MLEon trajectories that were inferred using the Waddington-OT (WOT) methodfrom [12], based on standard entropy regularized OT (4.7) with ϵ2 = Ï2ât.APPEX and WOT are distinguished by the fact that APPEX benefits fromfurther iterations, and allows the reference SDE to have anisotropic diffusionand non-zero drift. For each iteration of APPEX, we use Sinkhornâs algo-rithm to solve the anisotropic entropy-regularized optimal transport prob-lem for trajectory inference, and we use the closed form MLE solutions fromProposition 4.2. Due to time complexity and numerical stability, we uselinearized discretizations for the Gaussian transition kernels (unless d = 1).5.1 Revisiting non-identifiable SDEsWe consider three non-identifiable pairs of SDEs from the literature (summa-rized in Appendix A.1), but revise X0 to be non-auto-rotationally invariant.By Theorem 3.3, parameter estimation is feasible under the revised initialdistribution. 10 replicates of each experiment were performed, such that thereplicates featured data from different valid p0 initializations and differentinitial diffusivities Ï2. We track the mean absolute error (MAE), plotted inFigure 5.1, between the true drift/diffusion parameters and their estimates.The results are consistent with Theorem 3.3 and Lemma 4.4. In each ex-ample, APPEX is able to estimate the drift and diffusion by iterativelyimproving upon both estimates, as evidenced by decreasing MAEs as thenumber of iterations increases. In particular, we observe the worst perfor-mance from the first iteration, which is equivalent to performing parameterestimation following trajectory inference by WOT. The most difficult SDEto infer was dXt = â10Xtdt+â10dWt, due to the high noise scale. We alsonote that APPEX succeeded in estimating the degenerate diffusion from theSDEs in Example A.3, despite the fact that degenerate diffusion can resultin infinite KL divergence with respect to a misspecified reference SDE.5.2 Higher dimensional randomly generatedlinear additive noise SDEsWe now consider parameter estimation on a broad range of higher dimen-sional random linear additive noise SDEs. For each dimension d = 3, . . . , 10,305.2. Higher dimensional randomly generated linear additive noise SDEs0 5 10 15 20 25 30Iteration0.00.10.20.30.40.5MAEMAE between estimated A and true AMAE between estimated H and true H(a) dXt = âXt dt+ dWt from Ex. A.10 5 10 15 20 25 30Iteration123456MAEMAE between estimated A and true AMAE between estimated H and true H(b) dXt = â10Xt dt+â10dWt from Ex.A.10 5 10 15 20 25 30Iteration0.00.20.40.60.81.01.2MAEMAE between estimated A and true AMAE between estimated H and true H(c) dXt =[1 00 1]dWt from Ex. A.20 5 10 15 20 25 30Iteration0.00.20.40.60.81.01.2MAEMAE between estimated A and true AMAE between estimated H and true H(d) dXt =[0 1â1 â0]Xtdt+[1 00 1]dWtfrom Ex. A.20 5 10 15 20 25 30Iteration012345MAEMAE between estimated A and true AMAE between estimated H and true H(e) dXt =[1 21 0]Xtdt+[1 2â1 â2]dWtfrom Ex.A.30 5 10 15 20 25 30Iteration012345MAEMAE between estimated A and true AMAE between estimated H and true H(f) dXt=[134323 â 13]Xtdt+[1 2â1 â2]dWtfrom Ex. A.3Figure 5.1: The mean absolute error for estimates of A and H using APPEXis shown per iteration for all three pairs of SDEs from Section A.1we generate 10 random SDEs. To create each drift matrix A, we randomlysample each of its d2 entries from Unif(â5, 5), and verify that its principal315.2. Higher dimensional randomly generated linear additive noise SDEseigenvalue is less than 1. This allows us to consider non-stationary pro-cesses, while obeying reasonable growth conditions. To create each diffusionG, we randomly initialize each of its d2 entries via Unif(â1, 1) and thenset H = GGâ¤. For numerical stability, we perform Sinkhorn on the loga-rithmic scale for this experiment. We evaluate each estimateâs performanceusing MAE and correlation with the true parameters, where correlation iscomputed by vectorizing the corresponding matrices.d A Estimation GGâ¤ EstimationWOT APPEX WOT APPEX3 0.351± 0.04 0.237± 0.04 0.793± 0.205 0.147± 0.0304 0.730± 0.067 0.328± 0.041 1.549± 0.439 0.415± 0.0705 0.912± 0.060 0.602± 0.195 2.174± 0.702 0.362± 0.0396 1.43± 0.170 0.358± 0.020 18.010± 8.636 0.256± 0.0467 1.480± 0.132 0.360± 0.015 6.807± 1.724 0.345± 0.0378 1.862± 0.137 0.460± 0.015 5.472± 1.266 0.359± 0.0199 1.803± 0.222 0.487± 0.016 8.134± 3.024 0.454± 0.12210 1.670± 0.241 0.439± 0.019 35.122± 28.529 0.317± 0.025Table 5.1: Mean absolute error (MAE) of estimated drift and diffusion fordimensions 3-10.d A Estimation GGâ¤ EstimationWOT APPEX WOT APPEX3 0.996± 0.001 0.998± 0.001 0.837± 0.048 0.985± 0.0054 0.943± 0.015 0.987± 0.005 0.729± 0.039 0.865± 0.0315 0.921± 0.016 0.952± 0.030 0.728± 0.040 0.909± 0.0186 0.794± 0.040 0.986± 0.001 0.530± 0.056 0.961± 0.0077 0.792± 0.029 0.988± 0.001 0.595± 0.037 0.946± 0.0128 0.699± 0.035 0.981± 0.002 0.611± 0.042 0.949± 0.0069 0.740± 0.033 0.978± 0.002 0.615± 0.025 0.919± 0.03310 0.760± 0.049 0.983± 0.001 0.641± 0.041 0.960± 0.006Table 5.2: Correlation between estimated and true drift and diffusion fordimensions 3-10.The results demonstrate that APPEX estimates both SDE parameter setsrobustly across all settings. Importantly, APPEX can handle arbitrary ad-ditive noise structures H = GGâ¤, as evidenced by the fact that its drift anddiffusion estimates have low MAE and high correlation with respect to the325.3. Causal discovery0 5 10 15 20 25 30Iteration0.800.850.900.951.00CorrelationCorrelation between estimated A and true ACorrelation between estimated H and true H(a) d = 30 5 10 15 20 25 30Iteration0.600.650.700.750.800.850.900.951.00CorrelationCorrelation between estimated A and true ACorrelation between estimated H and true H(b) d = 10Figure 5.2: The correlation between the estimated and true SDE parametersis plotted per iteration across 10 random linear additive noise SDEs fordimensions 3 and 10true parameters. This is significant because previous literature has focusedon the setting of isotropic noise, but in practice, we expect correlated noisestructures, as well as unequal noise along the main diagonal. In contrast,although WOT estimates the drift somewhat decently, it is unable to esti-mate the diffusion accurately, particularly in higher dimensions, since it isconstrained to isotropic noise in its reference SDE. Figure 5.2 shows howAPPEX is able to re-orient incorrect diffusion priors to closely match thetrue diffusion of the underlying SDE.5.3 Causal discoveryWe conclude with experiments that test our methodâs ability to recover theunderlying causal graph. We first consider causal discovery from systemswithout latent confounders, and then consider causal discovery from systemswith latent pairwise confounders through correlated diffusion.In both settings, we generate random SDEs of dimension d = 3, 5, 10 whosecausal structure is determined by Erdo¨s-Renyi graphs G(d, p), such that eachof the d(dâ1) possible simple directed edges are included with probability p.As in Section 5.2, we ensure that the maximal eigenvalue of the randomlygenerated ground truth drift matrix is at most 1 to prevent unboundedgrowth. Proposition 3.5 states that simple directed edges e = iâ j â E inthe causal graph G = ([d], E, EË) are characterized by the condition Aj,i Ì¸= 0.To simulate ground truth edges, we follow standard convention, by simu-lating edge weights Aj,i uniformly via Unif(â5,â0.5)âªUnif(0.5, 5) [78, 79].335.3. Causal discoveryWe infer the presence of an edge i â j, if in the estimated drift, AËj,i > ϵ(corresponding to a positive edge weight) or AËj,i < âÏµ (corresponding to anegative edge weight). We choose our threshold ϵ = 0.5 to match the mini-mal edge weight magnitude from our simulations. Inference of simple edgesis then measured according to the structural Hamming distance between Gand the estimated graph Gest, which considers the classification of positivedirected edges, negative directed edges, and absent edges. The structuralHamming distance adds 1 for every misclassified simple edge:dsimple(G,G) =â(i,j)â[d]Ã[d]1{sgn(Aj,i) Ì¸= sgn(AËj,i)1|AËj,i|>Ïµ}. (5.1)If the observed system does not have any latent confounders from correlateddiffusion, we enforce this by setting the diffusion matrix G to be strictlydiagonal. This ensures that the only edges in G are simple edges from thedrift. As in the experiments on random matrices, we set each of the ddiagonal entries of G via Unif(â1, 1). Conversely, when modeling pairwiselatent confounders, we introduce correlated diffusion by selecting a randomsubset of columns in G to have precisely two nonzero entries. In particular,we initialize all entries along the diagonal of G to be 1, and for each of theselected columns, e.g. column i, we set one non-diagonal entry Gi,j = 1,to model the pairwise latent confounder Ui,j . All other entries are set to0. The number of selected columns with pairwise latent confounders issampled from Unif{1, · · · â2d/3â}, e.g. we sample from Unif{1, 2} for d = 3,Unif{1, 2, 3} for d = 5, and Unif{1, . . . 6} for d = 10. Proposition 3.5 statesthat bidirected edges (i, j) â EË from pairwise latent confounders Ui,j , arecharacterized by the condition Hi,j = Gi · Gj Ì¸= 0. We infer the presenceof Ui,j , if in the estimated diffusion, |HËi,j | > ϵ. We set ϵ = 1 to matchthe data construction (i, j) â EË ââ Hi,j = Gi · Gj â¥ 1. Inference ofbidirected edges is measured according to the structural Hamming distance,which adds 1 for every misclassified bidirected edge:dlatent(G,Gest) =â(i,j)â[d]Ã[d]:i Ì¸=j1{(|Hi,j | > ϵ â© |HËi,j | â¤ Ïµ)âª (|Hi,j | â¤ ϵ â© |HËi,j | > ϵ)}.(5.2)Table 5.3 reports the mean structural Hamming distances for simple edgesfrom causal graphs estimated by WOT and APPEX, in the setting withoutlatent confounders. APPEX outperforms WOT across all tested dimensionsd â {3, 5, 10} and edge probabilities p â {0.1, 0.25, 0.5}, with the discrepancy345.3. Causal discoverybecoming more pronounced as the causal structure grows in complexity athigher dimensions and edge densities. We observe that APPEX can consis-tently infer complex causal relationships, whereas WOTâs inference is oftenlimited, due to mis-specified diffusion. For example, the first row in Fig-ure 5.3 demonstrates that APPEX recovers the systemâs negative feedbackstructure, while WOT introduces additional cycles in order to compensatefor a poorly specified reference diffusion.Table 5.4 reports the mean structural Hamming distances for simple andbidirected edges (from latent pairwise confounders) from causal graphs es-timated by WOT and APPEX, in the setting with correlated diffusion.The probability p of connecting a simple directed edge i â j was fixedat p = 0.25. As in the setting with uncorrelated diffusion, APPEX out-performs WOT in all evaluated settings, particularly for higher-dimensionalsystems.Table 5.3: Average Structural Hamming Distance (lower is better) withvarying dimensions and random edge probabilities p.d p = 0.1 p = 0.25 p = 0.5WOT APPEX WOT APPEX WOT APPEX3 0.40± 0.40 0.00± 0.00 1.80± 0.76 0.30± 0.30 1.40± 0.65 0.00± 0.005 3.80± 2.13 0.10± 0.10 3.56± 1.70 0.00± 0.00 4.40± 1.14 0.60± 0.2210 38.30± 8.53 0.50± 0.40 48.60± 6.55 0.80± 0.39 37.4± 4.66 3.30± 0.70Table 5.4: Average Structural Hamming Distance (HD) for simple edges,and latent confounders (lower is better) with varying dimensions and randomedge probability p = 0.25.d HD for simple edges HD for latent confoundersWOT APPEX WOT APPEX3 1.20± 0.51 0.7± 0.40 1.00± 0.30 0.00± 0.005 4.80± 1.47 1.60± 0.87 3.50± 0.64 0.10± 0.1010 38.00± 6.30 3.00± 0.77 23.80± 3.86 2.90± 1.39355.3. Causal discoveryFigure 5.3: The true and estimated causal graphs by WOT and APPEX forfour random SDEs of dimensions are illustrated. Red edges i â j repre-sent negative edge weights, such that X(i) negatively regulates X(j). Greenedges iâ j represent positive edge weights, such that X(i) positively regu-lates X(j). Dashed edges represent the effect of a latent confounder Ui,j onobserved variables i, j through correlated diffusion.36Chapter 6DiscussionMany recent works have focused on inferring the drift or diffusion froman SDEâs temporal marginals, yet the problem of non-identifiability re-mained largely unaddressed [8, 13, 15, 21, 33, 34, 40]. To resolve this is-sue, we provide necessary and sufficient identifiability conditions for linearadditive noise SDEs in Theorem 3.3. With this theorem, we prove thatnon-identifiability of the drift and diffusion arises precisely from the auto-rotational invariance of the initial distribution X0. We also show that thedrift and diffusion parameters identify the underlying causal graph for gen-eral additive noise SDEs. This provides a theoretical foundation for learningthe precise dynamics and causal structure of a linear additive noise process,solely from its temporal marginals, without relying on the usual prior knowl-edge assumptions and model restrictions.To build on our theoretical contributions, we also introduce APPEX, amethod that simultaneously learns drift and diffusion from temporal marginals.APPEX iterates towards the optimal parameter set by alternating betweentrajectory inference, via anisotropic entropy-regularized optimal transport,and MLE parameter estimation. Since each of these subprocedures areasymptotically optimal (Lemma 4.4), given a finite set of the processâ truemarginals, APPEX improves its estimates with each iteration. When testedon hundreds of randomly generated linear additive noise SDEs between di-mensions 3 to 10, APPEX is able to effectively recover the processâ drift,diffusion, and causal graph, while significantly outperforming Waddington-OT, which is one of the most popular inference methods in single-cell biology.Our theoretical and algorithmic contributions offer a number of insights toadvance statistical inference methods from temporal marginals. First, ourframework removes common model restrictions to expand the set of dynam-ics that can be inferred from marginals. This includes negative feedbackloops and repressilator dynamics, which both arise from rotational drift,and correlated noise, which arises from anisotropic diffusion. These features376.1. Limitations and future workare common in real-world datasets, and are especially relevant for accu-rate trajectory and GRN inference from single-cell data [33, 37]. Second,we demonstrate that joint identification of the unique drift and diffusionparameters is theoretically and practically possible in the âmarginals-onlyâsetting. Previous and concurrent works either estimate only one parametergiven the other [8, 13, 15, 21, 33, 38] or estimate both without guaranteeingnon-identifiability [34, 40]. We note that even if a user is only interested inone parameter, APPEX offers the ability to estimate it, regardless of howthe other parameter is initialized. This is particularly significant given that amisspecified diffusion typically leads to a misidentified drift, and vice-versa.Third, we remark that the fundamental condition for identifibility is satisfiedby a large class of distributions and should therefore be widely applicable.Indeed, the only pathological distributions are auto-rotationally invariant,which fall into three categories (see Figure 3.1): 1) satisfying invarianceunder a rotation after scaling and shearing, 2) having a mean-0 Gaussiancross-section, 3) having rank-degenerate support. One can therefore checkthe initial distribution, and perform an intervention if needed, to guaranteesystem identifiability from marginals generated by a linear additive noiseSDE.6.1 Limitations and future workAlthough this paper presents several breakthroughs for the identificationof SDEs from marginal snapshots, we acknowledge several limitations andopportunities for future work.First, while Theorem 3.3 shows that identifiability is guaranteed when theinitial distribution X0 is not auto-rotationally invariant, and that there ex-ist non-identifiable SDEs if X0 is auto-rotationally invariant, we note thatthis is distinct from the stronger claim that given a set of observed tempo-ral marginals from an SDE, the true parameters are uniquely identifiablefrom the data if and only X0 is auto-rotationally invariant. For example,if X0 is rotationally symmetric, e.g. N(0, Id), then it would be impossibleto determine rotational drift. Similarly, if X0 is rank-degenerate, then itwould be impossible to distinguish between drifts that each contain X0 inits nullspace. However, if the underlying SDE has rotational drift, and X0 isdegenerate, then the ensuing marginals likely identify the SDE parameters.Similarly, a full rank rotationally symmetric X0 may identify the underly-ing SDE from the rank degeneracy counterexample. Importantly, the trueSDE is generally unknown, so it is difficult to assess whether identifiabil-386.1. Limitations and future workity holds despite X0 being auto-rotationally invariant. For this reason, thestrong guarantee of identifiability under non-auto-rotationally invariant X0is practical.Second, our identifiability result is restricted to time homogeneous linear ad-ditive noise SDEs. Although this is a popular model for real-world processes,including for GRN inference [8, 19], it would be beneficial to characterizeidentifiability conditions under more general models. We note that char-acterizing identifiability is particularly challenging for non-additive noisemodels, where the noise depends on the process Xt, since the same marginalobservations can be explained by a nonlinear additive noise SDE and alinear multiplicative noise SDE [65]. For this reason, it may be easier toconsider identifiability for parameteric families of nonlinear drift under theadditive noise model. We note that we have already proven results aboutcausal discovery for additive noise SDEs in Section 3.2 and that APPEXis also formulated for general additive noise SDEs. Although the MLE fornonlinear drift may not have a closed form solution, optimizers like BFGS[5, 13, 15] and neural networks [21] have been successfully used for SDEparameter estimation on real and synthetic data.Third, our theory considers the infinite data setting, given the processâ truemarginals {pi}Nâ1i=0 , and we have not studied asymptotic convergence rates.In particular, our identifiability result was established given continuously ob-served marginals, which corresponds to obtaining infinite samples from infi-nite marginals. Similarly, the theory surrounding the optimality of APPEXassumes that we have infinite samples from a finite set of marginals. We leavethe asymptotic analysis of APPEX for future work. Directions for futurework include the sample complexity of entropic optimal transport for trajec-tory inference, which requires asymptotically more samples per marginal aswe decrease the entropic regularization (decreasing time granularity in theSDE setting) [80][Thm. 3], and the fact that the MLE of the AR(1) processfrom multiple trajectories is slightly biased [81]. We note that APPEX couldbe combined with previous methods that relax the marginal constraints (see[13]). This could improve estimates in low data settings, particularly in thebiological setting where taking measurements destroys samples.Fourth, even in the asymptotic setting, the nonconvex nature of APPEXâsalternating projection optimization problem makes the prospect of provingconvergence difficult. Although our experiments indicate that APPEXâsestimates approach the true solution, regardless of the initial guesses forthe drift A and observational diffusion H, we have not yet proved conver-39Code Availabilitygence due to the non-convexity of the optimization space. A formal proof ofconvergence and consistency would further strengthen APPEXâs theoreticalfoundation.Finally, we note that the trajectory inference step of APPEX, which we solveusing Sinkhornâs algorithm, is computationally intensive. When solving thestandard EOT problem (4.7), the asymptotic computation of each iterationof Sinkhornâs algorithm is O(N2ϵ3), where N is the number of samples, and ϵis the scalar entropy regularization [82]. The number of iterations requiredis further related to the average error, which decays like O( 1âNÏµd) [80]. Sim-ilar asymptotics hold for the AEOT problem (4.8). We note that accuratetrajectory inference relies on a large number of samples (N), a small en-tropic regularization ϵ2 â¼ δt, and potentially high-dimensional data. Ourmethod is therefore impractical on very large datasets, which are prevalentin biological applications. It would therefore be advantageous to solve theassociated Schrodinger bridge problem using more efficient methods. Wenote that recent work [83] makes use of closed form solutions to the SBproblem when X0 is a degenerate Dirac distribution, to increase computa-tional efficiency. We also note that the SB problem can be expressed viaforwards and backwards drifts [22], and can therefore be solved via gener-ative diffusion models [84, 85]. Thus, when a closed form solution is notknown, generative diffusion models may be more efficient solvers for thetrajectory inference problem. Concurrent work estimates drift and diffusionusing other methods, including JKO energy minimization [34] and reproduc-ing Hilbert kernel spaces to minimize distance to the Fokker-Planck solution[40]. We have not yet compared the accuracy and computational efficiencyof these methods against APPEX, and a systematic comparative analysiswould be an interesting direction for future research.Code AvailabilityThe Python code for reproducing the experimental results and figures isavailable at https://github.com/guanton/APPEX.40Bibliography[1] Vincent Guan, Joseph Janssen, Hossein Rahmani, Andrew Warren,Stephen Zhang, Elina Robeva, and Geoffrey Schiebinger. Identifyingdrift, diffusion, and causal structure from temporal snapshots, 2025.[2] Albert Einstein. Investigations on the Theory of the Brownian Move-ment. Courier Corporation, 1956.[3] JM MacInnes and FV Bracco. Stochastic particle dispersion modelingand the tracer-particle limit. Physics of Fluids A: Fluid Dynamics,4(12):2809â2824, 1992.[4] Wei Fang and Michael B Giles. Adaptive eulerâmaruyama method forsdes with nonglobally lipschitz drift. The Annals of Applied Probabil-ity, 30(2):526â560, 2020.[5] Yuanyuan Wang, Xi Geng, Wei Huang, Biwei Huang, and MingmingGong. Generator identification for linear sdes with additive and mul-tiplicative noise. Advances in Neural Information Processing Systems,36, 2024.[6] Grigorios A Pavliotis. Stochastic processes and applications. Texts inApplied Mathematics, 60, 2014.[7] Johan du Buisson and Hugo Touchette. Dynamical large deviationsof linear diffusions. Physical Review E, 107(5):054111, 2023.[8] Stephen Y Zhang. Joint trajectory and network inference via referencefitting. arXiv preprint arXiv:2409.06879, 2024.[9] CH Waddington. How animal develop, 1935.[10] William Ogle et al. Aristotle: on the parts of animals. Kegan Paul,French & Company, 1882.41Bibliography[11] Aristotle. Generation of Animals, volume 366 of Loeb Classical Li-brary. Harvard University Press, Cambridge, MA, 1942.[12] Geoffrey Schiebinger, Jian Shu, Marcin Tabaka, Brian Cleary, VidyaSubramanian, Aryeh Solomon, Joshua Gould, Siyan Liu, Stacie Lin,Peter Berube, et al. Optimal-transport analysis of single-cell gene ex-pression identifies developmental trajectories in reprogramming. Cell,176(4):928â943, 2019.[13] Hugo Lavenant, Stephen Zhang, Young-Heon Kim, and GeoffreySchiebinger. Towards a mathematical theory of trajectory inference.arXiv preprint arXiv:2102.09204, 2021.[14] Toshiaki Yachimura, Hanbo Wang, Yusuke Imoto, Momoko Yoshida,Sohei Tasaki, Yoji Kojima, Yukihiro Yabuta, Mitinori Saitou, andYasuaki Hiraoka. scegot: Single-cell trajectory inference frameworkbased on entropic gaussian mixture optimal transport. bioRxiv, pages2023â09, 2023.[15] Le´na¨ıc Chizat, Stephen Zhang, Matthieu Heitz, and GeoffreySchiebinger. Trajectory inference via mean-field langevin in pathspace. Advances in Neural Information Processing Systems, 35:16731â16742, 2022.[16] Atte Aalto, Lauri Viitasaari, Pauliina Ilmonen, Laurent Mombaerts,and Jorge Gonc¸alves. Gene regulatory network inference from sparselysampled noisy data. Nature communications, 11(1):3493, 2020.[17] Philippe Brouillard, Se´bastien Lachapelle, Alexandre Lacoste, SimonLacoste-Julien, and Alexandre Drouin. Differentiable causal discoveryfrom interventional data. Advances in Neural Information ProcessingSystems, 33:21865â21877, 2020.[18] Alejandro Tejada-Lapuerta, Paul Bertin, Stefan Bauer, HananehAliee, Yoshua Bengio, and Fabian J Theis. Causal machine learningfor single-cell genomics. arXiv preprint arXiv:2310.14935, 2023.[19] Martin Rohbeck, Brian Clarke, Katharina Mikulik, Alexandra Pettet,Oliver Stegle, and Kai Ueltzho¨ffer. Bicycle: Intervention-based causaldiscovery with cycles. In Causal Learning and Reasoning, pages 209â242. PMLR, 2024.42Bibliography[20] Paramita Chakraborty, Mark M Meerschaert, and Chae Young Lim.Parameter estimation for fractional transport: A particle-tracking ap-proach. Water resources research, 45(10), 2009.[21] Yunyi Shen, Renato Berlinghieri, and Tamara Broderick. Learning avector field from snapshots of unidentified particles rather than parti-cle trajectories. In ICLR 2024 Workshop on AI4DifferentialEquationsIn Science, 2024.[22] Francisco Vargas, Pierre Thodoroff, Austen Lamacraft, and NeilLawrence. Solving schro¨dinger bridges via maximum likelihood. En-tropy, 23(9):1134, 2021.[23] Yunyi Shen, Renato Berlinghieri, and Tamara Broderick. Multi-marginal schr\u00E2 odinger bridges with iterative reference. arXivpreprint arXiv:2408.06277, 2024.[24] Gaspard Monge. Me´moire sur la the´orie des de´blais et des remblais.De lâImprimerie Royale, 1781.[25] Jean-David Benamou and Yann Brenier. A computational fluid me-chanics solution to the monge-kantorovich mass transfer problem. Nu-merische Mathematik, 84(3):375â393, 2000.[26] Luigi Ambrosio. Optimal transportation and applications. Lecturenotes from the CIME summer school, 2012.[27] Alessio Figalli et al. On the structure of networks of optimal transport.Se´minaire Bourbaki, 3:1â22, 2018.[28] Gabriel Peyre´, Marco Cuturi, et al. Computational optimal trans-port. Center for Research in Economics and Statistics Working Pa-pers, (2017-86), 2017.[29] Gabriel Peyre´, Marco Cuturi, et al. Computational optimal trans-port: With applications to data science. Foundations and Trends®in Machine Learning, 11(5-6):355â607, 2019.[30] Leonid Kantorovich. On the translocation of masses. Dokl. Akad.Nauk. SSSR, 1941.[31] Luigi Ambrosio, Nicola Fusco, and Diego Pallara, editors. OptimalTransportation and Applications: Lectures given at the CIME SummerSchool, volume 1813 of Lecture Notes in Mathematics. Springer, 2001.43Bibliography[32] Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimaltransport. In Advances in Neural Information Processing Systems 26(NIPS 2013), pages 2292â2300, 2013.[33] Caleb Weinreb, Samuel Wolock, Betsabeh K Tusi, Merav Socolovsky,and Allon M Klein. Fundamental limits on dynamic inference fromsingle-cell snapshots. Proceedings of the National Academy of Sciences,115(10):E2467âE2476, 2018.[34] Antonio Terpin, Nicolas Lanzetti, and Florian Do¨rfler. Learning dif-fusion at lightspeed. arXiv preprint arXiv:2406.12616, 2024.[35] Suryanarayana Maddu, Victor Charde`s, Michael Shelley, et al. Infer-ring biological processes with intrinsic noise from cross-sectional data.arXiv preprint arXiv:2410.07501, 2024.[36] Søren Wengel Mogensen and Niels Richard Hansen. Graphical mod-eling of stochastic processes driven by correlated noise. Bernoulli,28(4):3023â3050, 2022.[37] Augusto Santos, Diogo Rente, Rui Seabra, and Jose´ MFMoura. Learn-ing the causal structure of networked dynamical systems under latentnodes and structured noise. In Proceedings of the AAAI Conferenceon Artificial Intelligence, volume 38, pages 14866â14874, 2024.[38] Aden Forrow. Consistent diffusion matrix estimation from populationtime series. arXiv preprint arXiv:2408.14408, 2024.[39] Keith J Beven, D Ed Henderson, and Alison D Reeves. Dispersion pa-rameters for undisturbed partially saturated soil. Journal of hydrology,143(1-2):19â43, 1993.[40] Luc Brogat-Motte, Riccardo Bonalli, and Alessandro Rudi. Learn-ing controlled stochastic differential equations. arXiv preprintarXiv:2411.01982, 2024.[41] Jaya PN Bishwal. Parameter estimation in stochastic differential equa-tions. Springer, 2007.[42] Alexandros Beskos, Omiros Papaspiliopoulos, Gareth O Roberts, andPaul Fearnhead. Exact and computationally efficient likelihood-basedestimation for discretely observed diffusion processes (with discussion).Journal of the Royal Statistical Society Series B: Statistical Method-ology, 68(3):333â382, 2006.44Bibliography[43] Jan Nygaard Nielsen, Henrik Madsen, and Peter C Young. Parameterestimation in stochastic differential equations: an overview. AnnualReviews in Control, 24:83â94, 2000.[44] Shahab Torkamani and Eric A Butcher. Stochastic parameter esti-mation in nonlinear time-delayed vibratory systems with distributeddelay. Journal of Sound and Vibration, 332(14):3404â3418, 2013.[45] Kevin McGoff, Sayan Mukherjee, and Natesh Pillai. Statistical infer-ence for dynamical systems: A review. Statist. Surv., 2015.[46] SC Kou, Benjamin P Olding, Martin Lysy, and Jun S Liu. A mul-tiresolution method for parameter estimation of diffusion processes.Journal of the American Statistical Association, 107(500):1558â1574,2012.[47] Richard Nickl and Kolyan Ray. Nonparametric statistical inferencefor drift vector fields of multi-dimensional diffusions. The Annals ofStatistics, 48(3):1383â1408, 2020.[48] Peter Craigmile, Radu Herbei, Ge Liu, and Grant Schneider. Statis-tical inference for stochastic differential equations. Wiley Interdisci-plinary Reviews: Computational Statistics, 15(2):e1585, 2023.[49] Georg Manten, Cecilia Casolo, Emilio Ferrucci, Søren Wengel Mo-gensen, Cristopher Salvi, and Niki Kilbertus. Signature kernel condi-tional independence tests in causal discovery for stochastic processes.arXiv preprint arXiv:2402.18477, 2024.[50] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, andGeorge Em Karniadakis. Learning nonlinear operators via deeponetbased on the universal approximation theorem of operators. Naturemachine intelligence, 3(3):218â229, 2021.[51] Rich Pawlowicz, Charles Hannah, and Andy Rosenberger. Lagrangianobservations of estuarine residence times, dispersion, and trapping inthe salish sea. Estuarine, Coastal and Shelf Science, 225:106246, 2019.[52] Ahmed Nafidi and IËlyasse Makroz. A comparison of methods for es-timating parameters of the stochastic lomax process: through sim-ulation study. Hacettepe Journal of Mathematics and Statistics,53(2):495â505, 2024.45Bibliography[53] Cole Trapnell, Davide Cacchiarelli, Jonna Grimsby, Prapti Pokharel,Shuqiang Li, Michael Morse, Niall J Lennon, Kenneth J Livak, Tarjei SMikkelsen, and John L Rinn. The dynamics and regulators of cellfate decisions are revealed by pseudotemporal ordering of single cells.Nature biotechnology, 32(4):381â386, 2014.[54] Manu Setty, Michelle D Tadmor, Shlomit Reich-Zeliger, Omer An-gel, Tomer Meir Salame, Pooja Kathail, Kristy Choi, Sean Bendall,Nir Friedman, and Dana Peâer. Wishbone identifies bifurcating de-velopmental trajectories from single-cell data. Nature biotechnology,34(6):637â645, 2016.[55] Jeffrey A Farrell, Yiqun Wang, Samantha J Riesenfeld, KarthikShekhar, Aviv Regev, and Alexander F Schier. Single-cell reconstruc-tion of developmental trajectories during zebrafish embryogenesis. Sci-ence, 360(6392):eaar3131, 2018.[56] CH Waddington. The strategy of the genes: a discussion of someaspecs of theoretical biology. London: Allen and Unwin, 1957.[57] Chuanjian Man and Christina W Tsai. Stochastic partial differentialequation-based model for suspended sediment transport in surface wa-ter flows. Journal of engineering mechanics, 133(4):422â430, 2007.[58] Mo´nica Basilio Hazas, Francesca Ziliotto, Jonghyun Lee, MassimoRolle, and Gabriele Chiogna. Evolution of plume geometry, dilutionand reactive mixing in porous media under highly transient flow fieldsat the surface water-groundwater interface. Journal of ContaminantHydrology, 258:104243, 2023.[59] Peter Salamon, Daniel Ferna`ndez-Garcia, and JJ Go´mez-Herna´ndez.Modeling tracer transport at the made site: The importance of het-erogeneity. Water Resources Research, 43(8), 2007.[60] Amro MM Elfeki. Prediction of contaminant plumes (shapes, spatialmoments and macrodispersion) in aquifers with insufficient geologicalinformation. Journal of Hydraulic Research, 44(6):841â856, 2006.[61] E Eric Adams and Lynn W Gelhar. Field study of dispersion in aheterogeneous aquifer: 2. spatial moments analysis. Water ResourcesResearch, 28(12):3293â3307, 1992.46Bibliography[62] J Mark Boggs, Steven C Young, Lisa M Beard, Lynn W Gelhar, Ken-neth R Rehfeldt, and E Eric Adams. Field study of dispersion in aheterogeneous aquifer: 1. overview and site description. Water Re-sources Research, 28(12):3281â3291, 1992.[63] DM Mackay, DL Freyberg, PV Roberts, and JA Cherry. A natu-ral gradient experiment on solute transport in a sand aquifer: 1. ap-proach and overview of plume movement. Water Resources Research,22(13):2017â2029, 1986.[64] Vanessa Nenna, Adam Pidlisecky, and Rosemary Knight. Applicationof an extended kalman filter approach to inversion of time-lapse electri-cal resistivity imaging data for monitoring recharge. Water ResourcesResearch, 47(10), 2011.[65] Megan A Coomer, Lucy Ham, and Michael PH Stumpf. Noise distortsthe epigenetic landscape and shapes cell-fate decisions. Cell Systems,13(1):83â102, 2022.[66] Philip Boeken and Joris M Mooij. Dynamic structural causal models.arXiv preprint arXiv:2406.01161, 2024.[67] James Clerk Maxwell. Ii. illustrations of the dynamical theory of gases.The London, Edinburgh, and Dublin Philosophical Magazine and Jour-nal of Science, 20(130):21â37, 1860.[68] PETER McCULLAGH. Tensor notation and cumulants of polynomi-als. Biometrika, 71(3):461â476, 1984.[69] Jo´zef Marcinkiewicz. Sur une proprie´te´ de la loi de gauss. Mathema-tische Zeitschrift, 44(1):612â618, 1939.[70] Niels Hansen and Alexander Sokol. Causal interpretation of stochasticdifferential equations. Electronic Journal of Probability, 2014.[71] Yiheng Liu, Elina Robeva, and Huanqing Wang. Learning linear non-gaussian graphical models with multidirected edges. Journal of CausalInference, 9(1):250â263, 2021.[72] Lars Lorch, Andreas Krause, and Bernhard Scho¨lkopf. Causal model-ing with stationary diffusions. In International Conference on Artifi-cial Intelligence and Statistics, pages 1927â1935. PMLR, 2024.47Bibliography[73] Hirotogu Akaike. Information theory and an extension of the maxi-mum likelihood principle. In Selected papers of Hirotugu Akaike, pages199â213. Springer, 1998.[74] Halbert White. Maximum likelihood estimation of misspecified mod-els. Econometrica: Journal of the econometric society, pages 1â25,1982.[75] Jean-David Benamou, Guillaume Carlier, Simone Di Marino, and LucaNenna. An entropy minimization approach to second-order variationalmean-field games. Mathematical Models and Methods in Applied Sci-ences, 29(08):1553â1583, 2019.[76] Hicham Janati, Boris Muzellec, Gabriel Peyre´, and Marco Cuturi.Entropic optimal transport between unbalanced gaussian measures hasa closed form. Advances in neural information processing systems,33:10468â10479, 2020.[77] Promit Ghosal, Marcel Nutz, and Espen Bernton. Stability of en-tropic optimal transport and schro¨dinger bridges. Journal of Func-tional Analysis, 283(9):109622, 2022.[78] Jakob Runge. Necessary and sufficient graphical conditions for opti-mal adjustment sets in causal graphical models with hidden variables.Advances in Neural Information Processing Systems, 34:15762â15773,2021.[79] Alexander Reisach, Christof Seiler, and Sebastian Weichwald. Bewareof the simulated dag! causal discovery benchmarks may be easy togame. Advances in Neural Information Processing Systems, 34:27772â27784, 2021.[80] Aude Genevay, Le´naic Chizat, Francis Bach, Marco Cuturi, andGabriel Peyre´. Sample complexity of sinkhorn divergences. In The22nd international conference on artificial intelligence and statistics,pages 1574â1583. PMLR, 2019.[81] Johannes Ledolter. Estimation bias in the first-order autoregressivemodel and its impact on predictions and prediction intervals. Com-munications in Statistics-Simulation and Computation, 38(4):771â787,2009.48Bibliography[82] Jason Altschuler, Jonathan Niles-Weed, and Philippe Rigollet. Near-linear time approximation algorithms for optimal transport viasinkhorn iteration. Advances in neural information processing sys-tems, 30, 2017.[83] Hanwen Huang. One-step data-driven generative model via schr\u00E2odinger bridge. arXiv preprint arXiv:2405.12453, 2024.[84] Valentin De Bortoli, James Thornton, Jeremy Heng, and ArnaudDoucet. Diffusion schro¨dinger bridge with applications to score-basedgenerative modeling. Advances in Neural Information Processing Sys-tems, 34:17695â17709, 2021.[85] Yuyang Shi, Valentin De Bortoli, Andrew Campbell, and ArnaudDoucet. Diffusion schro¨dinger bridge matching. Advances in NeuralInformation Processing Systems, 36, 2024.[86] Tatsunori Hashimoto, David Gifford, and Tommi Jaakkola. Learn-ing population-level diffusions with generative rnns. In InternationalConference on Machine Learning, pages 2417â2426. PMLR, 2016.[87] Stephan Mandt, Matthew Hoffman, and David Blei. A variationalanalysis of stochastic gradient algorithms. In International conferenceon machine learning, pages 354â363. PMLR, 2016.[88] Simo Sa¨rkka¨ and Arno Solin. Applied stochastic differential equations,volume 10. Cambridge University Press, 2019.[89] Bernt Oksendal. Stochastic differential equations: an introduction withapplications. Springer Science & Business Media, 2013.[90] Miranda Holmes-Cerfon. Applied stochastic analysis, 2015.[91] Daniel W Stroock. Partial differential equations for probabalists [sic].Cambridge University Press, 2008.[92] Vanessa Didelez. Graphical models for marked point processes basedon local independence. Journal of the Royal Statistical Society SeriesB: Statistical Methodology, 70(1):245â264, 2008.[93] Søren Wengel Mogensen and Niels Richard Hansen. Markov equiva-lence of marginalized local independence graphs. The Annals of Statis-tics, 48(1):539â559, 2020.49[94] Kaare Brandt Petersen, Michael Syskind Pedersen, et al. The matrixcookbook. Technical University of Denmark, 7(15):510, 2008.[95] Richard EWilliamson and Hale F Trotter. Multivariable mathematics:linear algebra, differential equations, calculus. (No Title), 1974.[96] L Campbell. Reduction theorems for the strong real jacobian con-jecture. In Annales Polonici Mathematici, volume 110, pages 1â11.Institute of Mathematics Polish Academy of Sciences, 2014.[97] Hicham Janati. Advances in Optimal transport and applications toneuroscience. PhD thesis, Institut Polytechnique de Paris, 2021.[98] Marcel Nutz. Introduction to entropic optimal transport. Lecturenotes, Columbia University, 2021.[99] Thomas M Cover. Elements of information theory. John Wiley &Sons, 1999.[100] S Peiris, R Mellor, and P Ainkaran. Maximum likelihood estimationfor short time series with replicated observations: A simulation study.InterStat, 9:1â16, 2003.50Appendix AAppendixA.1 Non-identifiable linear additive noise SDEsWe summarize three classical examples of non-identifiability from the liter-ature [5, 13, 23, 33, 86]. In all examples, the processes Xt and Yt share thesame time marginals.Example A.1 (Starting at stationary distribution [13, 23]). This exampleis motivated by the fact that multiple SDEs share the same stationary distri-bution. Hence, if p0 is initialized at the stationary distribution, then theseSDEs will be non-identifiable from their marginals. We provide a simpleexample.dXt = âXt dt+ dWt, X0 â¼ N (0, 12) (A.1)dYt = â10Yt dt+â10dWt, Y0 â¼ N (0, 12) (A.2)Here, both SDEs have the same stationary distribution N (0, 12) despite hav-ing significantly different individual trajectories. Indeed, the stationary dis-tribution of a 1-dimensional 0-mean Ornstein-Uhlenbeck (OU) process withdrift âÎ»Xt and diffusion Ï is N (0, Ï22Î»), which depends only on the drift/diffu-sivity ratio Î»/Ï2 [13]. For multivariate OU processes with drift A and ob-servational diffusion H, the stationary distribution N (0,Σ) depends only onthe relationship ΣA + AΣ = âH [87]. Intuitively, the set of OU processes,which share the same stationary distribution N (0,Σ) are equivalent up torescaling time. We note that this example also generalizes easily to the casewhere one component of X0 is mean-0 Gaussian, i.e. (PX0)1 â¼ N (0, Ï2).Example A.2 (Rotation around process mean [23, 33, 86]). In this ex-ample, the first SDE governing Xt is 2-dimensional Brownian motion andthe second SDE governing Yt adds a divergence-free rotational vector field(x, y)â (y,âx) about (0, 0), which is undetectable if the initial distribution51A.1. Non-identifiable linear additive noise SDEsp0 is rotationally invariant.dXt = dWt, X0 â¼ N (0, Id) (A.3)dYt =[0 1â1 0]Yt dt+ dWt, Y0 â¼ N (0, Id) (A.4)Hence, if p0 is an isotropic distribution with mean (0, 0), then Xt and Ytare non-identifiable from one another [23]. This can also be shown directlywith the Fokker-Planck equation (1.4), since â·(Ax)p(x, t) = âp(x, t) Â·Ax+p(x, t)â · A(x) = âp(x, t) · Ax = 0, if p is parallel to the rotational vectorfield Ax =[0 1â1 0]x.Example A.3 (Degenerate rank [5]). In this example motivated by [5],trajectories (and therefore temporal marginals) of distinct SDEs are non-identifiable if the process is degenerate, limiting observations to a strict sub-space of Rd.dXt =[1 21 0]Xt dt+[1 2â1 â2]dWt, X0 =[1â1](A.5)dYt =[1/3 4/32/3 â1/3]Yt dt+[1 2â1 â2]dWt, Y0 =[1â1](A.6)The drift matrices[1 21 0],[1/3 4/32/3 â1/3]each have eigenvector X0 =[1â1]with eigenvalue â1. Moreover, the diffusion is rank-degenerate with col-umn space span([1â1]). Thus, both SDEs will have identical behaviouralong span([1â1]) and are only differentiated by their behaviour elsewhere.Given X0 =[1â1], the processes will stay within span([1â1]) and are non-identifiable from one another. We note that this non-identifiability holdseven when we observe trajectories [5], whereas the first two examples areidentifiable from trajectories but not identifiable from marginals.Theorem 3.3 proves that identifiability in the âmarginals-onlyâ setting isguaranteed for linear additive noise SDEs if X0 is not auto-rotationally in-variant. Figure A.1 demonstrates how non-auto-rotationally invariant initialdistributions X0 resolves these classical non-identifiable examples from theliterature.52A.2. Properties of the Linear Additive Noise SDE2 1 0 1 2X0.00.20.40.60.81.01.21.4DensityA=-1 G=[1.]2 1 0 1 2X0.00.20.40.60.81.01.21.4DensityA=-10 G=[ 10]Marginal at time 0.2 (a) This corresponds toExample A.1 with X0 â¼ δ1. TheSDE with drift A = H = 10(right) has converged to thestationary distribution (red) attime t = 0.2, while the SDE withA = H = 1 (left) has not yetsufficiently drifted or diffused.(b) This corresponds toExample A.2, withX0 â¼ Unif{(2, 0), (2, 0.1)}. Thefirst SDE is irrotational andmaintains mean (2, 0.05). Thesecond SDE (right) rotatesabout the origin. At timet = 1.5, the mean has shifted toapproximately (0,â2).(c) This corresponds toExample A.3, withX0 â¼ Unif{(1, 0)â¤, (0, 1)â¤}. Thefirst SDE (left) has eigenvalue 2for (2, 1), while the second(right) has eigenvalue 1, leadingto more extreme intercepts inthe left plot at t = 0.8.Figure A.1: Marginals are plotted at various times for the SDE pairs fromSection 3 of the Appendix, following various non-auto-rotationally invariantinitial distributions.A.2 Properties of the Linear Additive Noise SDELet Xt evolve according to a linear additive noise SDE (1.3). Then itssolution can be obtained using the integrating factor eâAt [88][Sec 4.3] andis given byXt = eAtX0 +â« t0eA(tâs)GdWs (A.7)By the independent increments of Wt, we may interpret Xt as being com-posed of two independent r.v.âs. The first term, eAtX0, is a time-dependent53A.3. Weak formulation of the Fokker-Planck equationfunction of its initial distribution. The second term, given by the Ito inte-gral,â« t0 eA(tâs)GdWs â¼ N (0,Σt), is a mean-0 Gaussian with time-dependentvariance Σt =â« t0 eA(tâs)HeAâ¤(tâs)ds. Indeed, since the integrand is deter-ministic, Gaussianity follows from [89][Theorem A.7], and the mean andcovariance respectively follow from the non-anticipating property and Itoâsisometry.It follows that the exact transition kernel of a linear additive noise SDE satis-fies Xt+ât|Xt â¼ N (µt,Σt), where µt = eAâtXt[8]. The linearized transitionkernel is given by N (Xt +AXtât,Hât).A.3 Weak formulation of the Fokker-PlanckequationWe may define the linear operator L(Lf)(x) = b(x) · f(x) + 12â · (Hâf(x)),L is the generator of the additive noise SDE (1.2) [90]. Its adjoint operatoris given by(Lâg)(x) = ââ · (b(x, t)g(x)) +â · (Hâg(x)).and satisfiesâ«RdLf(x)g(x) dx = â¨Lf, gâ© = â¨f,Lâgâ© =â«Rdf(x)Lâg(x) dx,for all smooth and compactly supported test functions f, g â Câc (Rd).Then, the Fokker-Planck equation can be formulated in the weak sense asfollows. For each t â¥ 0 let pt be a probability measure on Rd and assumethat t 7â pt is narrowly continuous. Then we say pt solves the Fokker-Planckequation in the weak sense provided thatddtâ«Ï(x) dpt(x) =â«LÏ(x) dpt(x) âÏ â C2c (Rd). (A.8)This weak formulation accommodates the case where b is not differentiableand pt is not necessarily absolutely continuous with respect to the Lebesguemeasure. For pt with sufficiently regular Lebesgue density p(x, t) this weak54A.4. Properties of auto-rotationally invariant r.v.ssolution concept is equivalent to the classical sense of the Fokker-Planckequation, i.e.âtp(x, t) = Lâp(x, t).For the measure-theoretic weak formulation of the Fokker-Planck equation,we quote the following existence and uniqueness results from [91]. Underour standing assumptions on the drift and diffusion coefficients (namely:the drift is linear, and the diffusion is constant but possibly degenerate),by [91] Theorem 1.1.9 weak solutions to the Fokker-Planck equation existwhenever the initial condition Ï0 has finite second moments; and moreover,by Theorem 2.2.9 in [91], the solution is unique provided that all momentsof Ï0 are finite.A.4 Properties of auto-rotationally invariant r.v.sLemma A.4. Let Σ â» 0 and A â RdÃd be such thatAΣ+ Î£Aâ¤ = 0.Then, B = Î£1/2AÎ£â1/2 is skew-symmetric, and for all t â R,eAt = Î£â1/2 eBtΣ1/2.In other words, the Σ-generalized rotation eAt produces a classical rotationeBt after the change of variables y = Σ1/2x, which can be interpreted as ashearing and scaling transformation of Rd.Proof. Since Σ â» 0, the unique symmetric positive definite square root Σ1/2exists and is invertible. DefineB = Î£â1/2AÎ£1/2, Bâ¤ = Î£1/2Aâ¤Î£â1/2.Then, since AΣ + Î£Aâ¤ = 0, we can multiply on the left by Î£â1/2 and onthe right by Î£â1/2 to yield0 = Î£â1/2AÎ£1/2 +Î£1/2Aâ¤Î£â1/2 = B +Bâ¤,which proves that B is skew-symmetric, and therefore eBt defines a classicalrotation in Rd, such thateBt = eÎ£â1/2AÎ£1/2 = Î£â1/2eAtÎ£1/2,55A.4. Properties of auto-rotationally invariant r.v.swhere the last equality can be seen by writing eBt =âân=0 tnBnn! andnoting that Bn = (Î£â1/2AÎ£1/2)n = Î£â1/2AnÎ£1/2. Rearranging yieldseAt = Î£â1/2 eBtΣ1/2 as desired.To explicitly see that the change of variable y = Σ1/2x amounts to a shear-ing and scaling, note that any positive definite matrix admits an LDLâ¤factorization:Î£ = LDLâ¤,where L is unit lower-triangular (so its off-diagonal entries describe a sheartransformation) and D is a diagonal matrix with positive entries (represent-ing scaling along the coordinate axes). Thus, the transformation x 7â y =Σ1/2x decomposes into a shear (via L) and a scaling (via D1/2).Lemma A.5. Let X be a r.v. with covariance Σ and suppose that eAtX â¼X. Then, AΣ+ Î£Aâ¤ = 0.Proof. Let C(t) = Cov(eAtX). eAtX â¼ X implies that for all t â¥ 0, C(t) =Cov(X) = Σ. Hence,C(t) = eAtÎ£eAâ¤t = ΣC â²(t) = AeAtÎ£eAâ¤t + eAtÎ£Aâ¤eAâ¤t = 0C â²(0) = AΣ+ Î£Aâ¤ = 0.Lemma A.6. Let Cov(X) = Σ â» 0 and suppose that eAtX â¼ X for someA Ì¸= 0. Then X admits a density, and the level sets of its probability densityfunction are Σ-rotationally invariant ellipsoids:Let E(k)Σ = {x â Rd : xâ¤Î£â1x = k}, then âx, y â E(k)Σ , p(x) = p(y) (A.9)Proof. Since eAtX â¼ X with A Ì¸= 0, it follows that p(x) = p(eAtx) for allx â Rd. Since eAtE(k)Σ = E(k)Σ , the probability distribution of X must beconstant within each ellipsoid E(k)Σ . Since Σ â» 0, then each set E(k)Σ is non-degenerate in Rd, thus inducing a probability density function on X.Lemma A.7. Let X be a d dimensional r.v. such that dim(span(X)) < dwith probability 1. Then,eAtX â¼ Xadmits a non-trivial solution A Ì¸= 0.56A.5. Additional proofs for Theorem 3.3Proof. We expandeAtX = X +( âân=1tnn!An)XHence, a sufficient condition for the claim is that X is in the nullspace of Awith probability 1. Since dim(span(X)) < d a.s., then by the rank-nullitytheorem, there exists A such that AX = 0 a.s.Proposition A.8. Let X0 â¼ p0 be a discrete random variable. Then X0 isnot auto-rotationally invariant if and only if the support of p0 spans Rd.Proof. Lemma A.7 provides the only if direction.Now, let span(X0) = Rd a.s. Since eAt represents a continuous transforma-tion, and X0 is supported on a discrete set of points, which cannot all bemapped to 0 by A, it follows that the support would shift for each t > 0.We conclude that eAtX0 Ì¸â¼ X0.A.5 Additional proofs for Theorem 3.3Lemma A.9. Let AΣ+ ΣAT = 0. Then, for all Ï â¥ 0â« t0eA(tâs)Î£eAT (tâs)ds = tΣ. (A.10)Proof. Let Ï = tâ s and set M(Ï) = eAÏÎ£eAT Ï . Then,ddÏM(Ï) = AeAÏÎ£eAT Ï + eAÏÎ£eAT ÏAT (A.11)= eAÏAÎ£eAT Ï + eAÏÎ£AT eAT Ï (A.12)= eAÏ (AΣ+ ΣAT )eAT Ï = 0 (A.13)Since M(0) = Σ, it follows that M(Ï) = Σ for all Ï . We therefore concludethat â« t0eA(tâs)Î£eAT (tâs)ds =â« t0Σ = tΣ (A.14)57A.5. Additional proofs for Theorem 3.3Let X0 beauto-rotationally invariantIs Cov(X0) rank degenerate?X0 â null(A)for A Ì¸= 0Example A.3(rank degeneracy)is non-identifiableYesDoes there exist A Ì¸= 0s.t. eAtX0 â¼ X0 ?X0 has a densitywith ellipsoidal level setsExample A.2(rotation)is non-identifiableYes(PX0)1 â¼ N (0, Ï2)Example A.1(stationary start)is non-identifiableNoNoFigure A.2: The three types of auto-rotational invariance is illustrated. Eachtype leads to non-identifiability from marginals, when initialized at the dis-tribution.Lemma A.10. Let Σ âª° 0 be a symmetric positive semidefinite matrix andsuppose that A is skew-symmetric with respect to Σ, i.e.AΣ+ Î£Aâ¤ = 0,then without loss of generality, we may define A with Tr(A) = 0.Proof.Tr(AΣ+ Î£Aâ¤) = Tr(0) = 0.58A.5. Additional proofs for Theorem 3.3Using the linearity of the trace and the cyclic property Tr(XY ) = Tr(Y X):Tr(AΣ+ Î£Aâ¤) = Tr(AΣ) + Tr(Aâ¤Î£)= Tr(AΣ) + Tr(ΣA)= 2Tr(AΣ).Therefore,2Tr(AΣ) = 0 â Tr(AΣ) = 0.If Σ â» 0, then we immediately conclude Tr(A) = 0. Otherwise, we maypick A such that Σi,i = 0 =â Ai,i = 0 to ensure that Tr(A) = 0.Corollary A.11. If A is skew-symmetric with respect to Σ, and Σ â» 0,then Ax is a divergence-free vector field:â Â·Ax = 0Proof. We know thatâ Â·Ax =dâi=1ââxi(Aix)=dâi=1ââxidâi=1(Ai,jxj)=dâi=1Ai,i = Tr(A),which is 0 by Lemma A.10.Proposition A.12. Let A Ì¸= 0 such that eAtX â¼ X and let H = GGâ¤satisfy AH +HAâ¤ = 0. Then, for all γ â R, the SDEsdXt = γGdWt X0 â¼ X (A.15)dYt = AYt dt+ γGdWt Y0 â¼ X (A.16)will have the same marginals Xt â¼ Yt ât â¥ 0.59A.6. SDEs and causalityProof. The Fokker-Planck equations for the two SDEs are given byââtp(x, t) =Î³22â · (Hâp(x, t)) p(x, 0) = p0(x),(A.17)ââtp(x, t) = ââ · ((Ax)p(x, t)) + Î³22â · (Hâp(x, t)) p(x, 0) = p0(x).(A.18)To show that the processes exhibit the same time marginals, we show thatthe equations (A.17) and (A.18) are equivalent. We argue in the case wherep(x, t) is a general probability measure, via the weak formulation of theFokker Planck equation, given in equation (A.8) in Section A.3. It sufficesto show that for all Ï â C2c (Rd), we have thatâ«(âÏ(x) Â·Ax) dpt(x) = 0since this corresponds to the weak formulation of the divergence term â Â·((Ax)p(x, t)). Integrating by parts, we see thatâ«(âÏ(x) Â·Ax) dpt(x) = ââ«Ï(x)â · (Ax) dpt(x).By Corollary A.11, it follows thatâ·(Ax) = 0. Hence, â« (âÏ(x)·Ax) dpt(x) =0, as desired.A.6 SDEs and causalityLet V be a set of d endogenous variables {X(1), X(d)} and W be a set ofexogenous variables {X(d+1), X(d+r)}, which are excluded from the prin-cipal model, but may influence the endogenous variables. Suppose thatall variables evolve according to a general SDE, which can be written inâintegrand-integratorâ formdXt = a(t,Xt) dh(t,Xt). (A.19)By integrating and considering each variable X(j)t in the endogenous setj â V , [66] showed that we can define a dynamic structural causal model(DSCM) [Definition 1]:X(j)t = X(j)0 +â« t0aj(s,Xα(j)s ) dhj(s,Xβ(j)s ), (A.20)60A.6. SDEs and causalitywhere α(j), β(j) â V âªW represent the variables that are used as argumentsfor aj(t,Xt) and hj(t,Xt) respectively. In particular, (A.20) defines a causalgraph G with vertices V âªW , and edges E = {iâ j : j â V, i â Î±(j)âªÎ²(j)}[66][Definition 2]. In other words, we only consider edges that point towardsvariables in the endogenous set, and we include the edge i â j if X(i)influences the evolution of endogenous variable X(j) through the integrandaj(t,Xt) or integrator hj(t,Xt). Additionally, in [66], self-loops j â j areexcluded to have a simple SCM, but we do not follow this convention in thispaper. Under this more general framework, we see that latent confounderscan affect endogenous variables via the integrand or integrator, which can bedecoupled as drift and diffusion. For instance, any SDE driven by BrownianmotiondXt = b(t,Xt) dt+ Ï(t,Xt)dWt, (A.21)can be rewritten in the form (A.19) by setting a(t,Xt) =[b(t,Xt) Ï(t,Xt)]and h(t,Xt) =[t Wt]â¤. Given that Xt â Rd, Wt â Rm, we have corre-sponding dimensions b(Xt) â Rd, Ï(Xt) â RdÃm, Zt â Rm+1, and a(Xt) âRdÃ(m+1). The DSCM can therefore also be written asX(j)t = X(j)0 +â« t0bj(s,Xα(j)s )ds+â« t0Ï(s,XÎ²(j)s )dWs. (A.22)DSCMs observe graphical Markov properties under local independence [36,92, 93] under the assumption of independent driving noise, also known asindependent integrators [66][Theorem 4]. This is formally defined as β(j) âW and β(i) â© β(j) = â for all i, j â V .Although the observational diffusion H = GGâ¤ informs the presence of la-tent confounders between pairs of variables, we note that it is impossiblefor H alone to determine general multidirected edges (i1, · · · , ip) over p > 2components, which indicate an unobserved confounder of X(i1)t , · · · , X(ip)t .Indeed, given only H = GGâ¤, multiple causal interpretations may be pos-sible, because structurally different matrices G, GË can obey H = GGâ¤ =GËGËâ¤.Example A.13 (Non-identifiability of higher order multidirected edges).Consider two matrices G1, G2, which share the same observational diffusionH:G1 =ï£®ï£°0 1 11 0 11 1 0ï£¹ï£» , G2 =ï£®ï£°43 13 13134313131343ï£¹ï£» , H = G1Gâ¤1 = G2Gâ¤2 =ï£®ï£°2 1 11 2 11 1 2ï£¹ï£»61A.6. SDEs and causalityGiven G1 as the additive noise of an SDE, each pair of variables shares anoise source, since each row shares a nonzero column entry with anotherrow. However, there is no common noise source that is shared among all 3variables, since each column contains a 0. Hence, the causal interpretationunder G1 consists in three bidirected edges 1â 2, 1â 3, 2â 3. In contrast,the causal interpretation under G2 consists in a single multidirected edge(1, 2, 3), since all components share noise sources.Since G1 and G2 admit different causal graphs, despite having the sameobservational diffusion H, this shows that H provides information about theexistence of unobserved confounders between pairs of variables, but cannotprovide further causal structure about the confounder with respect to theother endogenous variables [d] {i, j}.We note that if the noise is not additive, then this further complicates thecausal interpretation. In this case, the driving noise Ï(Xt) may be a functionof the endogenous variables, i.e. β(j) â© V Ì¸= â . Thus, unlike the additivenoise setting in Theorem 3.5, simple edges i â j may be informed by thedriving noise, via β(j), rather than just the drift, via α(j). However, sinceonly Ï(Xt)Ï(Xt)â¤ is observable rather than Ï(Xt) itself, observational dataunder such a model would lend itself to multiple interpretations of the causalgraph. The idea is similar to Example A.13, where we saw one interpretationfeature three bidirected edges and another feature a multi-edge, but undernon-additive noise, the ambiguity extends to simple edges i â j. This isillustrated in Example 5.5 in [70].The causal framework also highlights the impracticality of common modelassumptions within the scRNA-seq literature, namely isotropic diffusionH = GGâ¤ = Ï2Id and irrotational drift b(Xt) = âÏ. We have seen inLemma 3.5 that latent confounders from shared noise can only be modeledwith anisotropic diffusion. Moreover, while we expect GRNs to contain feed-back loops, i.e. cycles in G, imposing irrotational drift prevents most cyclesfrom being considered. Lemma A.14 shows that imposing irrotationality ona linear model is equivalent to imposing a symmetric drift, which determinesa symmetric causal graph. However, as illustrated in [33][Fig. 6], symmetricGRNs cannot capture important relationships, such as negative feedbackloops or repressilator dynamics. Indeed, consider the 2-cycle i â j. ByLemma 3.5, i â j if Aj,i Ì¸= 0 and j â i if Ai,j Ì¸= 0. However, if the modelis irrotational, then by Proposition A.14, Ai,j = Aj,i, which can only modela positive feedback loop.Lemma A.14. Let Xt evolve according to a linear additive noise SDE (1.3).62A.7. Entropy-regularized Optimal TransportThen, the drift AXt is irrotational if and only if A is symmetric, i.e. A =Aâ¤.Proof. First, we note that for any symmetric matrix A, the vector fieldAx can be expressed as the gradient of the scalar potential given by thequadratic form Ï(x) = xâ¤Ax/2 [94][(96)-(97)]. Hence, Ax is irrotational ifA = Aâ¤.For the opposite direction, we recall the fact that any gradient field fromRn â Rn has a symmetric Jacobian [95, 96]. Since the Jacobian of a linearvector field Ax : Rn â Rn is A, we conclude that all irrotational linearvector fields from Rn â Rn are given by symmetric matrices.A.7 Entropy-regularized Optimal TransportThe entropy-regularized optimal transport problem also admits a dual for-mulation, in terms of finding a pair of potentials (f, g) with respect to theGaussian transition kernel K(x, y) â e c(x,y)ϵ2 â eââ¥yâxâ¥22Ïµ2 [97, 98]:Ïâ(x, y) = efâ(x)+gâ(y)K(x, y) dµ(x)dν(y) (A.23)fâ, gâ = supfâL1(µ),gâL1(Î½)Eµ(f) + EÎ½(g)â(â«RdÃRdef(x)+g(y)K(x, y) dµ(x)dÎ½(y)â 1). (A.24)In particular, given marginals µ, ν and the transition kernel K, Sinkhornâsalgorithm uses the dual formulation (A.24) to find the ϵ-entropy regularizedoptimal transport solution Ïâ, by solving for fâ, gâ via iterative projections[29].Definition A.15 (Markov concatenation of couplings). Given Polish spacesX1, X2, X3 and couplings Ï1,2 â P(X1 Ã X2) and Ï2,3 â P(X2 Ã X3) withidentical marginals µ2 on X2, the Markov concatenation Ï1,2 â¦ Ï2,3 of Ï1,2and Ï1,3 is a multi-coupling in P(X1 ÃX2 ÃX3) given byÏ1,2 â¦ Ï2,3(dx1, dx2, dx3) = Ï1,2(dx1 | x2)µ2(dx2)Ï2,3(dx3 | x2).Here Ï12(dx1 | x2) is the disintegration of Ï12 with respect to µ2 and Ï2,3(dx3 |x2) is the disintegration of Ï2,3 with respect to µ2.63A.7. Entropy-regularized Optimal TransportThe interpretation of the Markov concatenation is as follows: a randomâtrajectoryâ according Ï12â¦Ï23 corresponds to taking the first two steps dis-tributed according to Ï12, then the third step distributed according to âÏ23conditional on the second marginal of Ï12â. The existence of the Markovconcatenation is guaranteed by the disintegration theorem, and Markov con-catenations appear naturally in the time-discretized version of trajectoryinference via Schrodinger bridges [13]. In particular, given Polish spacesX1, . . . , X4 and couplings Ï12, Ï23, Ï34, it holds that Markov concatenationis associative, and so we can unambiguously define the iterated Markov con-catenation Ï12 â¦ Ï23 â¦ Ï34, see [75] Section 3.2.Proposition A.16 (Optimality of anisotropic entropy-regularized optimaltransport). Let the reference SDE be a linear additive noise SDE with driftA and diffusion H (1.3). Given a set of marginals p0, . . . , pNâ1 over times{ti}Nâ1i=0 (not necessarily coming from the reference SDE) for which DKL(pi |Leb) < â for all 1 â¤ i â¤ N â 1, and DKL(pi | (pA,H)i) < â for all 1 â¤i â¤ N â 2, suppose there exists Ïi,i+1, the anisotropic entropic OT solution(4.8) obtained by Sinkhornâs algorithm with marginals µ = pi, ν = pi+1 andtransition kernelKiA,H(x, y) â exp(â12(y â eA(ti+1âti)x)â¤(Î£i)â1(y â eA(ti+1âti)x) (A.25)where Σi =â« ti+1tieA(ti+1âs)HeAâ¤(ti+1âs) ds. Let Ï denote the joint distribu-tion given by the Markov concatenation of couplings (see Definition A.15):Ï = Ï0,1 â¦ . . . â¦ ÏNâ2,Nâ1,Then, Ï minimizes relative entropy to the law of the reference SDE:Ï = argminÏâÎ (p0,...,pNâ1)DKL(Ïâ¥pNA,H),where pNA,H is the law of the reference SDE, discretized over time pointst0, . . . , tNâ1.Proof. For each i = 0, . . . , N â 1, by (1.5), we haveÏi = argminÏâÎ (pi,pi+1)DKL(Ïâ¥KiA,H)Since KiA,H â¼ N (eA(ti+1âti)Xt,Î£i) is the transition kernel of the SDE attime ti [88][Sec 6.1], this implies thatÏi = argminÏâÎ (pi,pi+1)DKL(Ïâ¥(pA,H)i,i+1),64A.7. Entropy-regularized Optimal Transportwhere (pA,H)i,i+1 denotes the joint distribution of the marginals of pA,H atthe two time points ti, ti+1. The details are as follows. Let Ki,i+1(x, y)denote the transition kernel from time ti to time ti+1 for the SDE dXt =AXt+GdWt. In particular, if (pA,H)i is the marginal of this SDE at time ti,then we have thatKi,i+1(x, y)d(pA,H)i(x)dy is equal to the joint distribution(pA,H)i,i+1. HenceDKL(Ï | (pA,H)i,i+1) = DKL(Ï | KA,H(x, y)d(pA,H)i(x)dy).We claim that for any Ï â Î (pi, pi+1),DKL(Ï | KA,H(x, y)d(pA,H)idy)= DKL(Ï | KA,H(x, y)dpi(x)dpi+1(y))+DKL(pi | (pA,H)i) +DKL(pi+1 | Leb).Indeed,DKL(Ï | KA,H(x, y)d(pA,H)idy)=â«log(dÏd(K · (pA,H)i â Leb)(x, y))dÏ=â«log(dÏd(K · pi â dpi+1)(x, y)d(K · pi â dpi+1)d(K · (pA,H)i â Leb)(x, y))dÏ=â«log(dÏd(K · pi â dpi+1)(x, y)dpid(pA,H)i(x)dpi+1dLeb(y))dÏ= DKL(Ï | KA,H(x, y)dpi(x)dpi+1(y))+DKL(pi | (pA,H)i) +DKL(pi+1 | Leb).Note that (pA,H)0 = p0. Hence, under the assumption that for all i =1, . . . , N â 1, DKL(pi+1 | Leb) < â, the minimizers for the following twominimization problems are identical:minÏâÎ (pipi+1)DKL(Ï | (pA,H)i,i+1) andminÏâÎ (pipi+1)DKL(Ï | KA,H(x, y)dpi(x)dpi+1(y)).Since Ï is constructed as a Markov concatenation, the conclusion followsfrom Lemma 3.4 of [75], which in this case tells us that: if pNA,H is theprojection of the law of the SDE onto the set of times {t0, . . . , tNâ1}, then65A.7. Entropy-regularized Optimal Transportfor any N -coupling Ï which is constructed as a Markovian contatenationÏ0,1 â¦ . . . ÏNâ2,Nâ1, and has ith marginal equal to pi, we haveDKL(Ï | pNA,H) =Nâ2âi=0DKL(Ïi | (pA,H)i,i+1)âNâ2âi=1DKL(pi | (pA,H)i).Hence minimizing over each DKL(Ïi,i+1 | (pA,H)i,i+1) (for Ï â Î (pi, pi+1))is equivalent to minimizing over Markovian Ïâs with ith marginal pi.The previous proposition assumes that DKL(pi | Leb) < â for all i Ì¸= 0.The following brief computation shows that this assumption is satisfied inthe case where the diffusion matrix is non-degenerate.We note that Xt â¼ eAtX0+â« t0 eA(tâs)GdWs. As X0 is independent from theBrownian noise Ws, andâ« t0 eA(tâs)GdWs â¼ N (0,Σt), whereΣt =â« t0eA(tâs)HeAâ¤(tâs)ds, it follows thatpt = p0 â N (0,Σt).Thus, we may compute [99, Theorem 9.4.1]DKL(N (0,Σt) | Leb) = 12log(2Ïe)d det(Σt).From the definition of Σt, we note that det(Σt) Ì¸= 0 as long as H = GGâ¤ isfull rank. Then, we use the fact that the DKL(Â· | Leb) is convex, and applyJensenâs inequality for measures, to deduce thatDKL(p0 â N (0,Σt) | Leb) â¤â«DKL(N (x,Σt) | Leb)dp0(x)= DKL(N (0,Σt) | Leb).This shows directly that pt has finite entropy whenever t > 0 as desired.Remark A.17 (Applying EOT for SDE trajectory inference). In the stan-dard ϵ-regularized OT problem (4.7), the cost is the squared Euclidean dis-tance, and K â¼ N (x, ϵ2Id) is an isotropic Gaussian kernel. As noted in[13], this implies that entropy regularized OT can be leveraged for trajectory66A.8. Maximum Likelihood Estimationinference from observed marginals, given a reference SDE. For example, tofind the discretized law on paths Ïâ â Î (µ, ν) satisfying P (Xt = µ,Xt+ât =ν), which minimizes relative entropy to the law of a pure diffusion processdXt = ÏdWt, one should set the entropic regularization ϵ2 = Ï2ât. Thiswould correspond to minimizing the KL divergence to K(x, y) = eââ¥yâxâ¥22Ï2ât â¼N (x, Ï2ât), which is the transition kernel of the reference SDE. Similarly,as done in [8], one can perform trajectory inference given an Ornstein-Uhlenbeck reference SDE dXt = âAXtdt + ÏdWt by approximating thetransition distribution Xt+ât|t via K(x, y) = eâ â¥yâeAâtxâ¥22Ï2ât . This would cor-respond to reweighting the squared Euclidean cost with the drift matrix A,such that c(x, y) = â¥yâeAâtxâ¥/2, and applying standard entropy regularizedOT with ϵ2 = Ï2ât.A.8 Maximum Likelihood EstimationProof of Proposition (4.2). We follow the standard procedure [6, 100] forderiving maximum likelihood estimators. Our likelihood function is givenbyL = Î Nâ2i=0(Î Mâ1j=0 p(X(j)i+1|X(j)i )),where we have denoted Xi = Xti for shorthand. As in [6], we consider thediscretized law Xi+1|Xi â¼ N (Xi +AXiât,Hât), which impliesp(Xi+1|Xi) = 1(2Ïât)d/2det(H)1/2exp(â 12(âX(j)i âAX(j)i ât)â¤(Hât)â1(âX(j)i âAX(j)i ât)),where âX(j)i = X(j)i+1 âX(j)i . Plugging this back into the likelihood expres-sion yields67A.8. Maximum Likelihood EstimationL = 1(2Ïât)dM(Nâ1)2 det(H)M(Nâ1)2. expï£«ï£­â12Nâ2âi=0Mâ1âj=0(âX(j)i âAX(j)i âtï£¶ï£¸â¤ (Hât)â1 (âX(j)i âAX(j)i ât)) .(A.26)log(L) = âM(N â 1)2(d log(2Ïât) + log(det(H)))â 12Nâ2âi=0Mâ1âj=0(âX(j)i âAX(j)i ât)â¤(Hât)â1(âX(j)i âAX(j)i ât)(A.27)We then derive the maximum likelihood estimators through matrix calculus[94]:d log(L)dA= â12Nâ2âi=0Mâ1âj=0ddA(âX(j)i âAX(j)i ât)â¤(Hât)â1(âX(j)i âAX(j)i ât))= â12Nâ2âi=0Mâ1âj=0â2ât ddA((âX(j)i )â¤(Hât)â1AX(j)i)+ât2ddA(X(j)â¤i Aâ¤(Hât)â1AX(j)i)= â12Nâ2âi=0Mâ1âj=0â2ât(Hât)â1(âX(j)i )X(j)â¤i + 2ât2(Hât)â1AX(j)i X(j)â¤i=Nâ2âi=0Mâ1âj=0(H)â1(âX(j)i )X(j)â¤i âât(H)â1AX(j)i X(j)â¤i68A.8. Maximum Likelihood EstimationWe can solve for the MLE linear drift A by setting d log(L)dA = 0:Nâ2âi=0Mâ1âj=0ât(H)â1AX(j)i X(j)â¤i =Nâ2âi=0Mâ1âj=0(H)â1(âX(j)i )X(j)â¤iât(H)â1ANâ2âi=0Mâ1âj=0X(j)i X(j)â¤i = (H)â1Nâ2âi=0Mâ1âj=0âX(j)i X(j)â¤iA =1âtï£«ï£­Nâ2âi=0Mâ1âj=0âX(j)i X(j)â¤iï£¶ï£¸ï£«ï£­Nâ2âi=0Mâ1âj=0X(j)i X(j)â¤iï£¶ï£¸â1We now estimate the diffusion H. For simplicity, we work with Hâ1 =(GGâ¤)â1d log(L)dHâ1=ddHâ1(M(N â 1)2log(det(Hâ1))â 12âtNâ2âi=0Mâ1âj=0(âX(j)i âAX(j)i ât)â¤Hâ1(âX(j)i âAX(j)i ât))=M(N â 1)2H â 12âtNâ2âi=0Mâ1âj=0(âX(j)i âAX(j)i ât)(âX(j)i âAX(j)i ât)â¤We can solve for the MLE additive noise H by setting d log(L)dHâ1 = 0H =1M(N â 1)âtNâ2âi=0Mâ1âj=0(âX(j)i âAX(j)i ât)(âX(j)i âAX(j)i ât)â¤)=1MTNâ2âi=0Mâ1âj=0(âX(j)i âAX(j)i ât)(âX(j)i âAX(j)i ât)â¤).Corollary A.18. Let Xt evolve according to a general additive noise SDE69A.8. Maximum Likelihood Estimation(1.2). Then the log-likelihood function is given byL = 1(2Ïât)dM(Nâ1)2 det(H)M(Nâ1)2. expï£«ï£­âNâ2âi=0Mâ1âj=012(âX(j)i â b(X(j)i )ât)â¤(Hât)â1(âX(j)i â b(X(j)i )ât)ï£¶ï£¸log(L) = âM(N â 1)2(d log(2Ïât) + log(det(H)))â 12Nâ2âi=0Mâ1âj=0(âX(j)i â b(X(j)i )ât)â¤(Hât)â1(âX(j)i â b(X(j)i )ât))Hence, given that the drift function b is parameterized by values α(k)b , themaximum likelihood solution for the drift function b obeys0 = â12Nâ2âi=0Mâ1âj=0ddÎ±(k)b(âX(j)i â b(X(j)i )ât)â¤(Hât)â1(âX(j)i â b(X(j)i )ât))The MLE for the diffusion H admits the closed solutionH =1MTNâ1âi=0Mâ1âj=0(âX(j)i â b(X(j)i )ât)(âX(j)i â b(X(j)i )ât)â¤).Proposition A.19 (Exact MLE estimators for drift and diffusion of 1 di-mensional SDE (1.3) from multiple trajectories). Given M different tra-jectory time series over N different times: {X(j)ti : i â 0, ..., N â 1, j â0, ...,M â 1} sampled from the linear additive noise SDE (1.3), the maxi-mum likelihood solution for linear drift is approximated byAË =1âtlogï£«ï£­âNâ2i=0 âMâ1j=0 X(j)i+1X(j)iâNâ2i=0âMâ1j=0 X(j)2iï£¶ï£¸ (A.28)and the maximum likelihood solution for diffusion is approximated byÏË2 =1M(N â 1)âtNâ2âi=0Mâ1âj=0(X(j)i+1 â eAâtX(j)i )2 (A.29)70A.8. Maximum Likelihood EstimationProof. We proceed as in the proof of Proposition (4.2). The exact log-likelihood for the one dimensional case islog(L) = âM(N â 1)2(log(2Ïât) + log(Ï2))â 12Ï2âtNâ2âi=0Mâ1âj=0(X(j)i+1 â eAâtX(j)i )2To solve for AË, we compute0 =â log(L)âAâNâ2âi=0Mâ1âj=0ââA(â2X(j)i+1X(j)i eAât + e2AâtX(j)2i)e2AâtNâ2âi=0Mâ1âj=0X(j)2i = eAâtNâ2âi=0Mâ1âj=0X(j)i+1X(j)ieAât =âNâ2i=0âMâ1j=0 X(j)i+1X(j)iâNâ2i=0âMâ1j=0 X(j)2iA =1âtlogï£«ï£­âNâ2i=0 âMâ1j=0 X(j)i+1X(j)iâNâ2i=0âMâ1j=0 X(j)2iï£¶ï£¸Similarly, to solve for ÏË2, we compute0 =â log(L)âÏ2â âM(N â 1)2Ï2+12ât1(Ï2)2Nâ2âi=0Mâ1âj=0(X(j)i+1 â eAâtX(j)i )2M(N â 1)2=12âtÏ2Nâ2âi=0Mâ1âj=0(X(j)i+1 â eAâtX(j)i )2Ï2 =1M(N â 1)âtNâ2âi=0Mâ1âj=0(X(j)i+1 â eAâtX(j)i )2Remark A.20. Previous works have predominately focused on the case ofobserving a single long trajectory rather than a collection of short trajecto-ries. Suppose that one observes a set of N observed trajectories, then driftestimation may also be performed by averaging the classical MLE estimatorfor a single trajectory across observations: EN [AËT ] = 1Nâ« T0 Xt dXtâ« T0 X2t dt. Note thatthis is distinct from the MLE estimator that we derived for multiple trajec-tories in (4.11). Indeed, we observe that the latter estimator converges at amuch faster rate than the averaged classical estimator.71A.9. Additional experimentsFigure A.3: The mean squared error of the classical MLE drift estimator(Classical) converges at a slower rate compared to the MLE drift estimator(4.11) derived for multiple trajectories (Trajectory).A.9 Additional experimentsA.9.1 Consistency experiments(a) WOT (b) APPEXFigure A.4: Average mean squared error (MSE) of A and GGâ¤ as the num-ber of observations per marginal increases across 50 random 2d systems.Note this is a log-log plot, and the dotted line in the APPEX plot showsthat the convergence follows a power law relationship with power â0.125.72A.9. Additional experiments(a) WOT (b) APPEXFigure A.5: Average mean squared error (MSE) of A and GGâ¤ as the num-ber of observations per marginal increases across 70 random 12d systems.Note this is a log-log plot, and the dotted line in the APPEX plot showsthat the convergence follows a power law relationship with power -1.73
Thesis/Dissertation
eng
Mathematics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Robeva, Elina
Identifying drift, diffusion, and causal structure from temporal snapshots
Text