PaperlayerView exact revision on arXiv

Source-linked AI summary

Picard-Based Acceleration of Newton Continuation for Mean Field Game PDE Systems

Mathieu Lauriere, Andrew Shi

arXiv:2609.11050v2math.NAcs.NA

TL;DR

Comparative evidence on the cost and robustness of Picard and Newton methods for MFG systems remains limited. This paper compares them across benchmarks and combines Picard iterations with Newton continuation; Picard is cheaper in favorable regimes, while Newton helps reach lower viscosities, and the hybrid reduces computational cost in a two-dimensional double-well problem.

  • Problem

    The comparative cost and robustness of Picard and Newton methods across MFG regimes remain insufficiently explored and documented.

    Source: §1 · p. 2

  • Method

    The paper compares Picard and Newton for the same finite-difference discretization and combines inexpensive Picard iterations with Newton continuation in a hybrid method.

    Source: §1 · p. 2, §5 · p. 26

  • Results

    Picard provides substantial savings in favorable regimes, while Newton continuation is more robust at lower viscosities; the one-dimensional discounted benchmark reports a roughly 137-fold speedup for Picard over the coupled space-time Newton solve.

    Source: §1 · p. 2, §4.1 · p. 14

  • Takeaways & Limitations

    The hybrid replaces expensive high-viscosity Newton solves with Picard solves, then transfers the last converged Picard solution to Newton to reach a smaller target viscosity on the same discretization.

    Source: §5 · p. 26

    A cost-based handoff criterion remains future work, and the convergence and cost of an intermediate Picard/Newton variant remain unassessed.

    Source: §5 · p. 27

Abstract

from arXiv · show

We develop a method that uses Picard iterations to accelerate Newton continuation for the semi-implicit finite-difference discretization of forward-backward partial differential equation (PDE) systems arising in mean field games (MFGs). We first investigate the computational properties of the Picard and Newton methods, which are widely used separately in the MFG literature but whose comparative cost and robustness across different regimes remain insufficiently explored and documented. The Picard method uses an outer fixed-point iteration that alternates a forward Fokker-Planck solve and a backward Hamilton-Jacobi-Bellman solve. The Newton method instead applies Newton's method directly to the coupled nonlinear space-time system. Across one- and two-dimensional MFG benchmarks, Picard offers substantial computational savings in favorable regimes, but may fail at sufficiently low viscosity or require strong damping under temporal shocks, making Newton continuation preferable. With parameter continuation in the viscosity parameter, the Newton method is more robust in these regimes, at the cost of larger coupled linear systems. We relate these trade-offs to the residuals, Jacobian blocks, and sparsity structures produced by separable, local nonseparable, and nonlocal Hamiltonians. We then demonstrate how to combine inexpensive Picard iterations with Newton continuation in a hybrid method to reduce the total computational cost for a two- dimensional double-well problem.

1. Introduction

The paper compares Picard and Newton solvers for the same finite-difference MFG discretization, then uses their complementary strengths to motivate a hybrid method. Its contributions include an implementation-level analysis of solver costs and a hybrid strategy combining Picard iterations with Newton continuation.

1. Introduction

  • The solver structures are not new; prior work includes full space-time Newton solvers and alternating forward-backward fixed-point iterations.Related comparisons include runtime on a one-dimensional benchmark and Newton linearization with an Armijo line search.Source: §1 · p. 2
  • Picard can provide inexpensive initial solves, while Newton continuation supplies the robustness needed to reach lower viscosities across one- and two-dimensional benchmarks.The comparison is for the same semi-implicit finite-difference discretization.Source: §1 · p. 2
  • An implementation-level account relates solver costs to residuals, Jacobian blocks, and sparsity structures for separable, local nonseparable, and nonlocal Hamiltonians.Source: §1 · p. 2
  • The hybrid continuation method uses inexpensive Picard iterations to reduce coupled Newton solves and total computational cost.The paper demonstrates the strategy on a two-dimensional double-well problem.Source: §1 · p. 2

2. Governing equations

The paper sets out the finite-horizon, discounted, and ergodic MFG systems used in its numerical experiments, along with their connections. It specifies assumptions on the feedback control and relates stationary reference problems to finite-horizon computations.

