Deterministic and stochastic particle methods for the Fokker-Planck equation in S--formulation with application to point set registration

1.1 Point set registration in bounded domains

• Registration in bounded domains must respect geometric constraints, often requiring a bijection that maps the domain onto itself.Point-cloud alignment also arises in model order reduction and image processing.Source: §1.1 · p. 1
• Unknown point correspondences make direct joint estimation of correspondences and transformations a challenging combinatorial optimization problem.Probabilistic formulations have been used to avoid that optimization.Source: §1.1 · p. 1
• The Fokker–Planck-based registration procedure requires density estimation, numerical approximation of the equation, and particle transport.The paper focuses on the latter two tasks and also examines sensitivity to the density choice.Source: §1.1 · p. 2

1.2 Contributions and relation to previous work

• Particle methods focus computation where density is non-negligible and directly provide transport trajectories useful for registration-based model order reduction.The paper presents two complementary particle methods for this purpose.Source: §1.2 · p. 2
• The deterministic method discretizes the S-formulation with MLS, advects particles by target-driven drift, resolves diffusion on an evolving cloud, and imposes no-flux boundaries with ghost points.The evolving cloud can undergo strong compression and expansion.Source: §1.2 · p. 2
• The stochastic method simulates reflected Langevin trajectories and reconstructs deterministic osmotic paths using an adaptive-radius estimate of the evolving density’s score.The source distribution is represented by the initial particle ensemble, while the target enters through its score.Source: §1.2 · p. 2

2.1 Fokker-Planck equation for point set registration

• Under the theorem’s assumptions, the Fokker–Planck equation has a unique weak solution that converges in L1(Ω) to the target density.This convergence provides the density-level basis for the registration formulation.Source: §2.1 · p. 3
• The registration formulation separates density estimation, numerical approximation of the Fokker–Planck equation, and transport of particles.The transport map is sought so the pushforward density approaches the target as time tends to infinity.Source: §2.1 · pp. 3–4

2.2 S-formulation of the Fokker-Planck equation

• The S-formulation sets S = log ρ and rewrites the evolution as a mixed Lagrangian–Eulerian transport equation for particle positions and carried values of S.The diffusion term ∆S is linear in this formulation and is described as more amenable to implicit treatment.Source: §2.2 · p. 4
• The Neumann boundary condition can be implemented in the particle scheme with ghost particles mirrored across the boundary.In the absence of a driving term, the density equation reduces to the linear heat equation and the velocity is the osmotic velocity.Source: §2.2 · p. 4

2.3 Stochastic formulation of the Fokker-Planck equation

• The reflected stochastic process stays in the bounded domain, with boundary local time supplying the minimal inward displacement that prevents escape.The reflection term acts only at the boundary.Source: §2.3 · p. 5
• The stochastic formulation replaces deterministic Fokker–Planck solution with independent reflected Langevin trajectories, requiring the target score and initial samples rather than an analytic source density.The target density need only be known up to a positive multiplicative constant.Source: §2.3 · p. 5
• Figure 1 depicts grid points x_j, each point’s neighbor set C_i, and the distance r_max,i to its farthest neighbor.These quantities illustrate the local neighborhoods used in the particle discretization.Source: Figure 1 · p. 6

3 Meshfree method for the Fokker-Planck equation

• MLS particles follow the characteristic equation while the population equation is solved on the evolving cloud, concentrating computation where the density is nonzero without truncating the domain.Particles are initially distributed according to ρ0.Source: §3 · p. 6

3.1 Moving Least Squares

• Quadratic MLS gives second-order first derivatives and first-order second derivatives, while higher-order shifted-monomial fits can be ill-conditioned.The authors use this quadratic discretization throughout; LABFM showed boundary instabilities in some test cases.Source: §3.1 · p. 7

3.2.1 First-order scheme

• The first-order scheme uses implicit Euler for diffusion and explicit Euler for characteristics, preserving the equilibrium S = S∞ at the discrete level.Its time and MLS spatial discretizations are first-order, while the characteristic equation is second-order in space.Source: §3.2.1 · pp. 7–8

3.2.2 Second-order scheme

