Doob's Lagrangian: A Sample-Efficient Variational Approach to Transition Path Sampling

Authors: yuanqi  du
michael  plainer
rob  brekelmans
chenru  duan
frank  noe
carla p gomes
alan  apsuru-guzik
kirill  neklyudov
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 21:02:08

Abstract: Rare event sampling in dynamical systems is a fundamental problem arising in the natural sciences, which poses significant computational challenges due to an exponentially large space of trajectories. For settings where the dynamical system of interest follows a Brownian motion with known drift, the question of conditioning the process to reach a given endpoint or desired rare event is definitively answered by Doob's $h$-transform. However, the naive estimation of this transform is infeasible, as it requires simulating sufficiently many forward trajectories to estimate rare event probabilities. In this work, we propose a variational formulation of Doob's $h$-transform as an optimization problem over trajectories between a given initial point and the desired ending point. To solve this optimization, we propose a simulation-free training objective with a model parameterization that imposes the desired boundary conditions by design. Our approach significantly reduces the search space over trajectories and avoids expensive trajectory simulation and inefficient importance sampling estimators which are required in existing methods. We demonstrate the ability of our method to find feasible transition paths on real-world molecular simulation and protein folding tasks.

Main paper: 
 introduction: 
 Conditioning a stochastic process to obey a particular endpoint distribution, satisfy desired terminal conditions, or observe a rare event is a problem with a long history (Schrödinger, 1932; Doob, 1957) and wide-ranging appli-Figure 1 : Given reference dynamics, transition path sampling seeks to capture the conditional or posterior distribution over paths which reach a terminal set xT ∈ B. However, simulating the reference dynamics (blue) can be wasteful since we rarely obtain paths (orange) which reach (the vicinity) of the terminal set B. This is a major challenge for techniques based on importance sampling or Monte Carlo estimation, even when adding a control term to the reference dynamics. By contrast, our approach optimizes a tractable variational distribution over transition paths with a parameterization satisfying initial and terminal conditions by design. cations from generative modeling (De Bortoli et al., 2021; Chen et al., 2021a; Liu et al., 2022; 2023c; Somnath et al., 2023) to molecular simulation (Anderson, 2007; Wu et al., 2022; Plainer et al., 2023; Holdijk et al., 2024) , drug discovery (Kirmizialtin et al., 2012; 2015; Dickson, 2018) , and materials science (Xi et al., 2013; Selli et al., 2016) .Transition Path Sampling. In this work, we take a particular interest in the problem of transition path sampling (TPS) in computational chemistry (Dellago et al., 2002; Weinan and Vanden-Eijnden, 2010) , which attempts to describe how molecules transition between local energy minima or metastable states under random fluctuations or the influence of external forces. Understanding such transitions has numerous applications for combustion, catalysis, battery, material design, and protein folding (Zeng et al., 2020; Klucznik et al., 2024; Blau et al., 2021; Noé et al., 2009; Escobedo et al., 2009) . While the TPS problem is often framed as finding the 'most probable path' transitioning between states (Dürr and Bach, 1978; Vanden-Eijnden and Heymann, 2008) , we draw explicit connections to the Doob's h-transform and seek to match the full posterior distribution over conditioned processes. Doob's h-Transform. For Brownian motion diffusion processes, conditioning is known to be achieved by Doob's h-transform (Doob, 1957; Särkkä and Solin, 2019) . However, solving this problem amounts to estimating rare event probabilities or matching a complex target distribution. Approaches which involve simulation of trajectories to construct Monte Carlo expectations or importance sampling estimators (Papaspiliopoulos and Roberts, 2012; Schauer et al., 2017; Holdijk et al., 2024) can be extremely inefficient if the target event is rare or endpoint distribution is difficult to match. Recent methods based on score matching (Heng et al., 2021) or nonlinear Feynman-Kac (Chopin et al., 2023) still require simulation within an inner optimization loop.Variational Formulation of Doob's h-Transform. In this work, we propose a variational formulation of Doob's htransform as the solution to an optimization on the space of paths of probability distributions (Thm. 1). We focus on solving for the Doob transform conditioning on a particular terminal point, which is natural in the TPS setting (see Fig. 1 ). Taking inspiration from recent bridge matching methods (Peluchetti, 2021; 2023; Liu et al., 2022; Lipman et al., 2022; Shi et al., 2023; Liu et al., 2023a) , we propose a parameterization with the following attractive features.1. Every Sample Matters. In contrast to most existing approaches, our method is simulation-free, thereby avoiding computationally wasteful simulation methods to estimate rare-event probabilities and inefficient importance or rejection sampling. We thus refer to our approach as being sample-efficient.2. Optimization over Sampling. To approximate the conditioned process, we propose an expressive variational family which is tractable to sample and can be optimized end-to-end using neural networks.3. Problem-Informed Parameterization. Our parameterization enforces the boundary conditions by design, thereby reducing the search space for optimization and efficiently making use of the conditioning information. 
transition path sampling: 
 Consider a forward or reference stochastic process with states x t and transition probability ρ(x t+dt = y | x t = x).Starting from an initial x 0 = A, the density of the path isρ(xT , . . . , x dt | x0 = A) = T -dt t=dt ρ(x t+dt | xt) • ρ(x dt | x0 = A)The problem of rare event sampling aims to condition this reference stochastic process on some event at time T , for example, that the final state belongs to a particular set x T ∈ B. We are interested in sampling from the entire transition path, i.e. the posterior distribution over intermediate statesρ(xT -dt, . . . , xdt | x0 = A, xT ∈ B) = ρ(x T ∈B,x T -dt ...,x dt | x 0 =A) ρ(x T ∈B | x 0 =A).Moving to continuous time, we focus on the transition path sampling problem in the case where the reference process is given by a Brownian motion. In particular, we are motivated by applications in computational chemistry (Dellago et al., 2002; Weinan and Vanden-Eijnden, 2010) , where the reference process is given by molecular dynamics following either overdamped Langevin dynamics,dx t = -(γM ) -1 ∇ x U (x t ) • dt + (γM ) -1 /2 √ 2k B T • dW t , (1)or the second-order Langevin dynamics with spatial coordinates xt and velocities vt ,dx t = vt • dt , (2) dv t = -M -1 ∇ x U (x t ) -γM -1 vt • dt + M -1 /2 √ 2γk B T • dW t .where W t denotes the Wiener process. Note, U is a potential energy, k B T is the Boltzman constant times temperature, M is the mass matrix, and γ is the friction coefficient.Doob's h-transform addresses the question of conditioning a reference Brownian motion to satisfy a terminal condition such as x T ∈ B, thereby providing an avenue to solve the transition path sampling problem described above. Without loss of generality, and to provide a unified treatment of the dynamics in (1)-( 2), we consider the forward or reference stochastic differential equation (SDE),EQUATIONwith drift b t : R N → R N and diffusion matrixΞ t ∈ R N ×N such that G t := 1 2 Ξ t Ξ T t is positive definite. 1We denote the induced path measure as P ref 0:T ∈ P(C([0, T ] → R N )). Remarkably, Doob's h-transform (Doob, 1957; Särkkä and Solin, 2019, Sec. 7.5) shows that conditioning the reference process (3) on x T ∈ B yields another Brownian motion. Proposition 1. [Jamison (1975, Thm. 2)] Let h(x, t) := ρ(x T ∈ B | x t = x) denote the conditional transition density with respect to the reference process in (3). LettingG t := 1 2 Ξ t Ξ T t , the SDE dx t|0,T = bt(x t|0,T ) + 2Gt∇x log h(x t|0,T , t) • dt + Ξt dWt (4)is associated with the following transition probabilitiesEQUATIONfor s < t < T , where we omit the dependence of h(x, t) on B for simplicity of notation.See App. C.1 for proof. The conditioned transition probabilities in (5) allow us to directly construct the transition path in Sec. 2.1. Using Bayes rule, we haveρ(x T -dt , . . . , x dt | x0 = A, xT ∈ B) ρ(x T -dt . . . , x dt | x0 = A) = h(x T -dt , T -dt) h(A, 0)Thus, we can solve the TPS problem by exactly solving for the h-function and simulating the SDE in (4).Theorem 1. The following Lagrangian action minimization has a unique solution which matches Doob's h-transform in Prop. 1, where the optimal q * t|0,T (x) and v * t|0,T (x) = ∇ x log h(x, t) satisfy the PDEs in App. C.1 Prop. 3,EQUATIONEQUATIONEQUATIONFinally, the h-process and marginal density of the conditioned process satisfy a set of forward and backward Kolmogorov equations. These are crucial for deriving our variational objective, but deferred to App. A.1 Prop. 3 for space. 