2.1. MFG PDE systems for different optimal control formulations

  • The numerical experiments use finite-horizon, infinite-horizon discounted, and ergodic MFG systems, assuming a unique maximizing feedback control and well-defined derivatives required for Newton linearizations.Source: §2.1 · p. 2

2.1.1. Finite-horizon optimal control

  • The finite-horizon control problem yields a coupled HJB-FP system, with initial population density and terminal cost conditions.The FP equation describes density evolution under the optimal drift.Source: §2.1.1 · p. 3

2.1.2. Infinite-horizon discounted MFG system

  • The discounted system has a stationary equilibrium formulation, but computation uses a finite-horizon forward-backward version with a computational horizon and arbitrary terminal cost.The stationary formulation provides an analytical benchmark in special cases.Source: §2.1.2 · p. 3, §2.1.3 · p. 4

2.1.3. Infinite-horizon ergodic MFG system

  • The ergodic formulation characterizes a long-run average-cost equilibrium; the value function is defined up to an additive constant unless boundary data fix it.The experiments also use ergodic references for validation.Source: §2.1.3 · pp. 3–4

2.2. Connections between the formulations

  • Under suitable assumptions, the discounted value function approaches the ergodic cost after scaling, while long-horizon finite-horizon trajectories can approximate ergodic states away from boundary layers.These links serve as numerical validation devices, including midpoint comparisons with stationary references.Source: §2.2 · p. 4

3. Numerical methods

This section presents the semi-implicit finite-difference discretization of the forward-backward MFG system and the Picard and Newton solvers, including viscosity continuation and their linear-algebra backends. It also describes a hybrid strategy that uses Picard iterations before handing off to Newton continuation when Picard fails.

3.1. Finite difference discretization

  • The semi-implicit scheme uses a Godunov isotropic upwind Hamiltonian and discretizes the FP equation as a flux divergence, with each equation implicit in its own unknown.The scheme’s four coordinate-direction flux pairs form the discrete divergence; the HJB solve treats the next mass as known, and the FP solve treats the current value function as known.Source: §3.1 · pp. 5–6

3.2. Newton and Picard methods for the forward-backward problem

  • Picard alternates backward nonlinear HJB solves and forward linearized FP solves, whereas Newton solves the fully coupled space-time system in a larger linear system.Picard uses Newton’s method within each timestep’s HJB solve and evaluates density-dependent transport coefficients at the previous Picard iterate for non-separable Hamiltonians; Newton has local quadratic convergence with a suitable initial guess.Source: §3.2 · pp. 6, 8, 9, Table 1 · p. 9

3.3.1. Continuation and initialization

  • Viscosity continuation improves Newton’s robustness by decreasing viscosity from an easier problem and backtracking attempted decreases when a solve fails.The method uses each converged solution to initialize the next, starting from a larger viscosity where the naive initialization is sufficient.Source: §3.3.1 · p. 9

3.3.2. Linear algebra backend

  • For a two-dimensional grid with 200 points in each spatial direction and 100 timesteps, Newton requires solving a coupled space-time Jacobian with 8,080,000 degrees of freedom, making direct factorization intractable due to fill-in.Source: §3.3.2 · p. 9
  • The two-dimensional experiments use a common direct sparse backend, while nonlocal congestion creates dense Jacobian cross-couplings that make explicit assembly intractable at relevant grid sizes.The reported experiments use NT = 20 in two dimensions and NT = 10 for hybrid timings, with spsolve/SuperLU.Source: §3.3.2 · p. 10

3.4. Hybrid method: Picard-based acceleration of Newton continuation

  • Picard descends in viscosity on the target grid until convergence fails; Newton then continues the descent from the last converged Picard state.The strategy remains on the target grid because the double-well sweep suggests further refinement does not lower Picard’s minimum reachable viscosity once the grid is sufficiently fine.Source: §3.4 · p. 10
  • Picard attempts use independent diffusion-based initializations and retry failed viscosity steps with progressively reduced damping, down to θ = 0.05.The experiments use the fixed viscosity schedule 1, 0.3, 0.1, 0.03, 0.01, 0.003, 0.001, 0.0003, 0.0001; damping starts at θ = 0.8 and decreases by 0.1 after failure.Source: §3.4 · p. 10
  • After Picard exhausts its damping schedule, the last converged state defines νlimit and Newton re-solves there before continuing to the target viscosity.The re-solve reduces the coupled global residual because Picard stops using relative error between successive iterates; empirically, the converged Picard state lies within Newton’s basin of attraction at the same discretization.Source: §3.4 · p. 10

