Solving Fredholm Integral Equations of the Second Kind via Wasserstein Gradient Flows

1. Problem Statement

The paper addresses the numerical solution of Fredholm integral equations of the second kind, defined as:
$$\pi(x) = \phi(x) + \lambda \int_{\mathbb{R}^d} k(x, y)\pi(y) dy$$
where:

• $\pi$ is the unknown probability density (solution).
• $\phi$ is an observed forcing function.
• $k$ is a positive integral kernel.
• $\lambda \in \mathbb{R}$ is a scalar parameter.

Key Challenges:

• Unbounded Domains: Traditional methods (e.g., Nyström method, basis function discretization) struggle when the solution is defined over unbounded domains or in high dimensions ($d$ is large).
• Ill-Posedness: The problem is often ill-posed, particularly when the integral operator induced by $k$ is not invertible or poorly conditioned.
• Second Kind Complexity: Unlike first-kind equations, the unknown distribution $\pi$ appears on both sides of the equation (inside the integral and as the solution), leading to non-linear dependencies that complicate standard particle approximations.

2. Methodology

The authors propose a method based on Wasserstein Gradient Flows (WGF) and Mean-Field Particle Systems.

A. Variational Formulation

Instead of solving the integral equation directly, the authors formulate it as a minimization problem. They define a target functional $F_\alpha(\pi)$:
$$F_\alpha(\pi) = \text{KL}\left(\pi \,\Big|\, \phi + \lambda \int_{\mathbb{R}^d} k(\cdot, y)\pi(y)dy\right) + \alpha \text{KL}(\pi | \pi_0)$$

• KL Divergence: The first term measures the discrepancy between the candidate solution $\pi$ and the right-hand side of the integral equation. Minimizing this drives $\pi$ toward the solution.
• Regularization: The second term involves the Kullback-Leibler divergence between $\pi$ and a reference measure $\pi_0$, weighted by a regularization parameter $\alpha > 0$. This ensures the existence of a unique minimizer and stabilizes the solution, especially when the problem is ill-posed.
• Smoothing: A hyperparameter $\eta \ge 0$ is introduced in the denominator of the KL term to ensure the functional is well-defined and the drift is Lipschitz continuous.

B. Wasserstein Gradient Flow

The authors derive the gradient flow of $F_\alpha$ in the space of probability measures equipped with the Wasserstein metric ($W_2$).

• The flow is described by a partial differential equation (Fokker-Planck equation) for the time-evolving density $\pi_t$.
• This PDE corresponds to a McKean-Vlasov Stochastic Differential Equation (MKVSDE):
  $$dX_t = \left( \int_{\mathbb{R}^d} b_\eta(X_t, z, \pi_t) d\pi_t(z) - \alpha \nabla U(X_t) \right) dt + \sqrt{2(1+\alpha)} dB_t$$
  where $\pi_t$ is the law of $X_t$, $B_t$ is Brownian motion, and $b_\eta$ is a drift term involving nested expectations with respect to $\pi_t$.

Key Distinction from First-Kind Equations:
In first-kind equations (previous work by the authors), the drift depends on $\pi$ linearly. In this second-kind setting, the drift depends on $\pi$ non-linearly and via nested integrals (integrating the kernel against $\pi$ twice). This creates a more complex dependency structure.

C. Particle Approximation

Since the MKVSDE cannot be solved analytically, the authors approximate it using an interacting particle system of size $N$:
$$dX^{i,N}_t = \left( \frac{1}{N}\sum_{j=1}^N b_\eta(X^{i,N}_t, X^{j,N}_t, \pi^N_t) - \alpha \nabla U(X^{i,N}_t) \right) dt + \sqrt{2(1+\alpha)} dB^i_t$$
where $\pi^N_t = \frac{1}{N}\sum \delta_{X^{i,N}_t}$.

• Algorithm: The system is discretized in time using an Euler-Maruyama scheme.
• Computation: A naive implementation of the drift would cost $O(N^3)$. The authors optimize this to $O(N^2)$ by precomputing kernel integrals.
• Solution Reconstruction: The final solution $\hat{\pi}$ is reconstructed by plugging the particle approximation back into the original integral equation: $\hat{\pi}(x) = \phi(x) + \frac{\lambda}{N}\sum k(x, X^{i,N}_n)$.

