# Presentation: 197241

## Towards Optimal Complexity Algorithms for Optimal Design

- Barry Smith
- Mathematics and Computer Science Division
- Argonne National Laboratory
- Inspired by the work of George Biros and Omar
- Ghattas and presentations by David Keyes

## PDE-constrained Optimization

- Optimization time should be a small multiple (say 4) of the best PDE solution time
- Little software change between PDE simulation and optimization

- Goals

## Lagrangian Formulation

- State equation

- Adjoint equation

- Control equation

## Standard Methods: Nonlinear Iteration

- Given initial guess for w

-  Uses “block-box” PDE and adjoint PDE solvers 
-  Slowly convergent (block nonlinear Gauss-Seidel)
- Requires many PDE solves 

## Standard Methods: Newton Reduced SQP

- Standard Methods: Newton Reduced SQP

- Very fast (quadratic) convergence 
- Only requires “linearized” PDE solves,                    
- Forming Schur complement too expensive 

## Standard Methods: Quasi-Newton Reduced SQP

- Standard Methods: Quasi-Newton Reduced SQP

- Solve “approximate” Schur complement system

- Slow convergence, too many Newton steps 
- Only requires “linearized” PDE solves,                    
- Avoid computing Schur complement 

## Standard Methods: Full Space ApproachSAND:Simultaneous Analysis and Design“All at Once”Biros and Ghattas, Keyes and Gropp 87, Batterman and Heinkenschloss, many others

- Solve full nonlinear system via Newton using  preconditioned
-  Krylov method

- Few Newton steps 

## “All at Once”

- “All at Once”

- Precondition full linear system with an approximate Schur
-  complement (Quasi-Newton RSQP)

- Few Newton steps 
- Only requires “linearized” PDE solves,                    
- Performance depends on approximation to Schur complement 

## Proposed Approach

- Solve with grid sequencing, matrix-free Newton-
- Krylov(FGMRES)-(full)multigrid(matrix-free Krylov smoother)

- Few Newton steps 
- Only requires “linearized” PDE,               

## VariationPartial Matrix-based Smoothing 

- VariationPartial Matrix-based Smoothing

## Loop over grids coarsefine [Grid sequencing]
Construct data structures for grid 
(interpolate solution to grid from coarser grid)
Loop over Newton steps
Construct 
Loop over FGMRES cycles
Restrict fine grid residual to coarser grids
Loop over grids coarsefine [Full multigrid]
(interpolate solution to grid from coarser grid)
Loop over grids finecoarse (then back) [V-cycle]
     * Loop over GMRES smoothing steps
            + Loop over grid points applying
     * Form residual and restrict to coarser grid

- Construct data structures for grid
- (interpolate solution to grid from coarser grid)
- Loop over Newton steps
- Construct
- Loop over FGMRES cycles
- Restrict fine grid residual to coarser grids
- Loop over grids coarsefine [Full multigrid]
- (interpolate solution to grid from coarser grid)
- Loop over grids finecoarse (then back) [V-cycle]
-      * Loop over GMRES smoothing steps
-             + Loop over grid points applying
-      * Form residual and restrict to coarser grid

## Same Loop/Control Structure as Multigrid PDE Solver

- Field variables now include
-       are grid based, just like
-       are either “small”, grid based, or both
- If small, then store redundantly on all processors
- Interpolation/Restriction
-      same as       or as appropriate
- “small” kept on all grids
- otherwise use appropriate grid-based

## Abstract* Multigrid PDE Solver

- What minimal operations/objects does one need (on each level) to write an “abstract” multigrid PDE solver?
- Vectors for the field variables, x, b, r
- Linearized operator, A, to form residual r = b – A x
- Interpolation (restriction) operator
- Smoother (in our case requires only A)
- Function evaluation on field variables

- * Abstract = “general” = (easily) “reusable” = shockingly simple 

## Object Oriented Techniques I

