{"id":"02f9dc8e-a083-4c27-adeb-8e0c062a454a","arxiv_id":"1908.00533","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"A Sinkhorn-type proximal recursion, proved contractive in the Thompson metric, computes transient Fokker-Planck and McKean-Vlasov densities on weighted point clouds without spatial discretization.","lead":"This paper introduces a numerical method for propagating joint probability densities through continuous-time stochastic systems. It replaces spatial grids with a moving cloud of weighted points, updated by an optimal transport recursion.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The central claim lacks a consistency proof: the Euler-Maruyama/Sinkhorn split is not shown to converge to the JKO scheme (or FPK PDE) as h→0, and ε is held fixed.","rationale":"The strongest claim of the paper is that the discretized-in-time but not-in-space proximal recursions are theoretically equivalent to the transport PDEs in the small-time-step limit. The proof infrastructure in Section III-C establishes only that the block-coordinate iteration in Algorithm 1 converges to a fixed point of the finite-dimensional stationarity conditions (33) for a given cost matrix and a fixed entropic regularization. This is an inner-loop result; it says nothing about the outer time stepping. The numerical scheme actually implemented is a splitting: positions are advanced by explicit Euler-Maruyama (39), and then weights are updated by an entropic Sinkhorn step. For the central claim to hold, this split must be consistent with the JKO variational scheme (22) as h→0. The paper contains no such consistency analysis. In addition, the entropic parameter ε is held fixed in all experiments (ε=5×10^{-2} while h=10^{-3}); because the entropic term in (27) is εH(M) rather than hεH(M), the balance between transport regularization and free energy changes as h→0, so it is not even clear that the fixed-ε scheme has the correct h→0 limit. The numerical examples show visual agreement with analytical solutions at a single small h, but report no error metrics or h-refinement; they therefore do not supply the missing consistency evidence. I agree with the reader's weakest assumption. (The reader's additional point about eq. (33b) appears to be a typographical issue in the manuscript text, since the intended fixed point z⊙Γ^T y = ξ⊙z^{-βϵ/h} is consistent with Algorithm 1.) I would keep the conditional verdict: the paper is promising, but the central claim is not yet backed by a proof or an adequate numerical convergence study.","tokens_in":22079,"tokens_out":11716,"duration_ms":113532,"concrete_test":"Run the Ornstein-Uhlenbeck benchmark of Section IV-A1 at a fixed final time T=1 (a=1, β=1, ρ0=N(5,0.04)), using Algorithm 1 with ε=5×10^{-2}, N=400, tolerance δ=10^{-3}, and h=10^{-1}, 10^{-2}, 10^{-3}, 10^{-4}, 10^{-5}. Reconstruct a smooth density from the weighted point cloud (e.g., kernel density estimate with a bandwidth rule held fixed across h) and compute the L1 error against the analytical Gaussian N(µ(T), σ²(T)) given by the formulas in Section IV-A1. If the L1 error does not decrease monotonically as h decreases (ideally O(h)), then the split scheme is not consistent with the FPK PDE in the small-h limit, and the paper's central claim fails.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The paper's central claim is that Algorithm 1, together with Euler-Maruyama position updates, computes transient joint PDFs of FPK/McKean-Vlasov flows, with the small-time-step limit being the underlying transport PDEs. What is proven (Theorem 3 and Corollary 4) is only that the inner Sinkhorn weight recursion converges to the fixed point of the discrete optimality conditions (33) for fixed h, ε and a fixed cost matrix C_k; the proof itself also asserts Perron-Frobenius contractivity of the positive linear map Γ_k without a demonstration, although this can be repaired since the composite map is strict due to θ1. The load-bearing gap is the outer recursion: the split scheme replaces the JKO variational step (22) by (i) an explicit Euler-Maruyama move of the support points (39) and (ii) a weight update that solves the entropic-regularized OT problem (27) with fixed ε. No theorem or error estimate connects this two-step update to the JKO scheme or to the FPK PDE as h→0. Moreover, ε is not sent to zero; since the entropic term εH(M) in (27) is not scaled by h, as h→0 with ε fixed the relative weighting between the transport regularization and hF(ϱ) changes, so the limit is not obviously the gradient flow of F. The numerical sections report only visual agreement at h=10^{-3}, ε=5×10^{-2}, with no error-versus-h study, so they do not establish consistency. Without this, the central claim that the algorithm approximates the FPK PDE rests on unverified heuristic.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops a computational framework for propagating joint probability density functions (PDFs) of stochastic systems described by Fokker-Planck-Kolmogorov (FPK) and McKean-Vlasov equations. The approach is based on the Jordan-Kinderlehrer-Otto (JKO) gradient-flow formulation, discretizing time while keeping the state space continuous. The proposed algorithm represents the density as a weighted point cloud, advances the support points via Euler-Maruyama, and updates the weights by solving an entropic-regularized optimal transport proximal recursion. The inner fixed-point iteration (Algorithm 1) is proved contractive in the Thompson metric (Theorem 3), and numerical examples for linear Gaussian, nonlinear non-Gaussian, McKean-Vlasov, CIR, and a satellite dynamics problem are compared against analytical solutions or stationary densities.","tokens_in":22439,"tokens_out":7282,"duration_ms":75087,"significance":"If the full algorithmic pipeline were proved convergent to the underlying PDEs as the time step tends to zero, the paper would offer a novel, non-parametric, grid-free method for uncertainty propagation that avoids spatial discretization and its curse of dimensionality. The paper contains explicit algorithmic details, benchmarks against known analytical solutions, and timing measurements. However, the theoretical support provided only covers the inner Sinkhorn fixed-point recursion; the consistency of the outer splitting scheme with the JKO flow is not established. This missing analysis currently limits the significance of the central claims, although the numerical evidence suggests the approach may have merit.","major_comments":[{"comment":"The paper's abstract and introduction claim that the proposed scheme is 'theoretically equivalent to solving the underlying transport PDEs' in the small time-step limit, but the only convergence result proved (Theorem 3 and Corollary 4) concerns the inner fixed-point iteration for fixed h, epsilon, and C_k. No theorem or error estimate connects the two-step update (Euler-Maruyama for positions, then Algorithm 1 for weights) to the JKO scheme or to the FPK/McKean-Vlasov PDE as h tends to zero. The numerical experiments use a single time step h = 10^{-3} and do not report any error-versus-h study, so they do not fill this gap.","section":"III-C and overall"},{"comment":"The entropic regularization parameter epsilon is taken fixed and is not scaled with the time step h. In the objective (27), the term epsilon H(M) is not multiplied by h, whereas hF(rho) is. As h tends to zero with epsilon fixed, the relative weighting changes, and the limit of the proximal step need not be the gradient flow of F. The paper never sends epsilon to zero (e.g., by choosing epsilon = epsilon(h) with a suitable rate) nor proves that the entropically regularized minimizer converges to the unregularized JKO step. Without this, the claimed equivalence of the proximal recursion to the original free-energy gradient flow is not established.","section":"III-A, Eq. (27)"},{"comment":"There is an inconsistency in the meaning of the weights rho^i_k. The text states that rho^i_k 'denotes the value of the joint PDF evaluated at that point,' but Eq. (23) treats rho as a probability vector in the simplex, so the weights are probability masses summing to one. The discrete free energy in (27), equal to <psi + beta^{-1} log rho, rho>, is the negative entropy of the atomic measure, not the continuous free energy integral rho log rho dx. No convergence analysis is provided to show that the discrete weighted point cloud approximates the continuous density as N grows, which is necessary to justify the claim that the algorithm computes the transient joint PDF.","section":"III, paragraph before Eq. (23)"}],"minor_comments":[{"comment":"Equation (33b) is typeset ambiguously; it should read z \\odot (Gamma_k^\\top y) = \\xi_{k-1} \\odot z^{-\\beta\\epsilon/h}. The exponent on z is not clearly displayed.","section":"Theorem 1, Eq. (33b)"},{"comment":"The statement that the positive linear map Gamma_k is contractive in the Thompson metric 'by Perron-Frobenius theorem' is imprecise; the correct result is Birkhoff's contraction theorem for strictly positive matrices. Since the map theta_1 is already strictly contractive, this claim is not needed for the conclusion.","section":"Proof of Theorem 3"},{"comment":"The initialization z0 is chosen randomly. The authors could note that the Banach contraction principle guarantees convergence from any initial point in the positive orthant, so the random choice is immaterial to the convergence.","section":"Algorithm 1"},{"comment":"The computational-time plots (Figs. 9 and 14) do not state the hardware and software environment, which makes the timing results difficult to reproduce or compare.","section":"Section IV and V"}],"recommendation":"major_revision","confidential_remarks":"The central theoretical gap is substantial: the paper proves only the inner fixed-point contraction, not consistency of the full split scheme. The authors should either provide a convergence proof for the splitting as h -> 0 (with a suitable scaling of epsilon) or substantially temper the claims of theoretical equivalence. The numerical evidence is suggestive, but without the analysis, the paper does not currently meet the standard for a math.OC publication. The idea is interesting and may be salvageable through a major revision."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this is a serious attempt to give JKO computational teeth, and the numerics are good enough that I believe the algorithm does something real. But the paper overclaims: the step from the inner Sinkhorn recursion to the PDE is not proved, and the missing piece is not cosmetic.\n\nWhat's new is the specific point-cloud proximal recursion (Euler-Maruyama for support, entropic OT for weights) with a Thompson-metric contraction proof for the inner loop. That combination is not in the JKO/Sinkhorn literature. The extensions to CIR and to a conservative-dissipative satellite problem show real scope.\n\nThe numerics are a genuine plus. The OU, multivariate LTI, stationary bimodal, and CIR examples all match analytic solutions visually, and the 6D satellite marginals are plausible. That is not easy to get from a toy method.\n\nSoft spots, in order of size. First, no consistency theorem for the split scheme: as h goes to zero, the support points move by Euler-Maruyama while the weights solve an entropic OT problem with fixed epsilon. Since epsilon is not sent to zero or scaled by h, the limit is not obviously the JKO gradient flow. The paper proves only that the inner fixed point exists and is unique for each step. The claim that the split approximates the FPK PDE rests on unverified heuristic. A numerical error-versus-h study would partially address this, but it is not in the paper.\n\nSecond, the alleged inconsistency in (33b). I checked the math. If the printed equation is read as z times (Gamma^T y) equals xi times z^{-beta*epsilon/h}, it is consistent with the proof and with Algorithm 1. The typesetting is ambiguous and should be fixed, but it is not a load-bearing flaw. The real flaw is the missing outer convergence.\n\nThird, the Thompson contractivity of Gamma_k is asserted via Perron-Frobenius without proof. This is repairable, since the power map theta_1 is already a strict contraction, so the composite is contractive even if Gamma_k is merely non-expansive.\n\nFourth, the numerical validation would be stronger with quantitative error metrics, a description of how weights relate to the continuous PDF (weights live in the simplex but are compared directly to density values), and wall-clock comparisons against standard PDE solvers. The timing plots are too thin to support the \"remarkably fast\" claim.\n\nWho this is for: researchers in uncertainty propagation, filtering, or computational optimal transport who want a practical alternative to grid-based FPK. The paper deserves a serious referee. I would send it out, with the expectation of major revision: add a consistency argument or a careful numerical convergence study, clarify the point-cloud representation, and clean up the typography. The inner-loop result is worth keeping.","headline":"Promising point-cloud Sinkhorn/JKO density propagation, but the central claim is unproven: no consistency for the outer split, so treat it as a strong heuristic until that is fixed.","tokens_in":22959,"tokens_out":10864,"would_cite":false,"duration_ms":104004,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["49Q22","35Q84","60H35","65C30"],"pacs":[],"model":"deepseek-v4-flash","headline":"A time-only discretization can compute transient joint probability densities of stochastic systems without gridding the state space.","keywords":["Fokker-Planck-Kolmogorov equation","Wasserstein gradient flow","JKO scheme","entropic regularization","Sinkhorn iteration","uncertainty propagation","McKean-Vlasov dynamics","Thompson metric"],"falsifier":"Run the Ornstein-Uhlenbeck benchmark of Section IV-A1 with the same $N=400$ and $\\epsilon=5\\times 10^{-2}$ while shrinking $h$ from $10^{-2}$ to $10^{-6}$, and compare the proximal weighted mean and variance against the exact formulas $\\mu_0 \\exp(-at)$ and $(\\sigma_0^2 - 1/(a\\beta)) \\exp(-2at) + 1/(a\\beta)$; if the error does not systematically decrease as $h \\to 0$, the claimed consistency of the split scheme is false. Separately, compute the Thompson contraction factor of the linear map $z \\mapsto \\Gamma_k z$ for random positive $z$ with a randomly generated cost matrix $C_k$; if any eigenvalue of the associated projective map has modulus at least one, the asserted Perron-Frobenius strict contractivity fails.","tokens_in":21882,"feed_emoji":"🌀","tokens_out":9226,"duration_ms":86552,"temperature":0.7,"pith_summary":"This paper claims that transient joint probability density functions for Itô stochastic differential equations can be computed by discretizing time alone, never the state space. The core idea is to view the Fokker-Planck-Kolmogorov PDE as a Wasserstein gradient flow and to replace it by the Jordan-Kinderlehrer-Otto proximal recursion, which minimizes Wasserstein distance plus a free-energy term at each small time step. The paper shows how to make that infinite-dimensional recursion numerically concrete: an entropic regularization turns each proximal step into a cone-preserving fixed-point iteration whose unique solution gives updated probability weights on a moving point cloud. For systems already in the JKO canonical form, and for several systems that can be transformed into it, the paper argues that this yields fast, nonparametric, grid-free propagation of joint densities. A sympathetic reader would care because this offers a route around the curse of dimensionality that limits conventional PDE solvers.","feed_headline":"Time steps alone can propagate probability densities with no grid","feed_subtitle":"A proximal fixed-point scheme propagates joint PDFs of stochastic systems using weighted particles only.","key_machinery":"The load-bearing object is the JKO proximal recursion, a variational replacement of the Fokker-Planck-Kolmogorov PDE in which each small time step solves $\\arg\\inf_{\\varrho} \\frac{1}{2} W^2(\\varrho_{k-1},\\varrho) + h F(\\varrho)$ with $W$ the 2-Wasserstein metric and $F$ the free energy. Entropic regularization of the inner optimal-transport coupling makes the optimal matrix have the Sinkhorn form $m_{\\mathrm{opt}}(i,j) = \\exp(\\lambda_0(i)h/\\epsilon) \\exp(-C_k(i,j)/(2\\epsilon)) \\exp(\\lambda_1(j)h/\\epsilon)$, algebraicizing the problem. The resulting fixed-point equations are solved by a block-coordinate iteration whose contraction is established in the Thompson metric, $d_T(z,\\tilde{z}) = \\log \\max_i \\max(z_i/\\tilde{z}_i, \\tilde{z}_i/z_i)$, on the positive orthant, using Perron-Frobenius positivity of the kernel $\\Gamma_k = \\exp(-C_k/2\\epsilon)$. The point-cloud locations advance independently by the Euler-Maruyama scheme, so the whole algorithm is a temporal recursion on weighted scattered data rather than a spatial discretization.","core_discovery":"The central claim is that the proximal recursion $\\varrho_k = \\arg\\min_{\\varrho} \\frac{1}{2} W^2(\\varrho_{k-1}, \\varrho) + h F(\\varrho)$ is computationally tractable in higher dimensions once the transport term is entropy-regularized and dualized. The optimality conditions collapse to the coupled equations $y \\odot (\\Gamma_k z) = \\varrho_{k-1}$ and $z \\odot (\\Gamma_k^\\top y) = \\xi_{k-1} \\odot z - \\beta\\epsilon/h$, with the updated density given by $\\varrho_k = z \\odot (\\Gamma_k^\\top y)$. The authors prove that the block-coordinate iteration for these equations is strictly contractive in the Thompson metric on the positive orthant, so Algorithm 1 has a unique fixed point and converges; combining this weight update with Euler-Maruyama moves of the particle locations gives the full density-propagation algorithm. They demonstrate the scheme on linear Gaussian systems, a bimodal nonlinear system, a McKean-Vlasov mean-field problem, a multiplicative-noise process via the Lamperti transform, and a six-state mixed conservative-dissipative satellite dynamics example.","pith_inferences":["If the split scheme is only consistent in the limit $h\\to 0$, practical accuracy will depend on the joint tuning of $h$, the regularization $\\epsilon$, and particle number $N$; the paper proves contraction of the inner fixed-point iteration but not consistency of the outer split, so the observed accuracy in examples may not hold uniformly across parameters.","Because the inner iteration is a positive-cone contraction, an adaptive time-stepping rule driven by the measured Thompson contraction rate could reduce the number of physical steps while preserving the fixed-point guarantee.","The same dual Sinkhorn structure can be applied to nonlinear filtering by replacing the free-energy term with a data-likelihood term, yielding an optimal-transport-based Bayesian update for the prior density propagated by this scheme.","The fixed-point equations match the structure of entropic optimal transport, so the transient densities computed here may coincide with Schrödinger-bridge interpolations in the small-noise limit, tying uncertainty propagation to stochastic control."],"forward_implications":["Transient joint PDFs for nonlinear Itô SDEs with gradient drift can be computed as weighted point clouds at a cost that scales with the number of particles, not with the state-space volume, so the method sidesteps the curse of dimensionality.","Each proximal weight update provably converges to a unique fixed point, giving a deterministic inner loop with a global contraction guarantee for every physical time step.","Systems outside JKO canonical form, including multiplicative noise and mixed conservative-dissipative drift, can be brought into the framework by coordinate transforms or modified transport cost functions, as shown for the CIR process and a satellite problem.","The McKean-Vlasov case with interaction potentials is handled by a semi-implicit free energy, so mean-field density evolution is available at the same proximal-scheme cost.","Reported runtimes near $10^{-6}$ seconds per proximal update at $N=400$, against a physical time step of $10^{-3}$ seconds, indicate the scheme is fast enough for repeated density propagation in filtering and control loops."],"supporting_citations":[{"why":"Establishes the JKO variational scheme that the paper's proximal recursion is built on.","marker":"[21]"},{"why":"Supplies the entropy regularization and dualization that turn the infinite-dimensional proximal problem into Sinkhorn-style fixed-point equations.","marker":"[32]"},{"why":"Introduces the fast Sinkhorn iteration for entropic optimal transport that motivates the block-coordinate solver in Algorithm 1.","marker":"[33]"},{"why":"Provides the metric gradient-flow background and convergence framework for proximal recursions on the space of probability measures.","marker":"[23]"},{"why":"Gives the time-varying coordinate transformations that make non-symmetric linear systems amenable to the JKO canonical form in the multivariate example.","marker":"[18]"},{"why":"Supplies the conservative-dissipative proximal scheme used to extend the algorithm to mixed drift dynamics.","marker":"[52]"},{"why":"Proves convergence of the semi-implicit scheme used for the McKean-Vlasov interaction case.","marker":"[47]"}],"fun_headline_variants":["Gradient flow propagates densities with no spatial grid","Time-discretized gradient flow beats curse of dimensionality","Proximal recursion solves density flow faster without grids","Wasserstein gradient flow: fast joint PDF propagation, no discretization","Block-coordinate gradient flow propagates densities in high dims"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that combining Euler-Maruyama point-location updates with the Sinkhorn-style weight updates reproduces the JKO scheme, and hence the Fokker-Planck-Kolmogorov PDE, as the time step $h$ tends to zero; the paper proves contraction of the inner fixed-point iteration but does not prove this outer consistency, and it also asserts rather than proves the Perron-Frobenius strict contraction of $\\Gamma_k$.","fun_headline_variants_meta":{"raw":{"variants":["Gradient flow propagates densities with no spatial grid","Time-discretized gradient flow beats curse of dimensionality","Proximal recursion solves density flow faster without grids","Wasserstein gradient flow: fast joint PDF propagation, no discretization","Block-coordinate gradient flow propagates densities in high dims"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000268,"raw_usage":{"total_tokens":1657,"prompt_tokens":1025,"completion_tokens":632,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":641,"completion_tokens_details":{"reasoning_tokens":550}},"tokens_in":641,"tokens_out":632,"duration_ms":5945,"temperature":1.0,"reasoning_tokens":550,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T15:49:38.391193+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the Ornstein-Uhlenbeck benchmark of Section IV-A1 with the same $N=400$ and $\\epsilon=5\\times 10^{-2}$ while shrinking $h$ from $10^{-2}$ to $10^{-6}$, and compare the proximal weighted mean and variance against the exact formulas $\\mu_0 \\exp(-at)$ and $(\\sigma_0^2 - 1/(a\\beta)) \\exp(-2at) + 1/(a\\beta)$; if the error does not systematically decrease as $h \\to 0$, the claimed consistency of the split scheme is false. Separately, compute the Thompson contraction factor of the linear map $z \\mapsto \\Gamma_k z$ for random positive $z$ with a randomly generated cost matrix $C_k$; if any eigenvalue of the associated projective map has modulus at least one, the asserted Perron-Frobenius strict contractivity fails.","supporting_citations":[{"cited_title":"The variational formulation of the Fokker–Planck equation,","cited_arxiv_id":null,"evidence_quote":"Establishes the JKO variational scheme that the paper's proximal recursion is built on."},{"cited_title":"Ambrosio, N","cited_arxiv_id":null,"evidence_quote":"Provides the metric gradient-flow background and convergence framework for proximal recursions on the space of probability measures."},{"cited_title":"Gradient Flows in Uncertainty Propagation and Filtering of Linear Gaussian Systems","cited_arxiv_id":"1704.00102","evidence_quote":"Gives the time-varying coordinate transformations that make non-symmetric linear systems amenable to the JKO canonical form in the multivariate example."},{"cited_title":"Conservative-dissipative approximation schemes for a generalized Kramers equation,","cited_arxiv_id":null,"evidence_quote":"Supplies the conservative-dissipative proximal scheme used to extend the algorithm to mixed drift dynamics."},{"cited_title":"On some nonlinear evolution systems which are pertur- bations of Wasserstein gradient flows,","cited_arxiv_id":null,"evidence_quote":"Proves convergence of the semi-implicit scheme used for the McKean-Vlasov interaction case."}],"review_version":1}