3. Key Contributions & Theoretical Results

1. Existence and Uniqueness:

   • Proved that the regularized functional $F_\alpha$ admits a unique minimizer $\pi^*_{\alpha, \eta}$ in the space of probability measures with finite second moments, provided $\alpha > 0$.
   • Established conditions under which the minimizer of the regularized problem converges to the solution of the unregularized problem as $\alpha, \eta \to 0$.

2. Well-Posedness of the MKVSDE:

   • Derived the specific form of the drift for second-kind equations.
   • Proved existence and uniqueness of strong solutions to the associated MKVSDE under smoothness assumptions on $k$ and $\phi$.
   • Showed that the law of the particle system converges to the unique minimizer of the functional as $t \to \infty$.

3. Propagation of Chaos and Error Bounds:

   • Established a propagation of chaos result: the interacting particle system approximates the non-linear MKVSDE with an error scaling as $O(N^{-1/2})$.
   • Derived a global error bound combining particle approximation error and time-discretization error:
     $$\text{Error} \leq \frac{c_1(T)}{\sqrt{N}} + c_2(T)\sqrt{\gamma}$$
     where $\gamma$ is the time step.
   • Provided guidelines for optimal tuning: For a fixed computational budget $B \approx N^2/\gamma$, the optimal scaling is $N \propto B^{1/3}$ and $\gamma \propto N^{-1}$, yielding a total error decay of $O(N^{-1/2})$.

4. Handling Ill-Posedness:

   • Demonstrated that the regularization term $\alpha \text{KL}(\pi|\pi_0)$ is crucial. Without it (or with an improper reference measure), the algorithm can become unstable or fail to converge to a unique solution.

4. Experimental Results

The authors tested the method (named FE2kind-WGF) on several benchmarks:

• Gaussian Example:

  • Compared against Reversible Jump MCMC (RJ-MCMC), a standard method based on von Neumann series expansion.
  • Result: FE2kind-WGF remained stable and accurate even as $\lambda \to 1$ (where the von Neumann series diverges). RJ-MCMC failed in these regimes.
  • Showed that the choice of reference measure $\pi_0$ significantly impacts accuracy; a concentrated $\pi_0$ can compensate for diffusion effects in the SDE.

• Karhunen-Loève Expansions (Eigenfunctions):

  • Solved for eigenfunctions of squared exponential and exponential kernels.
  • Result: Compared against the Nyström method (deterministic discretization). FE2kind-WGF achieved up to 5x lower error for large $N$ and did not require fixed spatial discretization, making it suitable for unbounded domains.

• Gaussian Process State Space Models (GP-SSM):

  • Solved for the equilibrium distribution of a non-linear SSM where $\phi \equiv 0$ and $\lambda = 1$.
  • Result: The method successfully recovered the invariant measure where von Neumann-based methods fail (as they require a non-trivial forcing term).

5. Significance and Impact

• Extension to Second Kind: This work successfully extends the Wasserstein Gradient Flow framework from first-kind to second-kind integral equations, a significant theoretical leap due to the nested non-linearity.
• Unbounded Domains: The particle-based approach naturally handles unbounded domains without the need for artificial truncation or basis function expansions, which are bottlenecks in high dimensions.
• Robustness: The method is robust to ill-conditioning and cases where standard series expansions (von Neumann) diverge.
• Scalability: While the cost is $O(N^2)$, it avoids the "curse of dimensionality" associated with grid-based methods, making it a promising candidate for high-dimensional inverse problems in statistics, econometrics, and physics.
• Regularization Insight: The paper provides a rigorous justification for using KL-regularization with a reference measure to stabilize inverse problems, offering practical guidance on selecting $\pi_0$ and $\alpha$.

In summary, the paper presents a robust, theoretically grounded, and empirically superior method for solving a class of difficult integral equations that are otherwise intractable for standard numerical techniques.