doob's lagrangian: 
 Consider reference dynamics in the form of either (1) or ( 2), with known drift b t or energy U . We restrict our attention to conditioning on a terminal rare event of reaching a given endpoint x T = B, along with an initial point x 0 = A.We approach finding Doob's h-transform via a least action principle where, in the Thm. 1, we define a Lagrangian action whose minimization yields the optimal q * t|0,T (x) = ρ t|0,T (x) and v * t|0,T (x) = ∇ x log h(x, t) from Prop. 1. This objective will form the basis for our computational approach, with proof of Thm. 1 deferred to App. C.2. We provide additional analysis of our objective in App. C.1.Challenges of Optimizing (6a). We highlight several distinctive features of our problem which inform the development of computational methods in Sec. 3.2. 1. First, (6a) requires that we are able to efficiently sample from the conditioned process in (13) or q t|0,T . This appears challenging due to the nonlinearity of both the reference and variational drifts, b t and v t|0,T .2. For a given q t|0,T , it can be difficult to solve for v t|0,T which satisfies the Fokker-Planck equation in (6b).3. Finally, we would like to strictly enforce the boundary constraints on q t|0,T or Q v 0:T to avoid simulating or wasting computation on trajectories for which x T ̸ = B.Our parameterization of q t|0,T will avoid simulation of the SDE (13) (Challenge 1), provide analytic solutions for v t|0,T satisfying (6b) with given q t|0,T (Challenge 2), and enforce the boundary constraints by design (Challenge 3).We now propose a family of Gaussian (mixture) path parameterizations q t|0,T which overcome the aforementioned computational challenges, while still maintaining expressivity. We present all aspects of our proposed method in the context of the first-order dynamics and defer extensions to mixture paths and second-order dynamics to App. B.Tractable Drift v t|0,T for Variational Doob Objective. Consider modifying the Fokker-Planck constraint in (6b), to absorb all drift terms into a single vector field u t|0,T ,∂q t|0,T (x) ∂t = -∇x, q t|0,T (x) u t|0,T (x) + ij (Gt)ij ∂ 2 ∂xi∂xj q t|0,T (x). (7)To address Challenge 2, we restrict attention to variational families of q t|0,T ∈ Q where it is analytically tractable to calculate a vector field u (q) t|0,T which satisfies (7). We first consider the family of Gaussian paths Q G , in similar fashion to flow matching methods (Lipman et al., 2022; Tong et al., 2023; Liu et al., 2023a) , with proof in App. B. Proposition 2. For the family of endpoint-conditioned marginals q t|0,T (x) = N (x | µ t|0,T , Σ t|0,T ),u (q) t|0,T (x) := ∂µt |0,T ∂t + 1 2 ∂Σt |0,T ∂t Σ -1 t|0,T -G t Σ -1 t|0,Tx -µ t|0,T (8)satisfies the Fokker-Planck equation (7) for q t|0,T and diffusion coefficientsG t = 1 2 Ξ t Ξ T t .Given u (q) t|0,T corresponding to q t|0,T , we can simply solve for the v t|0,T satisfying the Fokker-Planck euqation in (6b) in our variational Doob objective (Thm. 1). Since G t was assumed to be invertible and the base drift b t is known,EQUATIONWe may now evaluate terms involving v t|0,T in our Lagrangian objective in (6) using (9) directly, without solving an inner minimization over v t|0,T (addressing Challenge 2).Optimization over q t|0,T satisfying (6c). Given the ability to evaluate v (q) t|0,T for a given q t|0,T ∈ Q G as above, our variational Doob objective in (6a) reduces to a single optimization over the marginals q t|0,T of a conditioned process which satisfies the boundary conditions (6c).We consider parameterizing the mean µ t|0,T and covariance Σ t|0,T of our Gaussian path q t|0,T using a neural network. For simplicity, we consider a diagonal parameterizationΣ t|0,T = diag({σ 2 t|0,T ,d } D d=1 ). We parameterize a neural network NNET θ : [0, T ] × R D × R D → R D × R D which inputs (t, A, B). Its output is used to construct x t|0,T = µ (θ) t|0,T + Σ (θ) t|0,T ϵ, where ϵ ∼ N (0, ID).(10)µ (θ) t|0,T = (1 -t)A + t B + t(1 -t)NNET θ (t, A, B) [:D] Σ (θ) t|0,T = t(1 -t)diag NNET θ (t, A, B) [D:] + σ 2 min I.Crucially, our Gaussian parameterization addresses Challenge 1, in that we can easily draw samples x t|0,T ∼ q t|0,T from our variational conditioned process (6b) without simulating the corresponding SDE with nonlinear drift (13). Further, t(1 -t) coefficients in (10) ensure that the (smoothed) boundary conditions are satisfied by design (Challenge 3). We add σ 2 min to ensure invertibilty of Σ t|0,T (see ( 8)) as t → 0 or t → T , but preserve q 0 (x 0 ) = N (x 0 | A, σ 2 min I D ) ≈ δ(x 0 -A) and q T (x T ) = N (x T | B, σ 2 min I D ) ≈ δ(x T -B). Reparameterization Gradients. Since we have now shown that our parameterization satisfies the constraints (6b)-( 6c) by design, we can finally optimize our variational Doob objective with respect to q t|0,T ∈ Q G using the reparameterization trick (Kingma and Welling, 2013; Rezende et al., 2014) . In particular, for the expectation at each t in (6a),∇ θ E q (θ) t|0,T (x) v (q,θ) t|0,T (x), Gt v (q,θ) t|0,T (x) (11) = E N (ϵ|0,ID) ∇ θ v (q,θ) t|0,T (g(t, ϵ; θ) , Gt v (q,θ) t|0,T (g(t, ϵ; θ)where x = g(t, ϵ; θ) is the parameterized map in (10) and vEQUATIONt|0,T in ( 8)-( 9). In practice, we approximate gradients using a single sample of ϵ at uniformly sampled discrete time points 0 ≤ t ≤ T which represent physical time (e.g., femtoseconds). 