- 1996 I wrote an “object oriented” class (called MG) to “contain” these operations/objects
- An array of these MGs (one MG for each level) is  all that was needed by the multigrid code
- The “methods/classes” in the MG object were “filled in” manually, for each grid level, by the application code.

## Object Oriented Techniques II

- Later I introduced a “grid” class that could produce
- A (empty) field vector
- An interpolation(restriction) between two grids
- A linearized operator
- A refined “grid”
- Given any “grid” object the MG object can be “filled in” , for each grid level automatically

## Using the Multigrid PDE Solver Code for the Optimization Problem

- Introduce “composite”
- Field vector that can contain the    ,    , and
- Interpolation(restriction)
- Operator on the composite field vector
- Composite “grid” that generates the above
- Packer/Unpacker that moves between the composite field vector and its components (for use by the function evaluation).

- Composite “grid” is just a “grid”
-  so can be used directly by the multigrid code 

## Sample Code*: PDE Solver App

- DACreate1d(...,&da)   // Create “grid” object
- MGCreate(...,&mg)  // Holder of MG object
- MGSetDM(mg,da)
- MGSetUseMatrixFree(mg)
- MGSetSNES(mg,PDEFunction,...)
- MGSolve(mg)

- *ugh

## Sample Code: PDEFunction

- ... // Determine my part of grid – xs, xm, ...
- DAGlobalToLocal (da,U, vu) // Get local ghosted part of U
- DAVecGetArray(da,vu,&u) //Access underlying field values
- ...
- for (i=xs; i<xs+xm; i++) {
- ...
-    fu[i]   = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h);

- }

- // Local
- function
- evaluation

## Sample Code: Optimization App

- DACreate1d(...,&da)   // Create “grid” object
- PackCreate(...,&pack)  // Create composite grid object
-    PackAddArray(pack,n) // Location to hold
-    PackAddDA(pack,da)                               //
-    PackAddDA(pack,da)                               //
- MGCreate(...,&mg)  // Holder of MG object
- MGSetDM(mg,pack)
- MGSetUseMatrixFree(mg)
- MGSetSNES(mg, Function,...)
- MGSolve(mg)

## Sample Code*:Optimization Function

- ... // Determine my part of grid – xs, xm, ...
- PackGetLocalVectors(pack,&w,&vu,&vl) // Get local
- PackScatter(pack,UWL,w,vu,vl)          // ghosted parts
- ...
- call LocalPDEFunction(...,w,vu,fu)   //
- call LocalPDEJacobian(...,w,vu,&J)  //
- call MatMultTranspose(J,L,vfl)
- for (i=xs; i<xs+xm; i++) {       //       and
-    fl [i]   += 2.*h*u[i]
-    fw[i]  = ...

- *ugh, again

## Numerical Results: Toy Problem

## Solution and Lagrange Multipliers

- X

## Floating Point Work

| Number of Levels | 3 | 5 | 7 | 9 | 11 | 13 | 15 |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Ratio of floating point work | 5.0 | 3.6 | 3.1 | 2.5 | 2.3 | 2.3 | 2.3 |

- Amount of work for optimization problem should be small multiple of work for PDE.

## Note on Symmetry

- Symmetric form has “pairs” of eigenvalues of opposite signs

- Symmetric form requires twice as many smoothing steps 
- Non-symmetric form requires use of   “non-symmetric” Krylov
-     method, e.g. GMRES 

## Conclusion

- Presented an approach to PDE constrained optimization that
- Offers the possibility of O(n) work 
- Reuses much of the software/algorithm infrastructure that already exists for solving PDEs with multigrid 
- Does not treat PDE solver as black-box, rather incorporates it directly into optimization algorithm, while retaining “separation of concerns”
- Needed. Much more numerical experience

## Thanks

- Fred Howes and the Department of Energy for funding our work on PETSc.

- www.mcs.anl.gov/petsc