4. Numerical experiments

These numerical experiments assess Picard and Newton behavior across stationary-reference, vanishing-discount, discounted, nonlocal-congestion, and 2D double-well benchmarks, including how nonlocal structure, temporal shocks, and low viscosity affect feasibility. They also examine a Picard-to-Newton hybrid that switches at Picard’s lower-viscosity limit and reduces cumulative runtime versus pure Newton.

4.1. A discounted problem

  • Direct Newton attempts to reach nu=0.5 diverged until the current viscosity was sufficiently close; backtracking required 35 continuation attempts, of which 11 converged.The reported schedule used reduction and backtracking factors of 0.5; each converged viscosity value required 4–6 Newton updates.Source: §4.1 · p. 12, Figure 3 · p. 13
  • In the one-dimensional benchmark, Picard’s sequential temporal marching and reduced-order spatial Jacobians yielded a roughly 137-fold speedup over the coupled space-time Newton solve (17-fold excluding failed backtracks).Source: §4.1 · p. 14

4.2. An ergodic problem

  • The ergodic validation compared a finite-horizon solution with an analytic stationary profile; constant terminal data instead converged to a different solution, even when the horizon was extended to T=40.The experiment anchored the terminal value to the analytic profile because the coupling is not monotone.Source: §4.2 · p. 14, Figure 5 · p. 15
  • O(rho) scaling in ||rho u_rho(., T/2) - lambda_{Omega,h}||_infinity numerically supports the discrete vanishing-discount relation, with lambda_{Omega,h}=2.616867 as the grid benchmark.This is a numerical verification for the chosen grid and horizon, not an independent convergence theorem; the maximum error is attained on the boundary, while interior values also decrease at first order.Source: §4.2 · p. 15, Figure 6 · p. 16

4.3. A problem with nonlocal congestion effects

  • At Nx=Ny=200, assembled Newton required over 250 GB for the AU,M quadrant in the 2D nonlocal case, versus approximately 32 MB of dense-block overhead in 1D, favoring Picard in 2D.This limitation applies to the assembled backend; a matrix-free method could apply the convolution-based dense blocks without storing them explicitly.Source: §4.3 · pp. 17–18
  • In 2D, nonlocal density stays more concentrated than local throughout the horizon: averaging lowers congestion across the initial blob, so the square patch disperses more slowly, bifurcation is not pronounced until t≈0.75, and terminal well peaks exceed local 2D.Source: §4.3 · pp. 18–19

4.4. A “traffic-light” problem: an example of a temporal “shock”

  • For the temporal-shock traffic-light problem, Picard fails at θ = 0.3 but converges in 141 iterations at θ = 0.1; Newton takes 7.55 s to reach ν = 0.001 versus 25.74 s for Picard at ν = 0.02.Newton continuation robustly traverses the shock from higher viscosity, whereas Picard convergence remains brittle; the wall-clock comparison uses the same grid but different viscosities.Source: §4.4 · p. 20, Figure 10 · p. 20

4.5. A problem with a double-well potential

  • In the asymmetric double-well problem, density bifurcates toward both targets but preferentially concentrates in the deeper left well.The tilt makes the left well lower in terminal cost, and congestion encourages agents to spread out.Source: Figure 12 · p. 22, §4.5 · p. 22

4.5.1. Picard contractivity diagnostic

  • Picard loses local contraction as ν approaches zero when the spectral radius of its linearized iteration reaches or exceeds 1; this diagnostic is computed in 1D, not directly in 2D.For damped Picard, local contraction requires ρ(Pθ) < 1; the 2D failure is expected to share the mechanism, but its space-time Jacobian spectral radius was not computed.Source: §4.5.1 · pp. 22–23, Figure 13 · p. 23

