# Presentation: 352575

## Towards Optimal Solvers for PDE-Constrained Optimization

-  Volkan Akcelik    Omar Ghattas
- Carnegie Mellon University
- George Biros
- New York University

- TOPS Winter Meeting
- Livermore, CA     January 25-26, 2002

## Simulation vs. Optimization

- PDE model
- Simulation (forward) problem
- Given “data” *x* (e.g. material coefficients, domain or boundary sources, geometry), find state variables *u *(velocity, stress, temperature, magnetic or electric field, displacement, etc.)
- Optimization (inverse) problem
- Given desired goal involving *u*, find *x*
- Optimal design, optimal control, or parameter estimation

## PDE-Constrained Optimization 

## First-order Optimality Conditions

## Examples 

- Optimal design
- Shape optimization of viscous flow system
-  *N**d* *= O(1)*
- Optimal control
- Boundary flow control
*N**d** = O(N**u**2/3**)*
- Parameter estimation
- Heterogeneous inverse wave propagation
*N**d** = O(N**u**)*

## Optimal Design: Artificial Heart Shape Optimization 1

- Goal: Find pump geometry that minimizes blood damage

-  (J. Antaki & G. Burgreen, Univ Pittsburgh Medical Center)

- G. Burgreen

## Optimal Design: Artificial Heart Shape Optimization 2 

- Burgreen & Antaki, 1996

## Optimal Control:Boundary Flow Control 1

- MHD flow control using boundary magnetic field

## Optimal Control: Boundary Flow Control 2

- Goal: find boundary suction/injection that minimizes energy dissipation in a viscous incompressible flow

## Optimal Control: Viscous Boundary Flow Control 2

- (G. Biros, 2000)

## Optimal Control: Viscous Boundary Flow Control 3

- No control

- Optimal control

## Parameter Estimation: Seismic Inversion 1

** **Forward Problem**: **Given soil material and earthquake source parameters, find earthquake ground motion.
** **Inverse Problem**: **Given earthquake observations, estimate material and source parameters.
- Part of Quake KDI Project at CMU (J. Bielak, D. O’Hallaron), Berkeley (J. Shewchuk), and SDSU (S. Day, H. Magistrale)

## Parameter Estimation: Seismic Inversion 2

## Parameter Estimation: Seismic Inversion 3

- Simulation of 1994 Northridge Earthquake

## Parameter Estimation: Seismic Inversion 4

- Surface Visualization

## Lagrange-Newton-Krylov-Schur method

- Consider discrete form of problem:

## Optimality Conditions

- Lagrangian function:
- First order optimality conditions:
- Define:

## Newton’s Method for Optimality Conditions (SQP, Newton-Lagrange)

- A Newton step:

- Two possibilities for solution:
- Full space method: solve KKT system with a Krylov method—need good preconditioner
- Reduced space method: eliminate state variables and Lagrange multipliers

- (KKT system)

## Quasi-Newton Reduced SQP

- LU factorization of (permuted) KKT matrix
- Problems:
- Reduced Hessian **W***z* needs *m* (number of decision variables) PDE solves per optimization iteration
- 2nd derivatives needed
- Quasi-Newton RSQP:
- Replace reduced Hessian by quasi-Newton approximation
- Discard other 2nd derivative terms
- Properties
- Just 2 PDE solves required per optimization iteration
- 2-step superlinear convergence

## Use Schur-based approximate LU factorization as KKT preconditioner

Since preconditioner, can replace PDE solves by actions of PDE preconditioner
Properties:
No PDE solves per iteration (2 preconditioner applications required per inner Krylov iteration)
Newton convergence for outer iteration
Inspiration:
Schur complement domain decomposition preconditioners (Keyes & Gropp, 1987)

- Since preconditioner, can replace PDE solves by actions of PDE preconditioner
- Properties:
- No PDE solves per iteration (2 preconditioner applications required per inner Krylov iteration)
- Newton convergence for outer iteration
- Inspiration:
- Schur complement domain decomposition preconditioners (Keyes & Gropp, 1987)

- Approximate QN-RSQP Preconditioner

## Effectiveness of Schur Preconditioner

- Preconditioned KKT system (w/exact **A**)
- Schur preconditioner results in:
- Reduced condition number
- Clustering of eigenvalues
- But not a normal matrix—convergence rate bounds depend on condition number of eigenvector matrix
- In practice, works very well

## How to Precondition Reduced Hessian?