experiments: 
 We investigate the capabilities of our approach across a variety of different settings. We first illustrate features of our method on toy examples before continuing to real-world molecular systems, including a commonly-used benchmark system, alanine dipeptide, and a small protein, Chignolin. The code behind our method is available at the following link. Before diving into the experiments, we introduce the evaluation procedure and baseline methods.Evaluation metrics. In our evaluation, we emphasize two key quantities: accuracy and efficiency. Efficiency is evaluated by the number of calls to the potential energy function, which requires extensive computation and dominates the runtime of larger molecules. For accuracy, we evaluate the log-likelihood of each sampled path and the maximum energy point (saddle point/transition state) along each sampled path. A good method samples many probable paths (i.e., high log-likelihood) and an accurate transition state (i.e., small maximum energy). See App. E for further details.Baselines. We compare our approach against the Markov Chain Monte Carlo (MCMC)-based two-way shooting method with uniform point selection with variable or fixed length trajectories. We found that two-way shooting produced the most diverse path ensembles among possible baselines, although the acceptance probability can be relatively low for systems dominated by diffusive dynamics (Brotzakis and Bolhuis, 2016) and might be improved by learning shooting point selection. This baseline gives theoretical guarantees about the ensemble and thus can be considered as a proxy for the ground truth.Müller-Brown Potential. The Müller-Brown potential is a popular benchmark to study transition path sampling between metastable states. It consists of three local minima, and we aim to sample transition paths connecting state A and state B with a circular state definition. In Fig. 2 , we visualize the potential and the sampled paths and can see that the same ensemble is sampled for both our method and two-way shooting. Our method exhibits a slightly reduced variance for unlikely transitions. In Table 1 , we can observe that MCMC-based methods require many potential evaluations to achieve a good result, which comes from the low acceptance rate (especially when fixing the lengths of trajectories). Our method requires fewer energy evaluations (1 million vs. 1 billion) while finding paths with similar energy and likelihood. Note, the likelihood for variable approaches has been omitted, as it is governed by the number of steps in the trajectory and cannot be compared directly.Gaussian Mixture. We further consider a potential in which the states are separated by a symmetric high-energy barrier that allows for two distinct reaction channels. In Fig. 3 , we observe that a single Gaussian path cannot model a system with multiple modes of transition paths. Nevertheless, this issue can be resolved using a mixture of Gaussian paths, with slightly increased computational cost.Experiment Setup. We evaluate our methods on realworld high-dimensional molecular systems governed by the second-order dynamics (2): alanine dipeptide and Chignolin. Alanine dipeptide is a well-studied system of 22 atoms (66 total degrees of freedom), where the molecule can be described by two collective variables (CV): the dihedral angles ϕ, ψ. Chignolin is a larger system consisting of 10 residues with 138 atoms (414 total degrees of freedom) that cannot be summarized as easily. We use an AMBER14 forcefield (Maier et al., 2015) implemented in OpenMM (Eastman et al., 2017) but use DMFF (Wang et al., 2023) to backpropagate through the energy evaluations.Alanine Dipeptide. In Table 2 , we report results for four variants of our models, which either predict Cartesian coordinates or internal coordinates in the form of bond lengths and dihedral angles (compare App. E), either with or without Gaussian mixture. For our method, operating in internal coordinates takes more iterations to converge but generates better results compared to Cartesian coordinates, where internal coordinates have nicely distributed input and our network does need not learn equivariances (Du et al., 2022) . Similarly, Gaussian mixture paths perform slightly better than a single Gaussian path due to the extra expressiveness.We note that paths sampled with Gaussian mixture exhibit a larger variance in max energy as they represent multiple reaction channels.We find that prior-informed definitions of the desired initial and target states (i.e., CV) are necessary for MCMC to work efficiently with fixed-length trajectories. Finding these CVs in practice is challenging and only possible in this instance because the molecule is small and well-studied. For the larger system size in Table 2 , it becomes intractable to use MCMC for reaching precise states A, B ('exact') instead of larger regions ('relaxed'), or for computing fixed-length trajectories. Variable length MCMC with relaxed endpoint conditions and fixed-length MCMC with CV perform well on this task, but our method is competitive using fewer evaluations and more strict boundary conditions.Chignolin. The folding dynamics of Chignolin already pose a challenge and have not yet been well-studied compared to alanine dipeptide. We illustrate the qualitative experimental results for this system in Fig. 4 . Operating in Cartesian space, our model samples a feasible transition within 12.8M potential energy evaluation calls and a transition with a duration of T = 1ps, which is faster compared to T = 0.6µs in Lindorff-Larsen et al. (2011). 