4.5.2. Two-dimensional cost and hybrid continuation

  • 312.43 seconds were required to factor the 55,000-DOF Newton system, whose factors had 128,928,646 nonzeros and a 283.36× fill-in ratio.The representative system used Nh = 50, NT = 10, and ν = 0.5; its matrix had 455,000 nonzero entries.Source: §4.5.2 · p. 25
  • The smallest successful Picard viscosity rises from 10^-4 at Nh = 16 to 10^-2 at Nh = 64 and 128, motivating Newton continuation below the apparent floor.At Nh = 32, the smallest successful viscosity is 3 × 10^-3; the reported floor is under the prescribed damping schedule and stopping criteria.Source: §4.5.2 · p. 25
  • At Nh = 128 and νtarget = 10^-4, pure Picard fails and pure Newton costs 46,410 s; the hybrid uses Picard before switching to Newton continuation.The hybrid’s handoff follows Picard failure at the next viscosity; the supplied passages do not report its total runtime.Source: §4.5.2 · pp. 23, 26, Table 6 · p. 26
  • Roughly 62% of pure Newton’s total cost goes to high-ν continuation steps, which the hybrid replaces with cheaper Picard iterations to reduce cumulative runtime.The hybrid switches to Newton at Picard’s lower-viscosity limit, reaching the target where pure Picard fails.Source: §4.5.2 · p. 26

5. Conclusion and future work

The hybrid uses inexpensive Picard solves during initial viscosity descent, then transfers the last converged state to Newton continuation to reach a smaller target viscosity. Its benefits depend on the computational regime and reported linear-solver backend, with several handoff and globalization questions left for future study.

5. Conclusion and future work

  • The hybrid reduces total computational cost by replacing high-viscosity Newton solves with Picard solves, then handing off to Newton to reach a smaller target viscosity.It transfers the last converged Picard solution to Newton on the same discretization.Source: §5 · p. 26
  • Dense coupling makes full Newton Jacobian assembly impractical for nonlocal congestion at tested grid sizes, while strong damping erodes Picard’s advantage in the traffic-light problem.There, the recorded Newton run reaches a lower viscosity in less time than the recorded Picard run.Source: §5 · p. 26
  • The reported computational gains use a common direct sparse backend, so alternative linear solvers may change phase costs and the appropriate handoff point.The authors propose assessing block strategies for Newton and incomplete-LU preconditioning for Picard, distinguishing the number of linear solves from their individual cost.Source: §5 · p. 26
  • A cost-based handoff criterion and an intermediate method with outer HJB–FP decoupling but globally solved HJB remain to be studied.The intermediate variant’s convergence and computational cost have not yet been assessed.Source: §5 · p. 27
  • The authors propose studying viscosity continuation alongside Armijo line-search globalization, which controls Newton steps at fixed parameters, and developing computable convergence and schedule criteria.These are identified as next steps toward a more systematic hybrid method.Source: §5 · p. 27

Appendix A. Pseudocode

The pseudocode specifies the coupled Newton solve, damped Picard sweeps, viscosity continuation, and a hybrid handoff from Picard to Newton.

Appendix A. Pseudocode

  • Newton iteratively evaluates residuals, constructs the Jacobian, solves for a coupled update, and returns the converged value and density.The update solves Jφ(W^k)δW^k = −φ(W^k).Source: Appendix A · p. 30
  • Picard alternates a backward HJB solve using Newton at each time step with a forward Fokker–Planck solve, then damps the value and density updates.The algorithm initializes the value at zero and the density from the initial distribution.Source: Appendix A · p. 31
  • In each Picard sweep, the forward Fokker–Planck step uses the updated value and evaluates density-dependent coefficients at the previous density iterate.The density is marched forward from the initial condition before its damped update.Source: Appendix A · p. 31
  • Newton continuation reuses the previous solution as the next viscosity-level initial guess and backtracks the viscosity step after failed convergence.It aborts if the step size falls below tolerance.Source: Appendix A · p. 32
  • The hybrid descends in viscosity with Picard, retries failed levels with stronger damping, then hands off the last converged state to Newton for the remaining descent.Newton re-solves at the handoff viscosity before continuing toward the target.Source: Appendix A · p. 33
Loading 2609.11050v2…