• ARS(1, 2, 2) yields a second-order scheme with explicit midpoint updates for characteristics and implicit midpoint updates for the population equation.The scheme uses an IMEX Runge-Kutta framework.Source: §3.2.2 · p. 8

3.3 Boundary conditions

• Edge-particle MLS neighborhoods become upwinded, so the spatial domain need not be truncated; ghost particles with reflected log-density impose reflective wall behavior.At a boundary, stencil symmetry makes the particle velocity's normal component zero.Source: §3.3 · pp. 8–9

3.4 Evaluation of the score

• For Gaussian-mixture target densities, the log-sum-exp trick rescales component contributions before exponentiation to prevent simultaneous underflow in the score calculation.The common positive scaling factor cancels between the score's numerator and denominator.Source: §3.4 · p. 9

4.1 Monte Carlo particle method

• Euler–Maruyama advances Langevin particles with weak order one and, for constant diffusion, strong order one; proposed positions inside circular obstacles are reflected outward.The method uses a consistent approximation to the Skorokhod problem for reflective obstacles.Source: §4.1 · p. 10

4.2 Osmotic paths

• A local mean-shift estimate from stochastic particles approximates the density score and points along the density gradient, with larger neighborhoods reducing statistical error but increasing bias.The estimate uses particles inside a circular neighborhood centered at the path location.Source: §4.2 · pp. 11–12

4.2.1 Optimal radius for bias-variance balancing of the score estimator

• The adaptive radius balances O(r^2) truncation bias against Monte Carlo error and is updated using estimated local density.Its proportionality constant is calibrated empirically and may need recalibration for each problem; shared-ensemble radius and score estimates produced nearly identical paths to split ensembles in the registration experiment.Source: §4.2.1 · pp. 12–13

5.1 Calibration of the adaptive-radius constant; 5.2.1 MLS

• The calibrated adaptive-radius constant is κ = 1.1675547312210792; diffusion tests show first-order Euler and second-order IMEX midpoint convergence.For the Gaussian diffusion test, quadratic MLS reproduces the log-density derivatives up to numerical error, so observed errors are dominated by temporal discretization.Source: §5.1 · p. 14, §5.2.1 · p. 15

5.2.2 Stochastic method; 5.6 Discussion

• Adaptive radii significantly reduce Monte Carlo mapping noise, but mean-shift score estimates remain sensitive to timestep and region size where density variations are large.Meshfree stability is not guaranteed for every particle cloud, so adaptive time stepping remains necessary.Source: §5.2.2 · pp. 15–16, §5.6 · p. 21

5.3.1 MLS

• In the advection–diffusion test, MLS methods again achieve first- and second-order convergence for density and the x-coordinate of final particle positions.The midpoint method is visually more accurate than the Euler method in the plotted paths.Source: §5.3.1 · p. 16, Figure 5 · p. 17

5.3.2 Stochastic method

• Monte Carlo paths agree well with paths obtained by integrating osmotic velocity using the exact density distribution.The small-region radii are adapted during the simulation to balance bias and Monte Carlo error.Source: §5.3.2 · p. 17

5.4.2 MLS

• At the cylinder obstacle, MLS particles split around its top and bottom, and the ghost-point method correctly enforces reflective boundary conditions before the groups merge.The final particle distribution is slightly distorted; the reported upper/lower asymmetry reflects asymmetry in the initial particle distribution.Source: §5.4.2 · p. 18

5.4.3 Stochastic method

• Near the obstacle, Monte Carlo osmotic paths can become discontinuous when high particle density pushes small regions back against the drift.The particle cloud stretches around the obstacle and eventually detaches from its boundary.Source: §5.4.3 · p. 18

5.5.1 MLS

• The KDE-based MLS paths differ qualitatively from the GMM-based paths, showing that density estimation significantly affects point-set registration results.The KDE approximation of the initial point cloud is also reported as less accurate than the GMM approximation.Source: §5.5.1 · p. 21

5.5.2 Stochastic method

• Monte Carlo osmotic paths follow the GMM-based MLS routes qualitatively, but their final center locations differ appreciably.The MLS simulations map the initial cloud toward a GMM- or KDE-represented target distribution rather than explicitly mapping each point to a target sample.Source: §5.5.2 · p. 21, §5.5.1 · p. 21

