A Scalable Interior-Point Gauss-Newton Method for PDE-Constrained Optimization with Bound Constraints

TL;DR

The parallel scalability of the preconditioner, achieved via algebraic multigrid component solvers when applicable, and the aforementioned algorithmic scalability permits a parallel scalable means to compute solutions of a large class of PDE‐ and bound‐constrained problems.

Abstract

We present a scalable approach to solve a class of elliptic partial differential equation (PDE)-constrained optimization problems with bound constraints. This approach utilizes a robust full-space interior-point (IP)-Gauss-Newton optimization method. To cope with the poorly-conditioned IP-Gauss-Newton saddle-point linear systems that need to be solved, once per optimization step, we propose two spectrally related preconditioners. These preconditioners leverage the limited informativeness of data in regularized PDE-constrained optimization problems. A block Gauss-Seidel preconditioner is proposed for the GMRES-based solution of the IP-Gauss-Newton linear systems. It is shown, for a large-class of PDE- and bound-constrained optimization problems, that the spectrum of the block Gauss-Seidel preconditioned IP-Gauss-Newton matrix is asymptotically independent of discretization and is not impacted by the ill-conditioning that notoriously plagues interior-point methods. We propose a regularization and log-barrier Hessian preconditioner for the preconditioned conjugate gradient (PCG)-based solution of the related IP-Gauss-Newton-Schur complement linear systems. The scalability of the approach is demonstrated on an example problem with bound and nonlinear elliptic PDE constraints. The numerical solution of the optimization problem is shown to require a discretization independent number of IP-Gauss-Newton linear solves. Furthermore, the linear systems are solved in a discretization and IP ill-conditioning independent number of preconditioned Krylov subspace iterations. The parallel scalability of preconditioner and linear system matrix applies, achieved with algebraic multigrid based solvers, and the aforementioned algorithmic scalability permits a parallel scalable means to compute solutions of a large class of PDE- and bound-constrained problems.

Figures (4)

• Figure 1: Top row: noisy state observations $u_{d,\zeta}$ (left), state reconstruction $u^{\star}$ (middle), computed adjoint $\lambda^{\star}$ (right). Bottom row: true parameter $\rho_{\text{true}}$ (left), parameter reconstruction $\rho^{\star}$ (middle), computed bound constraint Lagrange multiplier $z_{\ell}^{\star}$ (right). The dimension of each of the discretized fields is $124\,609$.
• Figure 2: Strong scaling for the IPM framework on the example problem \ref{['sec:problemdescription']} using the preconditioned reduced-space IP-Gauss-Newton CG solver (dashed lines) described in Section \ref{['subsec:CGpreconditioner']} and the Gauss-Seidel preconditioned IP-Gauss-Newton GMRES solver (solid lines) described in Section \ref{['subsec:GMRESpreconditioner']}.
• Figure 3: The number of $\boldsymbol{W}_{\! \boldsymbol{\rho}, \boldsymbol{\rho}}^{\mu}$ preconditioned CG, block Gauss-Seidel ($\boldsymbol{\tilde{A}}$) preconditioned GMRES and central null ($\boldsymbol{\tilde{A}}_{\text{cen}}$) preconditioned GMRES iterations to solve the IP-Gauss-Newton linear system and the log-barrier $\mu$ at each step of the IP-Gauss-Newton method. The results are for the problem described in Section \ref{['sec:problemdescription']} where the dimension of the discretized fields is equal to $2\,362\,369$. The number of Krylov subspace iterations has some initial variation, but quickly becomes relatively constant as the log-barrier parameter $\mu$ varies over approximately four orders of magnitude, where the central null preconditioner has a significantly higher iteration count. This provides numerical evidence that the preconditioned matrices have spectral properties that are asymptotically independent of $\mu$.
• Figure 4: Left: spatial structure of a random sample $\zeta$. Right: Seminorm of discrepancy $(u^{\star}-u_{d,\zeta})$ and seminorm of noise $\zeta$ as functions of the regularization parameters $\gamma_{1}=\gamma_{2}$, plot obtained in order to apply the Morozov discrepancy principle for choosing numerical values of the regularization parameters with a $5\%$ noise level $\sigma_{\zeta}$.

Theorems & Definitions (12)

• proposition 1
• proof
• proposition 2
• proof
• proposition 3
• proof
• proposition 4
• proof
• proposition 5
• proof