conclusion, limitations and future work: 
 In this paper, we propose an efficient computational framework for transition path sampling with Brownian dynamics. We formulate the transition path sampling problem by using Doob's h-transform to condition a reference stochastic process, and propose a variational formulation for efficient optimization. Specifically, we propose a simulation-free training objective and model parameterization that imposes boundary conditions as hard constraints. We compare our methods with MCMC-based baselines and show comparable accuracy with lower computational costs on both synthetic datasets and real-world molecular systems. Finally, our method might be improved or extended by (1) accounting for conditioning on a set of terminal events, (2) amortizing over many state pairs or systems and finally learning an unconditioned process, and (3) accommodating variable length paths.  
 
 Appendix: EQUATION)Reparameterizing (12b) in terms of s(x, t) := log h(x, t), we can also writeEQUATION)We provide proof of this proposition in App. C.1. We recall the simple example of the Brownian bridge, although our methods will handle the case of nonlinear reference processes as in ( 1)-( 2).Example 1 (Brownian Bridge). Consider a reference process dx t = dW t . Conditioning on x T = B for a particular value of B, the Doob h-transform amounts to a Gaussian transition kernel h(x, t):= ρ(x T = B | x t = x) = N (x T | x t , T -t).Plugging the log gradient ∇ x log h(x, t) into (4) yields the conditioned SDE known as the Brownian bridge,dx t|0,T = (x T -x t|0,T )/(T -t)dt + dW t .A.2. Analysis of the Variational Objective in Thm. 1Unconstrained Dual Objective. Introducing Lagrange multipliers to enforce the constraints in (6b)-(6c) and eliminating v t|0,T , we obtain an alternative, unconstrained version of (6a).Corollary 1. The Lagrangian objective in Thm. 1 which solves Doob's h-transform is equivalent toS = min q t|0,T max s s(B, 1) -s(A, 0) - 1 0dt dx q t|0,T ∂s ∂t + ∇s, Gt∇s + ∇s, bt + ∇, Gt∇s if q t|0,T satisfies (6c). Note v t|0,T (x) = ∇ x s(x, t), with s * (x, t) = log h(x, t) at optimality.2This objective is similar to the objectives optimized by Action Matching methods (Neklyudov et al., 2023; 2024) . Notably, the objective in Cor. 1 is expressed directly in terms of the (log) of the h-function for fixed conditioning information x T = B.We also note that the Hamilton Jacobi-style quantity, whose expectation appears in the final term, is zero at optimality in (12c) of Prop. 3.Path Measure Perspective. We next interpret our variational objective in Thm. 1 as minimizing a KL divergence over path measures. Let P ref 0:T denote the law of the reference SDE in (3) with fixedP ref 0 = δ(x 0 -A). Let Q v 0:T denote the law of a controlled process similar to (4), but with a variational v t|0,T in place of ∇ x log h,EQUATION)Note that the density q t|0,T of Q v 0:T satisfies the Fokker-Planck equation in (6b) (Särkkä and Solin, 2019, Sec. 5 .2) Using the Girsanov Theorem, the objective in (6a) can then be viewed as a KL divergence minimization over path measures Q v 0:T which satisfy the boundary constraints. Corollary 2. The Lagrangian objective in Thm. 1 is equivalent to the following optimization of Q v 0:T S := minEQUATIONwhere the minimizing argument recovers the path measure P * 0:T associated with the SDE in (4).Our Lagrangian action minimization thus solves a Schrödinger Bridge (SB) problem (Schrödinger, 1932; Léonard, 2014) with Dirac delta functions as the endpoint measures. Our objective in (6a) particularly resembles optimal control formulations of SB (Chen et al., 2016; 2021b) . While it is well-known that the Doob h-transforms (and large deviations more generally) play a role in the solution to SB problems (e.g. Jamison (1975) ; Léonard (2014)), our interest in the transition path sampling problem leads to specific computational decisions below. See App. D for further discussion.We begin by proving Prop. 2 in the main text, before discussing how our algorithm extends to second order dynamics (App. B.1) and mixtures of Gaussians (App. B.2). Proposition. 3. For the family of endpoint-conditioned marginals q t|0,T (x) = N (x | µ t|0,T , Σ t|0,T ),EQUATIONsatisfies the Fokker-Planck equation (7) for q t|0,T and diffusion coefficientsG t = 1 2 Ξ t Ξ T t .Proof. Consider the following identities for the Gaussian family of marginals q t (x) = N (x|µ t , Σ t ), where we omit conditioning q t ← q t|0,T for simplicity of notation,EQUATIONEQUATIONEQUATIONWe begin by solving for a vector field u o t (x) that satisfies the continuity equation (where u o t denotes the drift of an ODE)EQUATION)The vector field satisfying this equation isEQUATIONwhich we can confirm using the identities in (16). Indeed, for the terms on the RHS of Eq. ( 17),-∇ x , u o t = - 1 2 tr Σ -1 t ∂Σ t ∂t , -∇ x log q t , u o t = Σ -1 t (x -µ t ), ∂µ t ∂t + 1 2 (x -µ t ) T Σ -1 t ∂Σ t ∂t Σ -1 t (x -µ t ) .Putting these terms and the time derivative from (16c) into Eq. ( 17) we conclude the proof.However, we are eventually interested in finding the formula for the drift u t that satisfies the Fokker-Planck equation in (7).That is, to describe the same evolution of density ∂qt(x) ∂t , the relationship between u t and u o t is as follows∂q t (x) ∂t = -∇ x , q t u o t = -∇ x , q t u t + ∇ x , G t ∇ x q t = -∇ x , q t u t + ∇ x , G t q t ∇ x log q t = -∇ x , q t (u t -G t ∇ x log q t ) u o tFinally, we use the identities in ( 16) to obtainu t = u o t + G t ∇ x log q t = ∂µ t ∂t + 1 2 ∂Σ t ∂t Σ -1 t (x -µ t ) -G t Σ -1 t (x -µ t ) =⇒ u t = ∂µ t ∂t + 1 2 ∂Σ t ∂t Σ -1 t -G t Σ -1 t (x -µ t )To handle the case of the second-order dynamics in (2), we can adapt our recipe from the previous section with minimal modifications by extending the state space x ∈ R D to include velocities v, with x = (x, v) ∈ R 2D . However, note that the dynamics in (2) are no longer stochastic in the spatial coordinates x. To ensure invertibility of G t and existence of the h-transform, we add a small nonzero diffusion coefficient in the coordinate space x, so that the reference process in Eq. ( 3) is given byEQUATION)All steps in our algorithm proceed in similar fashion to Sec. 3.2. We now parameterize q t|0,T (x, v) using NNET θ : [0, T ]×R 2D ×R 2D → R 2D ×R 2D , which outputs mean perturbations and per-dimension variances to calculate µ x t|0,T , µ v t|0,T and Σ x t|0,T , Σ v t|0,T and sample (x, v), as in ( 10). Note that we parameterize Σ x t|0,T , Σ v t|0,T separately, matching the block diagonal form of ( 19). We calculate v(q) t|0,T (x, v) := [v x(q) t|0,T , v v(q) t|0,T ] from u (q) t|0,T (x, v) = [u x(q) t|0,T , u v(q)t|0,T ] as in ( 8)-( 9), withG -1 t = ( 1 2 Ξ t Ξ T t ) -1given by ( 19). The Lagrangian objective in ( 6) minimizes the norm of the concatenated vector v (q) t|0,T (x, v), which depends on the reference drift b t (x, v) in ( 19).Note that the true Doob h-transform may not yield marginals which follow the unimodal Gaussian distributions in the previous sections. To increase the expressivity of our variational family of conditioned processes, we consider extending our parameterization to mixtures of Gaussians, q t|0,T ∈ Q K MoG . For a given set of K mixture weights w k and component Gaussian paths q k t|0,T , the following identity allows us to recover the drift u (q) t|0,T of the corresponding mixture distribution q t|0,T .Proposition 4. Given a set of processes q k t|0,T (x) and mixtures weights w k , the vector field satisfying the Fokker-Planck equation in (7) for the mixture q t|0,T (x) = k w k q k t|0,T (x) is given byEQUATIONwhere u (q,k) t|0,T (x) satisfies the Fokker-Planck equation in (7) for q k t|0,T (x). This identity holds for both first order dynamics in spatial coordinates only or second-order dynamics in x = (x, v).Proof. See Peluchetti (2023) Theorem 1 and its proof in their App. C. Finally, we can calculate v (q) t|0,T (x) by comparing u (q) t|0,T (x) for the mixture of Gaussian path q t|0,T ∈ Q K MoG to the reference drift b t (x) as in ( 9), and proceed to minimize its norm as in (6). We use Gumbel softmax reparamerization gradients (Maddison et al., 2016; Jang et al., 2017) to optimize the mixture weights {w k } K k=1 alongside the neural network parameters {θ k } K k=1 for each Gaussian component {µ(θ k ) t|0,T , Σ (θ k )t|0,T } K k=1 and either first-or second-order dynamics.C.1. Proofs from Sec. 2.2 (Doob's h-Transform Background)Proposition. 2. [Jamison (1975, Thm. 2)] Let h(x, t) := ρ(x T ∈ B | x t = x) denote the conditional transition density with respect to the reference process in (3).Letting G t := 1 2 Ξ t Ξ T t , the SDE dx t|0,T = b t (x t|0,T ) + 2G t ∇ x log h(x t|0,T , t) • dt + Ξ t dW t (21)is associated with the following transition probabilitiesEQUATIONfor s < t < T , where we omit the dependence of h(x, t) on B for simplicity of notation.Proof. See Jamison (1975) for a simple proof based on Ito's Lemma, assuming smoothness and strict positivity of h.Proposition 3. The following PDEs are obeyed by (a) the marginal density of the conditioned process ρ t|0,T (x) = ρ(x t = x | x 0 = A, x T ∈ B) and (b) the h-function h(x, t) (which implicitly depends on B),EQUATION)Reparameterizing (12b) in terms of s(x, t) := log h(x, t), we can also writeEQUATIONProof. Let p(x t+s = y | x t = x) denote the transition probability of a reference diffusion processEQUATIONwhere (G t ) ij = 1 2 Ξ t+s Ξ T t+s . Now we condition the process on the end-point value x T ∈ B, and we get another kernel, i.e.EQUATIONWe let h(x, t) = p(x T ∈ B | x t = x) denote the conditional probability over the desired endpoint condition given x t = x.According to laws of conditional probability, we can describe how h(x, t) changes in time using the unconditioned transition probabilityEQUATIONwe take the derivative ∂ ∂s on both sides, and we getEQUATIONUsing the FP equation for the transition probability and integrating by parts, we have0 = dy p(xt+s = y | xt = x)   ∂h(y, t + s) ∂s + ∇yh(y, t + s), bt(y) + ij (Gt) ij ∂ 2 ∂y i ∂y j h(y, t + s)   .Note that this holds ∀x, hence, we have∂h(y, t + s) ∂s + ∇ y h(y, t + s), b t+s (y) + ij (G t ) ij ∂ 2 ∂y i ∂y j h(y, t + s) = 0 ,without any loss of generality we can set t = 0 ∂h(y, s) ∂sEQUATION)as desired to prove the optimality condition in (12b).To prove (12a), denote p(y, s) = p(x s = y | x 0 = x) and differentiate p(EQUATION(y) -p(y, s) ∇ y , 2D∇ y h(y, s) h(x, 0) ± ∇ y p(y, s), 2D∇ y h(y, s) h(x, 0) + ij (G t ) ij ∂ 2 ∂y i ∂y j h(y, s) h(x, 0) p(y, s) ,Note that h(x, 0) can be pulled outside the differential operator because it is a function of x. The PDE for the new kernel p(y, s | B) = p(x s = y | x 0 = x, x T ∈ B) (conditioned on the end-point) becomesEQUATIONwhich matches the desired PDE in (12a) thereby proving the first two parts of Prop. 3.Finally, to show (12c), we index time using t in Eq. ( 27) and change variables h(x, t) = e s(x,t) ,∂e s(x,t) ∂t + ∇ x e s(x,t) , b t (x) + ij (G t ) ij ∂ 2 ∂x i ∂x j e s(x,t) = 0 . e s(x,t) ∂s(x, t) ∂t + e s(x,t) ∇ x s(x, t), b t (x) + ∇, G t ∇e s(x,t) = 0Next, we simplify ∇, G t ∇e s(x,t) = ∇, G t e s(x,t) ∇s(x, t) = ∇e s(x,t) , G t ∇s(x, t) + e s(x,t) ∇, G t ∇s(x, t) = e s(x,t) ∇s(x, t), G t ∇s(x, t) + e s(x,t) ∇, G t ∇s(x, t) to finally writee s(x,t)   ∂s(x, t) ∂t + ∇ x s(x, t), b t (x) + ∇s(x, t), G t ∇s(x, t) + ij (G t ) ij ∂ 2 ∂x i ∂x j s(x, t)   = 0which demonstrates (12c) since the inner term must be zero.We begin by proving Cor. 1, whose proof actually contains the initial steps needed to prove our main theorem Thm. 1. In both proofs, we omit conditioning notation q t ← q t|0,T for simplicity and assume q t (x)s t (x) → 0 vanishes at the boundary x → ±∞, which is used when integrating by parts in x.Corollary 1. The Lagrangian objective in Thm. 1 which solves Doob's h-transform is equivalent toS = min q t|0,T max s s(B, 1) -s(A, 0) - 1 0dt dx q t|0,T ∂s ∂t + ∇s, Gt∇s + ∇s, bt + ∇, Gt∇s if q t|0,T satisfies (6c). Note v t|0,T (x) = ∇ x s(x, t), with s * (x, t) = log h(x, t) at optimality.3 Proof. Consider the following action functionalS = min qt ,vt dt dx q t (x) v t (x), G t v t (x) , s.t. ∂q t (x) ∂t = -∇ x , q t (x)(b t (x) + 2G t v t (x)) + ij (G t ) ij ∂ 2 ∂x i ∂x j q t (x) , q 0 (x) = δ(x -A), q 1 (x) = δ(x -B) .The Lagrangian of this optimization problem isL = 1 0 dt dx   q t v t , G t v t + s t   ∂q t ∂t + ∇, q t (b t + 2G t v t ) - ij (G t ) ij ∂ 2 ∂x i ∂x j q t     ,where s t is the dual variable and we omit the optimization arguments, with S = min qt ,vt max st L. Swapping the order of optimizations under strong duality, we take the variation with respect to v t in an arbitrary direction h t . Using G t = G T t , we obtainEQUATIONSubstituting into the above, we haveEQUATIONIntegrating by parts in t and in x, assuming that q t (x)s t (x) → 0 as x → ±∞, yieldsEQUATIONwhere in the second line, we integrate by parts in x again. Enforcing q 1 (x) = δ(x -B) and q 0 (x) = δ(x -A) and recalling S = min qt max st L after eliminating v t , we recover the optimization in the statement of the corollary.Theorem. 1. The following Lagrangian action functional has a unique solution which matches the Doob h-transform in Prop. 3,EQUATIONEQUATIONNamely, the optimal q * t|0,T (x) obeys (12a) and the optimal v * t|0,T (x) = ∇ x log h(x, t) = ∇ x s(x, t) follows (12b) or (12c).Proof. The proof proceeds from (30) above, S = minEQUATION)We first show that the optimality condition with respect to s t yields the Fokker-Planck equation for q t in Prop. 3 (12a), before deriving the PDE in (12b) as the optimality condition with respect to q t .Optimality Condition for (32) recovers Prop. 3 (12a): The variation with respect to s t of ( 33 EQUATION)where in (i) we integrate by parts x.We are now ready to set the variation of (33) with respect to s t (in an arbitrary direction h t ) equal to zero. Using (34), we haveEQUATIONwhich matches the desired optimality condition for the conditioned process in Prop. 3 (12a) .Optimality Condition for (32) recovers Prop. 3 (12b): Starting again from (33), we take the variation with respect to q t . First, we repeat identical steps (integrate by parts in both x and t) to reach (31),L = dx q 1 s 1 -dx q 0 s 0 - 1 0 dt dx q t ∂s t ∂t + ∇s t , G t ∇s t + ∇s t , b t + ∇, G t ∇s twhere it is now clear that taking the variation with respect to q t and setting equal to zero yieldsδL δq t [h t ] = 0 = ∂s t ∂t + ∇s t , G t ∇s t + ∇s t , b t + ∇, G t ∇s t (36)which is the desired PDE for s(x, t) = log h(x, t) in (12c). To obtain (12b), we note an identity used to simplify the last termij (Gt) ij ∂ 2 ∂x i ∂x j log ht = ∇, Gt∇ log ht = ∇, 1 ht Gt∇ht = - 1 h 2 t ∇ht, Gt∇ht + 1 ht ∇, Gt∇ht .Now, substituting s(x, t) = log h(x, t) into Eq. ( 36), we obtainEQUATIONwhich matches (12b) as desired.The last equation defines the backward Kolmogorov equation for the diffusion process with the drift b t (x) and covariance matrix G t , i.e. the function h t (x) defines the conditional density h t (x) = p(x T ∈ B ′ | x t = x) for some set B ′ , which agrees with the forward process with the same drift and covariance. The boundary condition q T (x) = δ(x -B) together with the backward equation define the unique solution to this PDE. Since the PDEs and the boundary conditions are the same as in Doob's h-transform, we have ht (x) = p(x T = B | x t = x).Corollary 2. The Lagrangian objective in Thm. 1 is equivalent to the following optimization of Q v 0:T S := minEQUATIONwhere the minimizing argument recovers the path measure P * 0:T associated with the SDE in (4).Proof. We use the Girsanov theorem (Särkkä and Solin, 2019, Sec. 7 .3) to calculate the KL divergence between the following two Brownian diffusions with fixed initial condition x 0 = A,EQUATIONEQUATIONIn particular, noting the difference of drifts is bt (x t ) + 2G t v t|0,T (x t ) -b t (x t ) = 2G t v t|0,T (x t ), the likelihood ratio is given bydQ v 0:T dP ref 0:T = q t|0,T (x 0 , ...x T ) ρ(x 0 , ...x T ) = exp - 1 2 T 0 2G t v t|0,T (x t ), (G t ) -1 2G t v t|0,T (x t ) dt (40) -2 G t v t|0,T (x t ) T G -1 t dW tWe finally calculate the KL divergence, noting that, after taking the log, the expectation of the integral ( • )dW t in the final term vanishes,EQUATIONwhich matches (6a) up to a constant factor of 2 does not change the optimum. We finally compare to the constraints in Thm. 1. First, it is clear that the diffusion in (39) satisfies the Fokker-Planck equation in (6b) (Särkkä and Solin, 2019, Sec. 5.2) . We respect (6c) by optimizing over endpoint-constrained path measures, which yieldsEQUATIONas desired.(Aligned) Schrödinger Bridge Matching Methods. Many existing 'bridge matching' approaches (Shi et al., 2023; Peluchetti, 2021; 2023; Liu et al., 2022; Lipman et al., 2022; Liu et al., 2023b) for SB and generative modeling rely on convenient properties of Brownian bridges and would require calculating h-transforms to simulate bridges for general reference processes. Our conditional Gaussian path parameterization is similar to Liu et al. (2023a) ; Neklyudov et al. (2024) , where analytic bridges are not available for SB problems with nonlinear reference drift or general costs. Somnath et al. (2023) ; Liu et al. (2023b) attempt to solve the SB problem given access to aligned data x 0 , x T ∼ q data 0,T assumed to be drawn from an optimal coupling. While the method in Somnath et al. (2023) involves approximating an h-transform, their goal is to obtain an unconditioned vector field v t to simulate a Markov process. However, De Bortoli et al.(2023) use Doob's h-transform to argue the learned Markov process will not preserve the empirical coupling unless q data 0,T is the optimal coupling for the SB problem, and show that an 'augmented' v 0,t which conditions on x 0 can correct this issue.After training on a dataset of x 0 , x T ∼ q data 0,T pairs using our method, we could consider using an (augmented) bridge matching objective (Shi et al., 2023; De Bortoli et al., 2023) to distill our learned v (q) t|0,T into a vector field v t or v 0,t which does not condition on the endpoint. Our use of a Gaussian path parameterization with samples from a fixed endpoint coupling and no Markovization step corresponds to a simplified version of the conditional optimal control step in Liu et al. (2023a) .Transition Path Sampling. We refer to the surveys of Dellago et al. (2002) ; Weinan and Vanden-Eijnden (2010) ; Bolhuis and Swenson (2021) for an overview of the TPS problem. Least action principles for TPS have a long history, building upon the Freidlin-Wentzell (Freidlin and Wentzell, 1998) and Onsager-Machlup (Onsager and Machlup, 1953; Dürr and Bach, 1978) Lagrangian functionals in the zero-noise limit and finite-noise cases. In particular, the Onsager-Machlup functional relates maximum a posteriori estimators or 'most probable (conditioned) paths' to the minimizers of an action functional similar to Thm. 1, where example algorithms include (Vanden-Eijnden and Heymann, 2008; Sheppard et al., 2008) . By contrast, our approach targets the entire posterior over transition paths using an expressive variational family. While Lu et al. (2017) provide analysis for the Gaussian family, we draw connections with Doob's h-transform and extend to mixtures of Gaussians.Shooting methods are among the most popular for sampling the posterior of transition paths. From a path that satisfies the boundary conditions (obtained, e.g., using high-temperature simulations), shooting picks points and directions in which to propose alterations, then simulates new trajectories and accepts or rejects using Metropolis-Hastings (MH) (Juraszek and Bolhuis, 2008; Borrero and Dellago, 2016; Jung et al., 2017; Falkner et al., 2023; Jung et al., 2023) . While the MCMC corrections yield theoretical guarantees, shooting methods involve expensive molecular dynamics (MD) simulations and need to balance high rejection rates with large changes in trajectories. One-way shooting methods sample paths efficiently but yield highly correlated samples. Two-way shooting methods, which we compare against in Sec. 4, are more expensive but typically sample diverse paths faster. Recent machine learning approaches such as Plainer et al. (2023) ; Lelièvre et al. (2023) aim to reduce the need for MD. Holdijk et al. (2024) propose a stochastic optimal control method that simulates (13) with a learned drift, but can be inefficient unless the terminal state is sampled frequently.The main dilemma of molecular dynamics comes from the accuracy and efficiency trade-off-accurate simulation requires solving the Schrödinger equation which is computationally intractable for large systems, while efficient simulation relies on empirical force fields which is inaccurate. Recently, there has been a surge of work in applying machine learning approaches to accelerate molecular simulation. One successful paradigm is machine learning force field (MLFF) which leverages the transferability and efficiency of machine learning methods to fit force/energy prediction models on quantum mechanical data and transfer across different atomic systems (Smith et al., 2017; Wang et al., 2018) . More recently, increasing attention has been focused on building atomic foundation models to encompass all types of molecular structures (Batatia et al., 2023; Shoghi et al., 2023; Zhang et al., 2022) .Sampling is a classical problem in molecular dynamics to draw samples from the Boltzmann distribution of molecular systems. Classical methods mainly rely on Markov chain Monte Carlo (MCMC) or MD which requires long mixing time for multimodal distributions with high energy barriers (Rotskoff, 2024) . Generative models in machine learning demonstrate promises in alleviating this problem by learning to draw independent samples from the Boltzmann distribution of molecular systems (known as Boltzmann generator) (Noé et al., 2019) . Numerous methods have been developed to utilize generative models as a proposal distribution for escaping local minima in running MCMC methods (Gabrié et al., 2022) . However, one critical issue is that generative models rely on training from samples. Although recent advances have been developed to learn from unnormalized density (i.e., energy) function, the training inefficiency limits their applicability to solve high-dimensional molecular dynamics problems. To circumvent the curse of dimensionality for the sampling problem, another branch of work study to learn coarse-grained representation with neural networks (Sidky et al., 2020) . For broader literature of applying machine learning to enhanced sampling, we refer the reader to Mehdi et al. (2024) .See (19) in App. B.1 to write (2) in the form of (3).Again, note that we omit the dependence of s(x, t) and h(x, t) on the conditioning information B.Again, note that we omit the dependence of s(x, t) and h(x, t) on the conditioning information B. 


