{"id":"11ef8f23-4e4f-4331-99f7-ac2d07766a1f","arxiv_id":"2411.08518","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":4.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"Known Girsanov and Bismut-Elworthy-Li identities are recast as explicit grid-free Monte Carlo algorithms for coupled Fokker-Planck and Hamilton-Jacobi-Bellman systems, with a neural-network demonstration on a Schrödinger bridge problem.","lead":"This note turns two classical stochastic calculus tools, the Girsanov theorem and the Bismut-Elworthy-Li formula, into Monte Carlo algorithms for the Fokker-Planck and Hamilton-Jacobi-Bellman equations of stochastic optimal control. The methods are tested on Schrödinger bridge problems and on a machine-learning prototype, all without spatial discretization.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Underdamped Algorithm 2 samples the non-reversible bridge of (27) with a naive backward Euler, missing the score correction; the conditional expectation (32) is biased and the underdamped FP integration is unsupported.","rationale":"The reader's weakest assumption correctly flags the backward-time sampler, but the precise failure is confined to the underdamped case: the overdamped Algorithm 1 simulates a driftless Brownian auxiliary process, whose backward Euler is exact at grid points, so no score correction is needed there. In contrast, Algorithm 2 simulates the underdamped process (27) backward with the same drift, which samples neither the backward SDE (25) nor the true bridge of the forward process I. The divergence-free property of the drift preserves phase-space volume for the deterministic part but does not make the stochastic process reversible; the time-reversed bridge requires the score term. This is a load-bearing gap because the underdamped Fokker-Planck integration is one of the paper's two central advertised methods, and Figure 3, its only underdamped FP validation, relies on the biased sampler and on the same group's perturbative reference [55]. The machine-learning example and the overdamped BEL results may be unaffected, but the central claim of a working underdamped Girsanov algorithm is not established. A correct version would need to simulate the bridge with the density-gradient drift or use a different reference measure, and the numerics would need to be redone. Hence the verdict should move from conditional to reject for the current form.","tokens_in":27866,"tokens_out":31352,"duration_ms":298547,"concrete_test":"Compute the exact solution of (23) for (22) with ∂U=0, m=τ=β=1, and Gaussian initial density with unit variance; the solution is Gaussian with known covariance (Ornstein-Uhlenbeck momentum plus integrated position). Run Algorithm 2 from a fixed terminal x at t=1 with h=0.05, 0.01, and 0.002, using at least 10^6 trajectories, and compare the estimate to the exact density. If the error does not vanish as h→0, the backward Euler is not an unbiased surrogate for the bridge in (32), and the underdamped method fails.","verdict_should_be":"REJECT","load_bearing_attack":"Section 3.1 claims the path measure I generated by (27) is time-reversal invariant because the drift is divergence-free. This conflates Liouville volume preservation with stochastic reversibility. For a diffusion with constant diffusion coefficient, the reverse-time bridge drift is -b + σσᵀ∇log p_{T-t}; even for the free case ∂U=0, the backward update in Algorithm 2, p_{t-h} = p_t - sqrt(2m/(τβ) h) ϵ, is not the Brownian-bridge conditional law of p_{t-h} given p_t. The true conditional mean is not p_t and the conditional variance is reduced; the omitted contribution is exactly the density-gradient correction. Consequently, Algorithm 2 does not sample the I-bridge in (32), so the underdamped estimates in Fig. 3 are biased at any fixed h and do not converge to the claimed pt. The τ-scaling in (32), which puts τ in the numerator so the weight diverges as τ→∞, is a related but secondary defect; the missing score correction is the deeper problem.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes two grid-free Monte Carlo methods for PDEs arising in stochastic optimal control: (i) a Girsanov/Feynman-Kac representation of the Fokker-Planck equation for a diffusion driven by a mechanical potential, with separate treatments for overdamped (Section 2) and underdamped (Section 3) dynamics and corresponding Algorithms 1 and 2; and (ii) Bismut-Elworthy-Li (BEL) representations for the gradient of the solution of a Hamilton-Jacobi-Bellman equation, in both non-degenerate (Section 4) and degenerate (Section 5) settings, with Algorithms 3 and 4. The methods are demonstrated on Schrödinger bridge problems and on a neural-network-trained optimal control protocol (Algorithm 5, Section 6). The overdamped results are checked against analytic examples and against the independent Caluya-Halder iterative method; the underdamped FP integration is checked against a perturbative reference from the authors' own prior work.","tokens_in":28115,"tokens_out":33157,"duration_ms":316882,"significance":"If correct, the proposed representations would provide attractive, dimension-robust Monte Carlo alternatives to spatial discretization for the coupled FP/HJB systems that arise in stochastic thermodynamics and control. The overdamped Girsanov representation, the BEL formulas, and the analytic and numerical verifications in Sections 2.2, 4.3, 4.4, and 5.2 are useful and appear sound. The code availability and the external comparison with Caluya-Halder in the overdamped cases are strengths. However, the underdamped FP algorithm in Section 3 is not correctly sampling the bridge measure it claims to use, and the displayed formula (32) contains a τ-scaling error. As a result, the numerical underdamped demonstration in Fig. 3 does not validate the method as presented.","major_comments":[{"comment":"Algorithm 2 samples a backward Euler discretization of the forward process (27) starting from the terminal point x. This is not the conditional law of (27) given X_t = x. For a diffusion with generator b·∇ + (1/2)σσ^T:∇∇, the reverse-time bridge drift is -b + σσ^T ∇ log p_{T-t}; the divergence-free condition on the deterministic part only gives Liouville volume preservation, not time-reversal invariance of the path measure. Concretely, in the free case ∂U = 0, the update for p in Algorithm 2 has conditional mean p_t, whereas the true Brownian-bridge conditional mean at the earlier time is not p_t. The estimator therefore has a bias of order one at every fixed step size and does not converge to the conditional expectation in (32) as h → 0. This undermines the underdamped FP integration claim and the evidence in Fig. 3.","section":"Section 3.1 and Algorithm 2"},{"comment":"The Girsanov weight in (32) has an incorrect scaling in the parameter τ. The derivation immediately preceding (32) gives, for m = 1, an exponent with β/(4τ) ||p_s||^2 and √β/(2√τ) ⟨dw_s, p_s⟩, whereas (32) states τβ/4 ||p_s||^2 and √(τβ)/2 ⟨dw_s, p_s⟩. These two expressions agree only when τ = 1. Since the numerical example in Section 3.2 uses τ = 1, the error is not visible in Fig. 3, but the formula as written is false for general τ and needs correction.","section":"Equation (32)"}],"minor_comments":[{"comment":"The underdamped FP reference solution is taken from the authors' own perturbative work [55], so the agreement shown in Fig. 3 is not an independent validation of the Monte Carlo scheme. This is a limitation, particularly because the algorithm itself is in question.","section":"Section 3.2, Fig. 3 caption"},{"comment":"The phrase \"divergence-less\" should be \"divergence-free\"; the wording conflates a property of the deterministic flow (Liouville volume preservation) with a property of the path measure, which is precisely the source of the sampling error in Algorithm 2.","section":"Section 3.1"},{"comment":"The index notation \"∂U_{t_{n-1-i}}\" in the pseudocode is confusing; clarity would be improved by writing the potential derivative as a function of q_{t_{n-1-i}}.","section":"Section 2.2.3, Algorithm 1"},{"comment":"There are several typographical errors, including \"auxilliary\" (Section 2.2.3), \"Kullback-Liebler\" (Section 2.2.4), and \"divergence-less\" (Section 3.1).","section":"General"}],"recommendation":"major_revision","confidential_remarks":"The paper has two independent thrusts. The overdamped Girsanov and BEL parts are sound and well documented, with external benchmarks for the overdamped case. The underdamped FP section, however, has a load-bearing algorithmic flaw and a formula-level τ error. The authors could fix the underdamped section by either (i) using the backward process (25) directly, for which the Euler-Maruyama discretization is consistent with the P^♭ measure, or (ii) implementing a correct diffusion-bridge sampler for the I-measure with the required score correction. Without such a change, the underdamped FP claim should be withdrawn or substantially weakened. The self-referential reference solution in the underdamped example also weakens the paper's demonstration."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Quick take: the overdamped half of this note is a genuinely useful recipe, with honest external validation; the underdamped half has a real sampling error that invalidates Fig. 3 as it stands.\n\nWhat's actually new: They take Girsanov and BEL, both classical, and turn them into concrete grid-free Monte Carlo algorithms for mechanical-force Fokker-Planck equations and degenerate BEL gradients, then demonstrate on Schrödinger bridge problems and a small ML training loop. The code is on GitHub, the analytic examples check out, and the overdamped FP and BEL comparisons against Caluya-Halder are convincing. That is legitimate, useful within-subfield work.\n\nThe soft spot is Section 3. The claim in Section 3.1 that the path measure I of (27) is time-reversal invariant because the drift is divergence-free is wrong: divergence-free gives Liouville volume preservation, not stochastic reversibility. The correct reverse-time bridge drift needs the score correction. Algorithm 2 runs (27) backward with the same drift and noise, so it samples the unconditioned backward process, not the bridge conditional on x_t=x that (32) requires. In the free case the backward momentum step misses the Brownian-bridge mean reversion. So the underdamped estimates in Fig. 3 are biased at any fixed h, and the claim that the method integrates the underdamped FP is unsupported. There is also a τ-dependence problem in (32): the weight has τ in the numerator where the Girsanov derivation gives 1/τ, so it diverges as τ→∞; whoever reviews this should check that algebra carefully. The self-citation to [55] for the underdamped benchmark is fine as a reference solution, but it is not an independent check, so the numerics carry less weight.\n\nThe reader's weakest assumption — that reverse-time Euler-Maruyama needs a density-gradient correction — is actually wrong for the overdamped Algorithm 1, because there they sample the drift-free measure Q from the terminal point, which is a legitimate backward Markov transition. The bias accusation lands only on the underdamped algorithm.\n\nBEL part is fine; the degenerate construction with ℓ_u is plausible and the analytic example verifies. The ML prototype is a prototype, honestly labeled, no quantitative validation, so I don't hold that against it.\n\nBottom line: this deserves a serious referee. The overdamped half is sound and usable, and the underdamped error is fixable (sample P^♭ of (25) instead of pretending I is reversible, and fix the τ). I'd send it out, with a request to correct Section 3 before acceptance.","headline":"The overdamped Girsanov integration and BEL gradients are solid and useful, but the underdamped Fokker-Planck algorithm samples the wrong bridge measure and the τ-scaling in (32) is off, so that part needs rework.","tokens_in":28637,"tokens_out":15652,"would_cite":true,"duration_ms":154234,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["35Q84","49L20","60H30","65C05","65C30"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper derives exact path-average representations for Fokker-Planck densities and value-function gradients, turning stochastic optimal control into grid-free Monte Carlo.","keywords":["Fokker-Planck equation","Girsanov theorem","Bismut-Elworthy-Li formula","Hamilton-Jacobi-Bellman equation","Schrödinger bridge","Monte Carlo methods","underdamped Langevin dynamics","stochastic optimal control"],"falsifier":"For a harmonic potential, where the transition density is Gaussian and exactly known, simulate the backward scheme (15), reconstruct the density at $t_\\iota$, and compare with the analytic forward kernel; any systematic bias as the time step shrinks would indicate that the reverse-time sampler is not an unbiased surrogate for the true backward transition kernel. Alternatively, in the quartic example of Section 2.2.3, record the empirical variance of the Girsanov weight $e^{-g}$ across samples and over longer horizons; exponential growth of the variance would show that the Monte Carlo estimate is dominated by rare paths.","tokens_in":27609,"feed_emoji":"🎲","tokens_out":15667,"duration_ms":140929,"temperature":0.7,"pith_summary":"The paper's aim is to turn two PDEs that arise in stochastic optimal control into quantities that can be estimated by Monte Carlo sampling, with no spatial grid. For Fokker-Planck equations whose drift is the gradient of a mechanical potential, it derives an exact representation of the probability density as an expectation over paths of an auxiliary diffusion, weighted by a Girsanov factor, for both overdamped and underdamped (inertia-plus-friction) dynamics. For the Hamilton-Jacobi-Bellman equation, it uses the Bismut-Elworthy-Li formula to express the gradient of the value function as an expectation of the terminal payoff times a stochastic integral, which is the quantity the stationarity condition needs to fix the optimal protocol. Both representations are evaluated numerically by Euler-Maruyama sampling of the underlying diffusion, so the cost scales with the number of sampled trajectories rather than with grid points per dimension. The paper demonstrates the combined scheme on Schrödinger bridge problems (steering a distribution between assigned endpoints with minimal divergence from a free diffusion) and on a neural-network-trained optimal control protocol, checking against iterative and perturbative reference solutions.","feed_headline":"Monte Carlo replaces grids for Fokker-Planck and control gradients","feed_subtitle":"Fokker-Planck densities and value gradients become samples of auxiliary diffusions, with Schrödinger-bridge and neural-network demos.","key_machinery":"The machinery has two parts. First, the factorization $p_t(q) = e^{-\\beta U_t(q)} f_t(q)$ recasts the forward Fokker-Planck equation as a backward Kolmogorov equation for $f_t$ along a time-reversed diffusion; Girsanov's change of measure then replaces that backward path measure by the measure of the free Wiener process, so the solution becomes an average of the initial density times the exponential Radon-Nikodym weight, accumulated as a running cost on backward Euler-Maruyama trajectories. Second, the Bismut-Elworthy-Li formula, derived through Malliavin's integration-by-parts identity, writes a directional derivative of the value function as an expectation of the terminal payoff multiplied by the stochastic integral $\\int_t^{t_f}\\langle dw_u, A_u^{-1}(x_u) x^*_{u,t} e_i\\rangle$, where $x^*_{u,t}$ is the derivative of the stochastic flow with respect to its initial point; the algorithms accumulate this weight along forward trajectories of the controlled diffusion. The same weight appears in the degenerate-diffusion version, with an auxiliary variation field chosen to enforce the required boundary conditions in phase space.","core_discovery":"The central discovery is that a Fokker-Planck equation whose drift is the gradient of a mechanical potential can be solved pointwise by a conditional expectation over paths of an auxiliary diffusion. For the overdamped equation (3), the solution is\n\n$$p_t(q) = $E^{{\\mathbb{Q}}$}\\Bigl[ p_{t_\\iota}(q_{t_\\iota}) \\, $e^{{-\\frac{\\beta}}${2}\\int_{t_\\iota}^{t}\\bigl(\\langle dq_s,(\\partial U_s)(q_s)\\rangle + \\frac{\\mu}{2}\\|(\\partial U_s)(q_s)\\|^2 ds\\bigr)} \\Bigm| q_t = q \\Bigr],$$\n\nwith the expectation over paths of the backward free diffusion; equation (32) is the underdamped analogue, in which the Girsanov weight absorbs the force and friction terms. For the Hamilton-Jacobi-Bellman equation (34), the Bismut-Elworthy-Li formula expresses each component of the gradient of the value function as an expectation of the terminal data and the running cost multiplied by stochastic-integral weights, equation (46), and Section 5 extends this to degenerate diffusions such as inertia-plus-friction dynamics. These identities are exact; Algorithms 1-5 evaluate them by Euler-Maruyama sampling of the relevant diffusion, which is what makes the integration grid-free.","pith_inferences":["Beyond the paper: the reciprocal-process identity (49) suggests a way to estimate the score function itself by conditioning the Bismut-Elworthy-Li weight on both endpoints, which would give a particle scheme for probability-flow equations without training a score network.","Beyond the paper: the Girsanov weight is an exponential of a path action, so its variance can be expected to grow with the horizon and with dimension; importance sampling or weight clipping may be needed before the method beats spectral or tensor-grid solvers in high-dimensional problems.","Beyond the paper: a direct numerical stress test would apply Algorithm 1 to a harmonic potential with an analytically known Gaussian transition density and monitor the empirical variance of the Girsanov weight as $\\beta$, $\\mu$, or the horizon is increased.","Beyond the paper: the degenerate BEL construction of Section 5 could serve as the training signal for an underdamped machine-learning bridge, extending Algorithm 5 to inertial dynamics; the paper suggests this as a possible generalization but does not demonstrate it."],"forward_implications":["The pointwise Fokker-Planck estimate costs one reverse-time trajectory per sample, so the method has no spatial mesh and can in principle be used in high-dimensional state spaces where grid-based solvers become unaffordable.","The same Girsanov construction covers underdamped (inertia-plus-friction) Langevin dynamics, not just overdamped diffusion, giving a Monte Carlo route to joint position-momentum densities in stochastic thermodynamics.","The Bismut-Elworthy-Li gradient estimator makes the stationarity condition (19) and its underdamped analogue (31) directly evaluable from samples, so optimal control updates do not require solving a nonlinear PDE on a grid.","Combining the two estimators yields an iterative Schrödinger bridge solver that the paper validates against a proximal iteration in the overdamped case and against perturbative predictions in the underdamped case.","The neural-network protocol of Algorithm 5 is a direct corollary: at each gradient-descent step the BEL estimator supplies the gradient of the value function, and the Girsanov estimator supplies the final density used to update the Lagrange multiplier."],"supporting_citations":[{"why":"supplies the martingale, backward stochastic calculus, Girsanov change of measure, and Dynkin formula facts used throughout the derivations.","marker":"[33]"},{"why":"defines backward stochastic differential equations and their martingale property, used in Propositions 2.1 and 3.1.","marker":"[36]"},{"why":"justifies the post-point stochastic-integral convention used to put the Girsanov weight into pre-point form in equation (11).","marker":"[42]"},{"why":"originates the Bismut and Malliavin-calculus gradient formula that the paper adapts to value functions.","marker":"[5]"},{"why":"states the heat-semigroup derivative formulas behind the non-degenerate Bismut-Elworthy-Li representation (46).","marker":"[20]"},{"why":"supplies the Malliavin integration-by-parts identity used to derive the gradient formula in Section 4.1.","marker":"[48]"},{"why":"provides the Bismut formulas for degenerate stochastic Hamiltonian systems that Section 5 extends.","marker":"[66]"},{"why":"is the iterative half-bridge reference solution used to validate the overdamped Fokker-Planck and value-function integrations.","marker":"[9]"},{"why":"is the proximal-gradient solver used as the comparison baseline for the time-independent-drift Fokker-Planck example.","marker":"[8]"},{"why":"supplies the perturbative optimal protocol and Gaussian-boundary reference data used in the underdamped examples.","marker":"[55]"}],"fun_headline_variants":["Sampling replaces meshes for Fokker-Planck and control gradients","Grid-free integration via Girsanov and Bismut-Elworthy-Li formulas","Path sampling solves optimal control PDEs without discretization","Monte Carlo PDE solver taps diffusion paths for gradients","Exact probabilistic formulas yield grid-free optimal control"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the exponential path weight has finite variance and that the reverse-time Euler-Maruyama sampler is an unbiased representation of the backward transition kernel; if either fails, the Monte Carlo averages are inaccurate even though the analytic representations are exact.","fun_headline_variants_meta":{"raw":{"variants":["Sampling replaces meshes for Fokker-Planck and control gradients","Grid-free integration via Girsanov and Bismut-Elworthy-Li formulas","Path sampling solves optimal control PDEs without discretization","Monte Carlo PDE solver taps diffusion paths for gradients","Exact probabilistic formulas yield grid-free optimal control"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00024,"raw_usage":{"total_tokens":1536,"prompt_tokens":984,"completion_tokens":552,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":600,"completion_tokens_details":{"reasoning_tokens":466}},"tokens_in":600,"tokens_out":552,"duration_ms":5403,"temperature":1.0,"reasoning_tokens":466,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T21:33:32.916776+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"For a harmonic potential, where the transition density is Gaussian and exactly known, simulate the backward scheme (15), reconstruct the density at $t_\\iota$, and compare with the analytic forward kernel; any systematic bias as the time step shrinks would indicate that the reverse-time sampler is not an unbiased surrogate for the true backward transition kernel. Alternatively, in the quartic example of Section 2.2.3, record the empirical variance of the Girsanov weight $e^{-g}$ across samples and over longer horizons; exponential growth of the variance would show that the Monte Carlo estimate is dominated by rare paths.","supporting_citations":[{"cited_title":"Klebaner","cited_arxiv_id":null,"evidence_quote":"supplies the martingale, backward stochastic calculus, Girsanov change of measure, and Dynkin formula facts used throughout the derivations."},{"cited_title":"On backward stochastic differential equations","cited_arxiv_id":null,"evidence_quote":"defines backward stochastic differential equations and their martingale property, used in Propositions 2.1 and 3.1."},{"cited_title":"Géométriedifférentiellestochastique,II","cited_arxiv_id":null,"evidence_quote":"justifies the post-point stochastic-integral convention used to put the Girsanov weight into pre-point form in equation (11)."},{"cited_title":"Large Deviations and the Malliavin Calculus","cited_arxiv_id":null,"evidence_quote":"originates the Bismut and Malliavin-calculus gradient formula that the paper adapts to value functions."},{"cited_title":"The Malliavin calculus and related topics","cited_arxiv_id":null,"evidence_quote":"supplies the Malliavin integration-by-parts identity used to derive the gradient formula in Section 4.1."},{"cited_title":"DifferentialEquations.jl–a performant and feature-rich ecosystem for solving differential equations in Julia","cited_arxiv_id":null,"evidence_quote":"supplies the perturbative optimal protocol and Gaussian-boundary reference data used in the underdamped examples."}],"review_version":1}