- Requirements
- Cannot compute reduced Hessian, only its (approximate) action on a vector
- Possibilities
- Quasi-Newton formula (limited memory for large m)
- Fixed small number of stationary iterations
- Sparse approximate inverse preconditioner (SPAI)
- Discrete approximation of underlying continuous Wz
- Sherman-Morrison-Woodbury, if reduced Hessian is of form “compact operator + SPD”

## 2-step Stationary Iterations as Reduced Hessian Preconditioner 

* *Properties
-  Same complexity as CG but a constant preconditioner
-  Needs only Matvecs (here with approximate reduced Hessian)
-  Needs estimates of extremal eigenvalues
-  Use to initialize L-BFGS preconditioner (effectively H0)

## Inexact Solution of KKT System

- Solve KKT system iteratively to tolerance dependent on
- ** x Lagrangian gradient
- Eisenstat & Walker inexact Newton theory applies for reduction of norm of Lagrangian gradient
- But will it satisfy merit function criteria?
- Can show that in vicinity of local minimum (for augmented Lagrangian merit function):
- If ** small enough, Armijo condition satisfied with unit steplengths
- If ** is *O*(Lagrangian gradient), quadratic convergence preserved

## Putting it all together: Lagrange-Newton-Krylov-Schur Method

- Continuation loop
- Optimization iteration (*Lagrange-Newton*)
- Estimate extremal eigenvalues of (approximate) reduced Hessian using Lanczos (retreat to QN-RSQP if negative)
- Inexact KKT solution via symmetric QMR (*Krylov*)
- Quasi-Newton RSQP preconditioner (*Schur*)
- 2-step stationary iterations+ L-BFGS approximation of inverse reduced Hessian
- PDE solve replaced by PDE preconditioner
- Backtracking line search on augmented Lagrangian or *l*1 merit function
- If no sufficient descent use QN-RSQP
- Compute derivatives, objective function, and residuals
- Update solution, tighten tolerances

## Application: Optimal Boundary Control of Viscous Flow

- Goal: find boundary suction/injection that minimizes energy dissipation in a viscous incompressible flow

## Fixed-size efficiency

- Veltisto + PETSc implementation

- algorithmic

- parallel

- overall

- (n=120,000)

## Isogranular algorithmic efficiency 

## Veltisto: A Library for PDE-Constrained Optimization  

- PETSc extension
- Object oriented
- Nonlinearly constrained optimization
- Targets steady PDE constraints
- Supports domain-based parallelism
- Matrix-free implementation
- In use at Sandia-Albuquerque

## Veltisto Features  

- Optimizers: LNKS, LNKS-IP, QN-RSQP, N-RSQP
- Preconditioners: LM-BFGS, SR1, Broyden, 2-step stationary
- Merit functions: augmented Lagrangian, L1, L1+second order correction
- Symmetric QMR (allows indefinite preconditioning)

## Veltisto Example 

## Veltisto Example 

## Conclusions	

- Schur-complementing the optimality system permits capitalizing on investment in PDE solvers
- Lagrange-Newton-Krylov-Schur methods can exhibit mesh-independent outer iterations for mesh-based control spaces
- Inner iteration potentially mesh independent
- Depends on effective PDE and reduced Hessian preconditioners
- Very good scalability
- >10x speedup over quasi-Newton
- Optimization in small multiple of simulation cost

## Future Work		

- Approximate block elimination preconditioners
- Other reduced Hessian preconditioners
- Direct preconditioning of KKT system
- DD preconditioners directly on KKT system
- Multigrid on KKT system (Barry)
- Multigrid preconditioner on block diagonal approximation of KKT system
- Reformulation of KKT system to make it SPD
- FOSLS?
- Augmented Lagrangian?
- Time dependence and inequalities remain a challenge

## Acoustic inversion

## First order optimality conditions

## First order optimality conditions, cont.

## Summary: first order optimality conditions (Karush-Kuhn-Tucker)

## Newton solution of the KKT conditions

## Schur complement method

-  Schur complement is self adjoint & positive definite near minimum
-  Action of Schur complement on q requires 1 forward & 1 adjoint wave prop

## Slide 41

## Gauss-Newton-Schur method

-  GN Schur complement is self adjoint & positive definite everywhere
-  Quadratic convergence for good fit problems, linear otherwise
-  Action of Schur complement on *q* requires 1 forward & 1 adjoint wave prop

## Spectrum of GN Schur complement (just least squares term; no regularization)

- Spectrum of GN Schur complement (just least squares term; no regularization)