Decision: This is a strong manuscript that proposes a novel variational technique for Monte Carlo sampling of random dynamical systems which is justified theoretically as well as experimentally. The reviewers are unanimous in its technical merit both conceptually and in practice. After the authors make a few improvements to the manuscript in terms of presentation and references, it will be an excellent addition to the NeurIPS conference program.

Reviews: review: paper_summary: paper_summary: This paper proposes a variational formulation of the Doob's h-transform to change the problem from the expensive simulations of trajectories to an optimization problem over possible trajectories from given initial point and end point. The model parameterization introduced imposes the desired boundary conditions by design and uses a mixture of Gaussian as the variational family. The variational formulation offers an alternative to expensive simulations and MCMC which is hard to scale. The proposed framework is applied on a simulated case study and real world problem in material science/computational chemistry of transition path sampling  which is a study of molecule transitions betweeb local energy minima and metastable states under random fluctuations. The results are along expected lines where it is possible to scale the solution on computationally intractable problems with MCMC and other sampling based methods.
strength_weakness: strength_weakness: Minor things
- Increase font sixe for Table1, bold the important results.
- I would have liked to see a bit more dicussion on the dimensionality, what are the typical values of D ..
questions: questions: 1. I am curious,  similar to Figure 3, the authors can show an experiment with familiy of dtistributions other than Gaussian. How does the expressivity then come to play ? 
2. as for any SDE problem solution, what is the effect of discretization values and parameters in practice as given in line 184-185.
3. How big is the constraint of fixed length(T) on transition paths in practice as solved in this method ? 
4. What is the behaviour of the problem when there are multiple basins of attraction ?
limitations: limitations: The authors have addressed the challenges quite well and what the future research directions could be taken.
paper_summary: paper_summary: The submitted manuscript presents a variational formulation of Doob’s h-transform, leading to a novel (simulation-free) computational approach for rare event sampling in transition paths. The task of interest involves conditioning a dynamical system driven by Brownian motion with a known drift term to reach a given endpoint. In theory, this terminal condition can be addressed using Doob’s h-transform, resulting in an associated SDE for the conditional dynamic. However, practical implementation requires knowledge of the h-transform. The authors propose a variational problem formulation that provides the necessary information about the h-transform as its solution. Solving this variational problem results in a computational approach for simulating the conditional dynamic. Compared to existing methods, it is claimed that the proposed approach avoids importance sampling estimators and expensive trajectory simulations. The method is tested on both synthetic and real datasets.
strength_weakness: strength_weakness: - The proposed approach in Section 3.2.1 introduces certain issues. The connection to Doob’s h-transform established in Theorem 1 is only guaranteed if an exact solution to the optimization problem (9) or (11) is found. However, by introducing the Gaussian parametrization of $q_{t|0,T}$, this guarantee is lost. There is a lack of discussion regarding the implications and potential effects of this parametrization on the overall accuracy and validity of the method.
- Moreover, the formulation of the task of solving equation (12) given $q_{t|0,T}$ is misleading. In reality, you are solving an inverse problem here. Given the solution of the PDE in (12), your goal is to reconstruct $u_{t|0,T}$, which is generally an ill-posed problem. For instance, introducing uncertainties into the model description via the drift $b_t$ or diffusion $\Xi_t$ can lead to significant challenges. This scenario is likely when applying the approach in practice under model misspecification. A critical question to address is the robustness of the proposed framework against model misspecification, such as small perturbations in $b_t$.
- The claim that the approach is simulation-free is not entirely clear. In line 173, it is stated that the Gaussian parametrization allows for the generation of arbitrary samples. However, this assumes that the states $x_{t|0,T}$ are independent for all $t$, which is not true when $x_{t|0,T}$ is a solution of the SDE in (10). The dependence between the states needs to be accounted for in the sampling process, and this aspect seems to be overlooked. 
- I am wondering what the measured quantities in the numerical experiments actually reveal about the correctness of the proposed approach. Evaluating the maximum energy and the likelihood might not provide sufficient information about the accuracy of the estimated distribution. If the sampled paths only follow high likelihood regions, this does not necessarily indicate that the correct distribution has been captured. The objective should be to estimate quantities of interest related to the conditioned SDE in (6). However, it is not clear whether the simulated paths correspond to accurate simulations of (6). For instance, how does the method perform when estimating rare event probabilities or other Monte Carlo estimators with respect to (6)? This aspect needs to be thoroughly addressed to validate the effectiveness of the proposed approach.
questions: questions: - To enhance the practical applicability of the proposed framework, it is crucial to study the impact of inexactly solving the variational problems (9) or (11). Specifically, how do errors propagate to the conditioned dynamical system driven by the SDE (6) when an approximate solution of (9) or (11) is used? Understanding this error propagation is essential for assessing the robustness and reliability of the proposed method in practical scenarios.
- Is the Optimization problem in (9) well-posed? This means, is the optimal solution unique?
- In several instances, there are missing commas in the notation for the inner product (e.g., equations (9b) and (12)).
- Many important mathematical details and assumptions are missing. Most importantly, what are the assumptions on the drift vector field $b_t$ to ensure the well-posedness of the proposed scheme?
- Regarding Section 3.2.2: Instead of introducing $\xi_{\min}$ I suppose that one could also directly work with a pseudoinverse when defining $G_t^{-1}$.
- In Proposition 4, should it be $u_{t|0,T}^{(k)}$ on the right hand side of the equation for $u$?
- When referring to probability density functions, it's important to use proper notation. Instead of $\rho(x_t=x)$, it would be clearer to use $\rho_t(x)$ depending on the context.
- In the abstract, you claim that no "inefficient" importance sampling estimators are required. However, it would be beneficial to compare this approach with such estimators to demonstrate its advantages more clearly.
limitations: limitations: - There is no discussion about error propagation arising due to the Gaussian parametrization.
- The approach is limited to conditioning the sample path on point sets $x_T = B$.
- Limited statistical evidence is presented in the numerical experiments.
paper_summary: paper_summary: The Authors of this paper tackle the problem of sampling conditioned SDEs with a specific interest in "transition path sampling", i.e. sampling a Langevin-type SDE undergoing a transition between an initial state (or set of states) $A$ and a final target set $B$. Sampling transition paths efficiently can provide a tremendous boost to research in catalysis or drug design. In this paper, the Authors model the transition path distribution with either (i) a parametrized Gaussian process or (ii) a parametrized mixture of Gaussians processes. These families of priors are then optimized by leveraging a variational formulation of the problem developed by the Authors starting from Doob's $h$-transform.
strength_weakness: strength_weakness: I find the paper overall clearly written, but I had to go through section 3.2 (Computational Approach) several times to grasp how the method can be deployed in practice. Specifically, I think that it can be improved the explanation of the fact that, because of the Gaussian prior, you only need to model the transition path probability and not the $h$-transform itself, as well as the explanation of the actual optimization step (Reparametrization of Gradients).

