A generalized alternating NGMRES method for PDE-constrained optimization problems governed by transport equations

TL;DR

The paper tackles accelerating PDE-constrained optimization problems governed by transport equations by casting the first-order update as a fixed-point iteration and injecting a generalized alternating NGMRES scheme. The proposed GA-NGMRES(w;σ,τ) periodically combines NGMRES updates with standard FP steps, yielding large runtime and iteration-time savings over state-of-the-art preconditioned gradient descent and Newton–Krylov methods. Extensive 2D numerical experiments across advection, mass-preserving, and incompressible transport models demonstrate speedups up to over 5x while maintaining accuracy, and show GA-NGMRES frequently outperforms Anderson acceleration variants. Limitations include a Matlab 2D prototype and sensitivity to small regularization α, with clear paths for 3D extension and deeper convergence analysis in future work.

Abstract

In this work, we propose a generalized alternating nonlinear generalized minimal residual method (GA-NGMRES) to accelerate first-order optimization schemes for PDE-constrained optimization problems governed by transport equations. We apply GA-NGMRES to a preconditioned first-order optimization scheme by interpreting the update rule as a fixed-point (FP) iteration. Our approach introduces a novel periodic mixing strategy that integrates NGMRES updates with FP steps. This new scheme improves efficiency in terms of both iteration count and runtime compared to the state-of-the-art. We include a comparison to first-order preconditioned gradient descent and preconditioned, inexact Gauss--Newton--Krylov methods. Since the proposed optimization scheme only relies on first-order derivative information, its implementation is straightforward. We evaluate performance as a function of hyperparameters, the mesh size, and the regularization parameter. We consider advection, incompressible flows, and mass-preserving transport (i.e., optimal transport-type problems) as PDE models. Stipulating adequate smoothness requirements based on variational regularization of the control variable ensures that the computed transport maps are diffeomorphic. Numerical experiments on real-world and synthetic problems highlight the robustness and effectiveness of the proposed method. Our approach yields runtimes that are up to 5x faster than state-of-the-art Newton--Krylov methods, without sacrificing accuracy. Additionally, our GA-NGMRES algorithm outperforms the well-known Anderson acceleration for the models and numerical approach considered in this work.

Figures (10)

• Figure 1: We show exemplary results for the baseline model ($H^2$ regularization; compressible velocity). The results correspond to run 8 in \ref{['t:h2-nk+rpgd']} (top row; NK) and run 12 in \ref{['t:nirep-300x300-na06-t0-na01-ngmres']} (bottom row; GA-NGMRES). Top row (from left to right): ($i$) the template image $m_0$ (image to be transported); $(ii)$ the reference image $m_1$, $(iii)$ the residual differences between $m_0$ and $m_1$ (white: small difference; black: large difference); and $(iv)$ the residual differences between the terminal state $m$ at $t=1$ and $m_1$ after solving for the optimal $v$. Bottom row (from left to right): $(i)$ final state $m$ at $t=1$; $(ii)$ optimal control variable $v$ (color indicates orientation); $(iii)$ determinant of the deformation gradient (the values are all positive, illustrating that the computed map $y$ is a diffeomorphism); and $(iv)$ computed mapping $y$.
• Figure 2: Convergence results for different optimization schemes. We plot the trend of the relative $\ell^\infty$-norm of the gradient $g^{(k)}$ and the mismatch (data fidelity term) as a function of the outer iteration count $k$. The results are for the nirep data. We show the plots for RPGD and our NK solver. For the NK method we consider three different preconditioners: the spectral (regularization) preconditioner (ireg); the two-level preconditioner (2lrpcsym), and the zero-velocity preconditioner (h0rpc).
• Figure 3: Convergence plots for GA-NGMRES. We consider the nirep dataset (native resolution: $300\times300$). We show the reduction of the relative norm of the gradient $g^{(k)}$ (top block) and the relative mismatch (bottom block) as a function of the iteration count $k$ for the hyperparameters $w$ and $p =(\sigma, \tau)$. The results shown here correspond to those reported in \ref{['t:nirep-300x300-na06-t0-na01-ngmres']} and \ref{['t:nirep-300x300-na06-t0-na01-ngmres_cont']}, respectively.
• Figure 4: We show exemplary results for modeling an incompressible transport map. The results correspond to run 7 in \ref{['t:cinfty-256x256-ngmres-stokes']} (GA-NGMRES). We consider an $H^3$-seminorm as a regularization model. Top row (from left to right): ($i$) the template image $m_0$ (image to be transported); $(ii)$ the reference image $m_1$, $(iii)$ the residual differences between $m_0$ and $m_1$ (white: small difference; black: large difference); and $(iv)$ the residual differences between the terminal state $m$ at $t=1$ and $m_1$ after solving for the optimal $v$. Bottom row (from left to right): $(i)$ final state $m$ at $t=1$; $(ii)$ optimal control variable $v$ (color indicates orientation); $(iii)$ determinant of the deformation gradient (the values are all positive, illustrating that the computed map $y$ is a diffeomorphism); and $(iv)$ computed mapping $y$. Notice that the determinant of the deformation gradient is equal to 1 to high accuracy ($(\text{min},\text{mean},\text{max},\text{std}) = (1.00e+00, 1.00e+00,1.00e+00,2.61e-13)$).
• Figure 5: We show exemplary results for mass preserving transport. The results correspond to run 13 in \ref{['t:otgauss-256x256-ngmres']} (GA-NGMRES). Top row (from left to right): ($i$) the template density $\pi_0$ (probability density to be transported); $(ii)$ the target density $\pi_1$, $(iii)$ the residual differences between $\pi_0$ and $\pi_1$ (white: small difference; black: large difference); and $(iv)$ the residual differences between the terminal state $\pi$ at $t=1$ and $\pi_1$ after solving for the optimal $v$. Bottom row (from left to right): $(i)$ final state $\pi$ at $t=1$; $(ii)$ optimal control variable $v$ (color indicates orientation); $(iii)$ determinant of the deformation gradient (the values are all positive, illustrating that the computed map $y$ is a diffeomorphism); and $(iv)$ computed mapping $y$.