- rough
- eigenvectors

- smooth eigenvectors

- 400 eigenvalues

- first 12 eigenvalues

## -(Gauss)-Newton-Schur-CG Solver

- Multiscale continuation over *u*, **, *q* grids and source frequency
- (Gauss)-Newton nonlinear iteration
- Conjugate gradient solution of Schur-decomposed linear system
- Currently no preconditioner

## 2D/3D Acoustic Wave Equation Example

- Piecewise bi(tri)linear approximation of pressure *u*, adjoint **, and square of wave speed q
- Explicit time integration
- Up to 512x512 and 64x64x64 grids, 2000 time steps
- 50 receivers, 1 harmonic source
- Homogeneous initial guess
- PETSc (_[www.petsc.anl.gov](http://www.petsc.anl.gov/)_) implementation
- Performance (Tikhonov regularization)
- 6 hours on 64 processors of T3E-900 at PSC
- Number of inner conjugate gradient iterations 10~20
- Number of outer Newton iterations 10~20

## Effect of different types of regularization

## Total variation regularization, finer grids

**Notes:**
 Hi Omar,
   I am sending you some images I obtained using TV regularization.
   I have solved three different images. Results look very nice. The
largest problem I have solved has 513x513 grid.
   I am sending four ps files. One of them is quite big. I hope you will
not have any problem opening it.
   19_1.ps has 9 figures. The first 7 of them are the result of 2x2, 3x3,
5x5, 9x9, 17x17, 33x33, 65x65 optimization grid with 65x65 wave
propagation grid. The 8th one uses 129x129 optimization grid and wave
propagation grid. The last one is the target. As you can see, one can not
tell the difference between the target and the solution. I have used 63
receivers and sources. The source is the Ricker wave. 65x65 grid has a
source with 0.2L wavelength, and 129x129 grid has source with 0.1L
wavelength. All solutions satisfy the convergence criteria.
   19_2.ps has 9 figures. Everything is similar to 19_1.ps. The 8th figure
-which is the final solution- has 513x513 grid.
   19_3.ps shows two sections showing velocity profiles from 8th and 9th
figures of 19_2.ps. Unlike 19_1.ps this time you can tell the difference
between the solution and the target.
   19_4.ps has 4 figures. The third figure has 257x257 grid. The solution
is very close to the target.
   These examples illustrates the power of total variation. I must
admit that this method is much petter that second derivative operator. The
main problem is the selection of optimization parameter. I have used trial
and error for the selection.
   Volkan

## Effect of different regularization weights(Tikhonov regularization)

## Multilevel vs. single level continuation(Tikhononv regularization)

- multilevel solution

## Effect of data noise on solution

## 3D Inversion

## Performance: linear and nonlinear iterations

## Performance, objective function value

## Objective function iteration/grid history                                           

## Conclusions	

- Multilevel continuation forces successive velocity estimates to remain within basin of attraction of global minimum
- Outer iterations are mesh-independent for Newton method
- Inner iterations weakly mesh-dependent for Tikhonov regularization (but scale poorly for TV regularization)
- Preconditioner in progress
- Limited memory quasi-Newton for least squares term
- Multigrid for TV term
- Use Sherman-Morrison-Woodbury to combine

## Computation of terms like

- Adjoint problem is a final value problem
- Forward problem is an initial value problem
- But integral requires simultaneous availability of ** and *u*
- Options
- Store everything: *O(n)* storage and *O(n)* work
- Store nothing: *O(1)* storage, *O(n**2**)* work
- Direct sensitivities: *O(1)* storage, *O(n**4**)* work
- Metrics:
- Work unit = work per time step
- Storage unit = storage per time step
*n* = number of time  steps (assume explicit method)
- Griewank’s checkpointing algorithm
- Requires *O(log n)* storage, *O(n log n)* work

## Checkpointing Algorithm for Time-dependent Adjoints

*n/2*

*log n *levels* ** n/2 log n *work

## Behavior of Least Squares Objective

**Source frequency**

**Frequency of perturbation**
**of wave speed**

**Multiple minima**

**Ill-conditioning**

## A Multiscale Algorithm

- Determine material frequency components separately
- Solve inversion problem using source frequency and grid resolution that are consistent with material frequency
- Ill-conditioning: Regularize against material frequency changes greater than appropriate for current grid
- Multiple minima: Rely on good initial guesses from coarser grids to remain in basins of attraction of global minima
- Refine source and material frequency components and grid resolution