The subscript $0,T$ is introduced for the first time in Eq. (6) without explanation and used throughout the manuscript to label variables related to the conditioned process. It might be useful to clarify this notation explicitly. 

The experimental evaluation is somewhat limited, especially concerning the Chignolin protein. Over the years, many transition path sampling strategies have been developed, as well as very much related enhanced sampling techniques. It would be nice to have a comparison also to different baselines, as well as a discussion on the computational complexity and actual running times of the baselines.
questions: questions: Is the $h$-transform related to the committor function in some way? If so, how does your Thm 1 relates to the variational formulation of the committor (see e.g. Eq. 20 of "Transition Path Theory and Path-Finding Algorithms for the Study of Rare Events" by Weinan E and Eric Vanden-Eijnden)
limitations: limitations: Limited experimental evaluation
paper_summary: paper_summary: This paper proposes a variational formulation of Doob's h-transform, which characterizes the distribution over paths with a given endpoint. Instead of relying on potentially wasteful sampling approaches, the authors propose directly optimizing a tractable variational distribution over transition paths which satisfy the initial and terminal conditions by design. This approach reduces the search space over trajectories and avoids trajectory simulation. Experiments on real-world molecular simulation and protein folding tasks demonstrate the applicability of the approach.
strength_weakness: strength_weakness: My main concern comes from trying to interpret the experimental results and in particular judging the performance compared to MCMC. It seems quite clear from Tables 1 and 2 that the variational approach requires fewer calls to the potential energy function than MCMC, however their performance differences are harder to judge in my opinion. 
More specifically, in Table 1 the standard deviations are so large that there is basically no meaningful statistical difference between the shown results especially for the Max Energy, but also the Log-likelihood. Now this might mean that MCMC and the presented method both perform well, but its surprising to me that there is that much inherent variation. Similarly in Table 2, the variance for the Max Energy of MCMC in the first line is huge. 
I'm also a bit confused by the Max Energy increases when using a mixture in Table 2.
To better understand how the performance of the variational approach improves during training, it would be nice to see a plot that shows the Max Energy as a function of the training epochs.