6 Conclusion and Outlook

• The Fokker–Planck formulation recasts registration as probability-density evolution toward a target, yielding a continuous flow that respects mass conservation and boundary constraints without explicit point matching.Source: §6 · p. 21
• Figures 10 and 11 depict the evolving particle cloud in a cylinder example and compare osmotic paths from MLS and Monte Carlo for different initial and target approximations.Figure 11 also plots MLS particles; for visibility, it shows only paths from small regions for Monte Carlo.Source: Figure 10 · p. 22, Figure 11 · p. 23
• The MLS method resolves diffusion on a deforming particle cloud and uses ghost points for reflective boundaries, capturing strong compression, expansion, and splitting.The reported effects arise from drift and domain geometry.Source: §6 · p. 24
• Monte Carlo reconstructs macroscopic osmotic paths from reflected Langevin trajectories using a local score approximation, trading deterministic accuracy for simplicity and scalability.The method is presented as especially attractive for high-dimensional problems when many particles and parallel computing resources are available.Source: §6 · p. 24
• Future work includes adaptive MLS stencils, time stepping and particle refinement; adaptive stochastic time stepping; density-estimator comparisons; and Schrödinger bridges in bounded domains [45].Source: §6 · p. 24

A Density estimation

• The experiments consider kernel density estimation and Gaussian mixture models for density estimation in two- and three-dimensional domains with N, M ≲ 10^3.Source: A Density estimation · p. 25
• The density-estimation method depends on the number of data points, dimension, and task of interest.Source: A Density estimation · p. 25
• The paper selects KDE and GMM because it targets moderate-size datasets in two- and three-dimensional domains.Source: A Density estimation · p. 25

A.1 Kernel density estimation

• For bounded domains, KDE normalization generally depends on sample locations because boundary truncation affects the kernels.Source: §A.1 · p. 25
• For sufficiently smooth densities, KDE converges in probability under λ → 0 and Nλ^d → ∞ as N → ∞; its asymptotically optimal bandwidth scales as λ ∝ N^(-1/(d+4)).Source: §A.1 · p. 25
• The experiments use leave-one-out and least-squares cross-validation to choose bandwidth, and these approaches yield similar experimental results.Source: §A.1 · p. 25
• Kernel choice affects approximation and stability: the kernel gradient should be Lipschitz continuous, the kernel positive within its support, and tail behavior matters for Fokker–Planck approximation.For rapidly decaying kernels, a unique nearest sample can dominate the far-field sum.Source: §A.1 · p. 25
• Common kernels show markedly different far-field behavior in their score function Sϕ.The passage gives Gaussian Sϕ(r) = −2r and begins an inverse multiquadric expression, but does not supply its full result.Source: §A.1 · pp. 25–26
• For the heavy-tailed kernel discussed, the nearest-sample approximation need not hold, although expansion of the full KDE shows its score vanishes in the far field.The paper leaves detailed comparisons of kernel functions and their effects on registration to future work.Source: §A.1 · p. 26

A.2 Gaussian mixture models

• A GMM represents a density compactly using K multivariate Gaussian components.Source: §A.2 · p. 26
• The mixture has weights summing to one, component means in R^d, and symmetric positive-definite covariance matrices.Source: §A.2 · p. 26
• The paper fits GMMs with MATLAB fitgmdist, selects K using AIC, and sets RegularizationValue to 10^-4 for numerical robustness.GMM variants with flexible far-field tails are noted, while the primary focus is simplifying approximation of the Fokker–Planck solution [50].Source: §A.2 · p. 26

• The initial and target distributions in Section 5.5 are each approximated by a GMM.Source: B GMM parameters · p. 26
• The listed component weights, means, and covariance matrices are rounded to four decimal places.Source: B GMM parameters · p. 26
• Simulations use full-precision parameters from the accompanying code; displayed mixture weights may not sum exactly to one and should be renormalized if used directly.Source: B GMM parameters · p. 26

• The initial distribution weights are w0 = [0.2507, 0.2559, 0.2700, 0.2235]⊤.Source: §B.1 · p. 26

• The target distribution weights are w∞ = [0.2424, 0.2439, 0.2347, 0.2789]⊤.Source: §B.2 · p. 26