Finally, I would have liked to see a discussion of the limitations of the approach. Section 6 has Limitations in the title, but does not actually discuss them in any way beyond extensions of the proposed method.
questions: questions: - How do you choose the number of mixture components in practice?
- Why does the MaxEnergy increase when using a mixture distribution as a variational approximation in Table 2?
- How do you explain the huge variance of the Max Energy of MCMC (variable length) in the first row of Table 2?
- What are the main limitations of the current approach?
limitations: limitations: Section 6 has "Limitations" in the title, but limitations are not discussed, only future work.
paper_summary: paper_summary: The paper is concerned with sampling trajectories with a terminal condition. For stochastic processes governed by a Brownian Motion, Doob's h-transform gives an posterior SDE that leads to samples with the final condition. However, estimating the h-function that is needed for the posterior SDE usually involves simulating trajectories, which is inefficient if the terminal condition is rarely reached.
The authors propose a simulation-free variational optimization method to estimate the h-function based on a least action principle and Gaussian approximations to the marginal densities.
strength_weakness: strength_weakness: 1. I found the path histograms in Figure 2 to be too cluttered. I would propose adding less samples.

2. The paper misses a learning curve. It would in general be interesting to have more training details.

Minor Weaknesses:

In Chapter 3 and in the appendix, the authors change from trajectory length T to the unit interval. There should be a sentence that explains this change.
questions: questions: I could not totally follow the derivation in the appendix. Can you explain how to get to Equation 25 from the previous equation in line 548?
I do not see why second term in equation 25 gets subtracted, while in line 548 only the last term has a minus.
limitations: limitations: The authors adequately adressed the limitations of their work.

Comments: title: Raising the score
Clarifying Misunderstanding of Test vs. Train Sampling
Thanks for the response
References
comment: Thank you for your response. I will increase my original score by one.
Thank you for your prompt response and thorough clarification of my concerns. I am happy to raise my score to 6. I did not give a higher score because I believe further effort is needed to enhance the presentation of the proposed framework by incorporating the points discussed in the rebuttal.
**Testing vs. Training**
> I might be misunderstanding the entire simulation-free approach, but don’t you sample independently at each time 𝑡 from the marginal distributions to generate paths of the conditioned SDEs?  

Simulation-free refers to our *training* method, where our objective in Thm 1 only requires samples from the time-marginals $q_{t|0,1}$ of the conditioned SDE. This justifies neglecting full trajectory information and sampling directly from $q_{t|0,1}$ (without SDE simulation) during training.

See Algorithm 1 of the general-response PDF or cell 20 of the [anonymized code](https://anonymous.4open.science/r/TPS-Doob-843E/notebooks/tps_gaussian.ipynb).

> … don’t you sample independently at each time $𝑡$ from the marginal distributions to generate paths of the conditioned SDEs?  Generally, how do you ensure the continuity of the sampled paths?

We **do not** sample each marginal independently to generate paths of the conditioned SDE.   Instead, our Algorithm 2 in the general-response PDF outlines our sampling approach.   

For a given set of time marginals $q_{t|0,1}$, we can calculate the appropriate drift $u_{t|0,1}$ for an SDE with given diffusion coefficients.   **At generation (test) time, we simulate this SDE**: $dx_{t|0,1} = u_{t|0,1}(x_{t|0,1}) dt + \Xi_t dW_t$ starting from $x_0 \sim \mathcal{N}(A, \sigma_{\text{min}} \mathbb{I})$.   Thus, our sampled paths are continuous (up to the SDE solver error).

Note that $u_{t|0,1}$ implies a control drift term $v_{t|0,1}$ (in Eq. 10) as in Eq. 14, but we do not need to evaluate the costly drift term $b_t$ at generation time.

See cells 24-25 of [anonymized code](https://anonymous.4open.science/r/TPS-Doob-843E/notebooks/tps_gaussian.ipynb).

> (Please correct me when I am misunderstanding something) 

Thank you for pinpointing this misunderstanding!   We hope that the above points help to address these concerns.

Based on your feedback, we agree that the following changes will greatly improve the final manuscript:
- Discussion of the sampling procedure (and Algorithm box) will be included in the main text
- Distinctions between the time marginals $q_{t|0,1}$ of the conditioned process and the full path measure or SDE will be carefully emphasized throughout.
- We will visualize individual trajectories in Figure 2 to show that the paths are continuous and similar to MCMC trajectories.
- We will modify line 47 to “our training method is simulation-free”.   It appears that all other usage of “simulation-free” explicitly refers to training or optimization, rather than generation.

**Evaluating Sampling Trajectories**
> I would be less critical of my concern if there were clear empirical experiments to support the correct sampling distribution.

Thank you for this suggestion. For the Müller-Brown experiment in Figure 2, we report Wasserstein-1 distances between (i) samples from expensive MCMC simulation (which we treat as ground-truth) and (ii) sample trajectories generated using Algorithm 2 with our learned model.    We report the mean and std of the W1 distance across discrete $t$:

| Wasserstein W1| Value |
|-----------------------|-------|
| Mean                  | 0.1251 |
| Std                   | 0.0392 |
| Median                | 0.1130 |
| Min                   | 0.0393  |
| Max                   | 0.2115 |

We will report further numbers and plots in the revised version of the manuscript. 

 

**Conditioning on Dirac Deltas**

> …challenges associated with conditioning on point sets. Doesn’t it impose technical challenges when applying Doob‘s h-transform for conditioning on null sets?

In the case of the conditioning on a point mass $\delta(x-B)$, the h-function becomes a density, i.e.
$h(y,t) := \rho(x_T = B | x_t = y)$ is the density of the transition probability $\mathbb{P}(x_T \in dx| x_t = y)$.
Conditioning on point sets is commonly used [1 Thm 7.11,2,3], and indeed, all the derivations in Appendix A hold if we take $h(y,t):=\rho(x_T = B | x_t = y)$.  We apologize for the confusion, we will clarify it in the next version.

[1] Särkkä, Simo, and Arno Solin. Applied stochastic differential equations. Vol. 10. Cambridge University Press, 2019.

[2] Heng, J., De Bortoli, V., Doucet, A. and Thornton, J., 2021. Simulating diffusion bridges with score matching. arXiv preprint arXiv:2111.07243.

[3] Liu, Xingchao, Lemeng Wu, Mao Ye, and Qiang Liu. "Let us build bridges: Understanding and extending diffusion generative models." arXiv preprint arXiv:2208.14699 (2022).
Thank you very much for the detailed response. 

I might be misunderstanding the entire simulation-free approach, but don’t you sample independently at each time $t$ from the marginal distributions to generate paths of the conditioned SDEs? (Please correct me when I am misunderstanding something) While this sampling scheme might be consistent under certain assumptions on the drift and diffusion coefficients, it certainly isn’t universally applicable. Generally, how do you ensure the continuity of the sampled paths? This concern remains relevant regardless of whether you employ Gaussian or Gaussian mixture parameterizations.

I would be less critical of my concern if there were clear empirical experiments to support the correct sampling distribution. In particular, the synthetic data example presents an opportunity to conduct a statistically robust case study.

Additionally, I remain concerned about the lack of explicit mathematical assumptions and the challenges associated with conditioning on point sets. Doesn’t it impose technical challenges when applying Doob‘s h-transform for conditioning on null sets? For example, the application of Jamison (1975) requires a positive function h (as also noted in line 514).

While I find the proposed framework very promising, I am maintaining my score due to the concerns I've outlined.
I thank the Authors for their thorough rebuttal and clarifications. I will keep my evaluation to "Accept".
Thank you for your response. I feel my original score is still appropriate.
[1] Blei, Kucukelbir, McAuliffe.  “Variational Inference:  A Review for Statisticians”, 2016. 

[2] Jónsson, H., Mills, G. and Jacobsen, K.W., 1998. Nudged elastic band method for finding minimum energy paths of transitions. In Classical and quantum dynamics in condensed phase simulations (pp. 385-404).

[3] Weinan, E., Ren, W. and Vanden-Eijnden, E., 2004. Minimum action method for the study of rare events. Communications on pure and applied mathematics, 57(5), pp.637-656. 

[4] C. Dellago, P. G. Bolhuis, and P. L. Geissler, 2006, “Transition Path Sampling Methods,” Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology Volume 1. Springer Berlin Heidelberg, pp. 349–391. 

[5] Brekelmans, Rob, and Kirill Neklyudov, 2023, "On Schrödinger Bridge Matching and Expectation Maximization." In NeurIPS 2023 Workshop Optimal Transport and Machine Learning.

[6] B. Jamison, 1975, „The Markov processes of Schrödinger“, Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, Bd. 32, Nr. 4. Springer Science and Business Media LLC, S. 323–331. 

[7] Holdijk, L., Du, Y., Hooft, F., Jaini, P., Ensing, B. and Welling, M., 2024. Stochastic optimal control for collective variable free sampling of molecular transition paths. Advances in Neural Information Processing Systems, 36.