Iterative linear solvers for PDE-constrained optimization

Document Sample
Iterative linear solvers for PDE-constrained optimization Powered By Docstoc
					         Iterative linear solvers for
    PDE-constrained optimization involving
                  fluid flow

                     Andy Wathen
                    Oxford University




joint work with
Tyrone Rees, Martin Stoll and Sue Thorne
                                           Dave Griffiths Meeting, 2010 – p.1/32
Let A ∈ Rn×n be symmetric and positive definite

                           1 T
                    minx   2
                             x Ax   − xT b

equivalent to solving
                           Ax = b




                                                 Dave Griffiths Meeting, 2010 – p.2/32
Constrained Optimization and Saddle-point problems:

                         1 T
                  minx   2
                           x Ax   − xT b

                   subject to Bx = c




                                               Dave Griffiths Meeting, 2010 – p.3/32
Constrained Optimization and Saddle-point problems:

                                 1 T
                      minx       2
                                   x Ax   − xT b

                        subject to Bx = c

Lagrangian: L(x, λ) = 1 xTAx − xT b + (Bx − c)T λ
                      2

λ : Lagrange multipliers.

 minx     ⇒          Ax + B ⋆ λ      =b
 maxλ     ⇒          Bx              =c


                       A BT           x            b
                 ⇒                          =
                       B 0            λ            c
(Benzi, Golub & Liesen (2005))
                                                       Dave Griffiths Meeting, 2010 – p.3/32
Classic problem of this type: the Stokes problem
‘Minimise energy subject to conserving mass’ ⇒


                          −∇2 y + ∇p = u
                 subject to     ∇·y = 0

y: velocity,   p: pressure,   u: body forces




                                                   Dave Griffiths Meeting, 2010 – p.4/32
Classic problem of this type: the Stokes problem
‘Minimise energy subject to conserving mass’ ⇒


                          −∇2 y + ∇p = u
                 subject to     ∇·y = 0

y: velocity,   p: pressure,       u: body forces

Mixed finite elements/MAC finite difference approximation:

     K BT         y           u               K 0
                       =           ,   K=              in R2
     B 0          p           g               0 K

K : discrete Laplacian


                                                    Dave Griffiths Meeting, 2010 – p.4/32
                       K BT
                       B 0
is symmetric indefinite → MINRES Krylov subspace
iterative solver (Paige & Saunders (1975))

Preconditioning?




                                             Dave Griffiths Meeting, 2010 – p.5/32
                          K BT
                          B 0
is symmetric indefinite → MINRES Krylov subspace
iterative solver (Paige & Saunders (1975))

Preconditioning?

Block diagonal/triangular preconditioners




                                             Dave Griffiths Meeting, 2010 – p.5/32
Block diagonal/triangular preconditioners:
based on observation (Murphy, Golub & W (2000), Korzak(1999))

                           A BT
                           B 0

preconditioned by
                                                    √
     A 0
 •           has 3 distinct eigenvalues (1, 1 ±
                                            2
                                                       5
                                                      2 )
     0 S
      A BT
 •              has 2 distinct eigenvalues
      0 S

where S = BA−1 B T (Schur Complement)
⇒ MINRES /GMRES terminates in 3 / 2 iterations
⇒ want approximations A, S ⇒ 3 / 2 clusters
⇒ fast convergence                                   Dave Griffiths Meeting, 2010 – p.6/32
For Stokes problem:
A = K is discrete Laplacians: use A =multigrid cycles

 • geometric mutigrid: relaxed Jacobi smoothing, standard
   grid transfers
 • algebraic multigrid: HSL routine HSL MI20
   (Boyle, Mihajlovic & Scott (2007))




                                                Dave Griffiths Meeting, 2010 – p.7/32
S for Stokes problem:
Babuska-Brezzi stability ⇒ Schur Complement spectrally
equivalent to the finite element identity matrix, the mass
matrix, M ie. the Gram matrix of the finite element basis
functions {φj , j = 1, . . . , n} in L2 (Ω), Ω ⊂ Rd :

          BA−1 B T ≈ ∇ · (∇2 )−1 ∇ ≈ Id → M
                       √
Specifically ∃γ > 0, Γ ≤ d such that

                            xT BA−1 B T x
                      γ≤                       ≤Γ
                                 xT M x

          where         M = {mi,j },        mi,j =          φi φj
                                                        Ω


(Bank, Welfert & Yserentant (’90), Silvester & W (’94), Griffiths (’95))
                                                                Dave Griffiths Meeting, 2010 – p.8/32
Mass matrix is effectively preconditioned by its diagonal
D =diag(M )     (W (1987))
eg. for Q1 (bilinear) elements in 2-dimensions (rectangles)
eigenvalues of D −1 M all in [1/4, 9/4]
eg. for Q1 (trilinear) finite elements in 3-dimensions (bricks)
eigenvalues of D −1 M all in [1/8, 27/8]
independently of mesh size h (ie. independently of discrete
problem dimension)
Could use S = D but better a few iterations of a diagonally
preconditioned iteration for M :
    diagonally scaled Conjugate Gradients: leads to a
    nonlinear preconditioner
    diagonally scaled Chebyshev (semi-)iteration: is linear
    and we have precise eigenvalue inclusion intervals
(W & Rees (2008))                                   Dave Griffiths Meeting, 2010 – p.9/32
           4
          10
                                               CG
                                               Chebyshev
           2
          10
  −1
    A




           0
|| r ||




          10
    k




           −2
          10



           −4
          10
               0   5            10        15             20
                       Number of iterations
                                                       Dave Griffiths Meeting, 2010 – p.10/32
For Stokes problem:
      A 0             AAM G       0
P =             =
      0 S               0   T20 (D −1 M )
Here: Q2-Q1 mixed finite element (Taylor-Hood),
      A = AAM G : 1 AMG V-cycle

         h               Iterations   CPU time (s)
         2−2   (187)         19         0.015
         2−3   (659)         24         0.073
         2−4   (2,467)       26         0.082
         2−5   (9,539)       28          0.21
         2−6   (37,507)      29          3.80
         2−7   (148,739)     29          15.5

                                                     Dave Griffiths Meeting, 2010 – p.11/32
             PDE-constrained Optimization
General problem:
Given Ω ⊂ R2 or R3 , y ∈ L2 (Ω) as some desired state
then for some (regularisation) parameter β

                     1         2            β       2
              min        y−y   L2 (Ω)   +       u   L2 (Ω)
               y,u   2                      2
subject to

               Ly = u     in Ω,   y=y           on ∂Ω


where L represents a partial differential operator



                                                             Dave Griffiths Meeting, 2010 – p.12/32
can also include simple bounds on the control:

                         u≤u≤u

and even bounds on the state (via Moreau-Yoshida

regulrization)
                         y≤y≤y



also boundary control. . .




                                                 Dave Griffiths Meeting, 2010 – p.13/32
Simple sample problem:
desirable y, controlable body force u

                    1         2            β       2
              min       y−y   L2 (Ω)   +       u   L2 (Ω)
              y,u   2                      2
subject to

             −∇2 y = u    in Ω,        y=y         on ∂Ω




                                                            Dave Griffiths Meeting, 2010 – p.14/32
                           1             2       β        2
                    min            y−y       +       u
                    y,u    2                     2
   subject to         − ∇2 y = u             in Ω,        y=y     on ∂Ω
Discretisation: finite elements
yh =    yj φj , y = (y1 , y2 , . . . , yn )T
uh =    uj φj , u = (u1 , u2 , . . . , un )T

                      1    T             T           β
             min          y My + y b +                   uT M u
              y,u     2                              2
              subject to               Ky = M u + d


M = {mi,j }, mi,j =            Ω   φi φj — mass matrix
K = {ki,j }, ki,j =       Ω   ∇φi · ∇φj — stiffness matrix
as before                                                         Dave Griffiths Meeting, 2010 – p.15/32
so, Lagrangian:
    1    T        T    β
        y My + y b +       uT M u + λT (Ky − M u − d)
    2                  2
stationarity ⇒ Saddle point system
                                   
              M   0     K T      y   −b
             0 βM −M   u  =  0 
              K −M       0       λ   d


                                            M  0
   Note      B=   K −M          and   A=
                                            0 βM

in usual saddle point form

                             A BT
                             B 0
                                                Dave Griffiths Meeting, 2010 – p.16/32
                                             M  0
  Recall      B=     K −M         and   A=
                                             0 βM

so S ?       S = BA−1 B T (Schur Complement)


                            M −1      0      KT
         =        K −M             1
                             0     β M −1    −M
              1
         =        M + KM −1 K T
              β

Unless approx β < 10−6 dominant part is S = KM −1 K T




                                               Dave Griffiths Meeting, 2010 – p.17/32
Hence preconditioner for
                                       
      M     0     K T       M 0    0
A =  0 βM −M  is P =  0 βM      0      
       K −M        0        0 0 KM −1 K T

Eigenvalues ν of P −1 A

                                 ν = 1,
                            
1                  2α1 h4                  1                2α2
     1 +     5+             ≤ ν ≤             1+    5+
2                    β                     2                   β
                                                                               
     1               2α2                   1                2α1 h4
or       1−    5+                ≤ ν   ≤       1 −   5+                        ,
     2                   β                 2                       β

where α1 , α2 are positive constants independent of h.
                                                           Dave Griffiths Meeting, 2010 – p.18/32
                          
      M  0    0
P =  0 βM    0      ,         β = 10−2
      0  0 KM −1 K T



   1.5




    1




   0.5




    0




  −0.5


         0   5   10   15   20   25
                                     Dave Griffiths Meeting, 2010 – p.19/32
But                                        
                   M  0    0
             P =  0 βM    0      
                   0  0 KM −1 K T
still expensive to use in practice so employ approximations
                 M ≃M       and K ≃ K
giving
                                     
               M  0    0
                                           A 0
         P =  0 βM    0      =
                                             0 S
               0  0 KM −1 K T



                                                  Dave Griffiths Meeting, 2010 – p.20/32
Important subtlety:
                              K≃K
does not imply that

                   KM −1 K T ≃ KM −1 K T

or indeed that
                            K K T ≃ KK T
without further conditions which are satisfied in this case
(Braess & Peisker (1986))




                                                   Dave Griffiths Meeting, 2010 – p.21/32
                                    
             M  0    0
                                          A 0
       P =  0 βM    0      =
                                            0 S
             0  0 KM −1 K T


so M ?: T20 (D −1 M ) as before

For L = −∇2 , K is a discrete Laplacian: use multigrid
cycles as before

In our examples:

K is the action of 2 V-cycles

M is the action of 20 Chebyshev semi-iterative steps



                                                 Dave Griffiths Meeting, 2010 – p.22/32
Example problem: Ω = [0, 1]d , d = 2, 3

                    1          2                 β       2
              min       y−y    L2 (Ω)      +         u   L2 (Ω)
              y,u   2                            2
subject to

             −∇2 y = u       in Ω,         y=y           on ∂Ω

Q1 (bilinear) finite elements, β = 10−2

                                      1


                        y:           0.8


                                     0.6


                                     0.4


                                     0.2


                                      0
                                      1

                                                                                       1
                                               0.5                              0.8
                                                                         0.6
                                                                   0.4
                                                             0.2
                                                     0   0


                                                                         Dave Griffiths Meeting, 2010 – p.23/32
CPU times (MINRES iterations) in 2D, tol 10−6
  h       3n backslash    MINRES MINRES
                          (PAM G )   (PM G )
2−2       27     0.0003    0.02 (7) 0.13 (7)
2−3      147      0.002    0.03 (9) 0.16 (9)
2−4      675       0.01    0.05 (9) 0.21 (9)
2−5     2883       0.08    0.14 (9) 0.41 (9)
2−6    11907       0.46    0.61 (9) 1.29 (9)
2−7    48387       3.10    2.61 (9) 5.09 (9)
2−8   195075       15.5   15.0 (11) 23.6 (9)
2−9   783363         —    75.6 (11) 136 (9)




                                         Dave Griffiths Meeting, 2010 – p.24/32
Control:


              0


           −0.05


            −0.1


           −0.15


            −0.2
               1
                   0.8                                                     1
                         0.6                                         0.8
                               0.4                             0.6
                                                         0.4
                                     0.2           0.2
                                           0   0




                                                                               Dave Griffiths Meeting, 2010 – p.25/32
CPU times (MINRES iterations) in 3D, tol 10−6
  h       3n backslash MINRES MINRES
                       (PAM G )  (PM G )
2−2       81     0.001 0.02 (7) 0.14 (8)
2−3     1029     0.013 0.13 (9) 0.26 (8)
2−4    10125      25.1 1.89 (8) 1.69 (8)
2−5    89373        — 22.1 (8) 15.9 (8)
2−6   750141        —        — 230 (10)




                                         Dave Griffiths Meeting, 2010 – p.26/32
                          Stokes Control
            1         2            1         2            β         2
    min           ˆ
                y−y   L2 (Ω)   +         ˆ
                                       p−p   L2 (Ω)   +       u     L2 (Ω)
    y,p,u   2                      2                      2


                subject to         − ∇2 y + ∇p = u
                                          ∇·y = 0

y: velocity,    p: pressure.
Mixed finite elements for (forward) Stokes problem:

      K BT            y            u             K 0
                           =           , K=                         in R2
      B 0             p            g             0 K




                                                                  Dave Griffiths Meeting, 2010 – p.27/32
Cost functional
     1    T       T       1    T       T       β
         y My y − y b +       p Mp p − p d +       u T Mu u
     2                    2                    2
combined with constraint via the Lagrangian ⇒
                                                          
     My 0         0       K     BT      y               b
    0 Mp         0       B       0  p              d      
                                                          
    0      0 βMu −Mu 0   u  =                      0      .
                                                          
    K B     T −M
                    u      0      0  λ                 h
      B     0     0        0      0     µ               k




                                                      Dave Griffiths Meeting, 2010 – p.28/32
Schur Complement:

       K BT        My−1  0          K BT             1
                          −1                   +         Mu
       B 0          0   Mp          B 0            β

and again ignore 2nd term for moderate β

So overall block diagonal preconditioner requires:
  My , Mp , Mu → Chebyshev
  and 2 Stokes approximations




                                                     Dave Griffiths Meeting, 2010 – p.29/32
Stokes preconditioners:

                   K BT            K 0
                              =
                   B 0             B Mp

and
                   K BT            K BT
                              =
                   B 0             0 Mp

on left and right respectively where K is multigrid cycles for
each discrete scalar Laplacian as before
(Silvester & W (1993), Klawonn (1998))
Gives symmetric Schur complement approximation


                                                   Dave Griffiths Meeting, 2010 – p.30/32
Here: Q2-Q1 mixed finite elements for cavity flow
      4 AMG V-cycles to approx each K
      20 Chebyshev semi-iterations for each M

         h               Iterations   CPU time (s)
         2−2   (344)         26          0.48
         2−3   (1,512)       31          1.05
         2−4   (6,344)       33          3.69
         2−5   (25,992)      33          18.0
         2−6   (105,224)     34          84.2
         2−7   (423,432)     34          342


recall: 29 MINRES iterations for forward Stokes solve (but
only 1 AMG V-cycle) ⇒ approx 10 times more work to solve
the control problem than a single PDE solve.
                                                     Dave Griffiths Meeting, 2010 – p.31/32
                             References
Rees, T., Dollar, H.S. & Wathen, A.J., 2009, ‘Optimal solvers for
PDE-constrained Optimization’, to appear in SIAM J. Sci. Comput.,
Oxford University NA report 08/10.

Murphy, M.F., Golub, G.H. & Wathen, A.J., 2000, ‘A note on
preconditioning for indefinite linear systems’, SIAM J. Sci. Comput.
21(6), pp. 1969-1972.

Silvester, D.J. & Wathen, A.J., 1994, ‘Fast iterative solution of stabilised
Stokes systems Part II: Using general block preconditioners’, SIAM J.
Numer Anal. 31, 1352-1367.

Wathen, A.J. & Rees, T., 2009, ‘Chebyshev semi-iteration in
preconditioning’, Victor Pereyra special issue of ETNA.

Benzi, M, Golub, G.H. & Liesen, J., 2005, ‘Numerical Solution of Saddle
Point Problems’, Acta Numerica 2005.

                                                                 Dave Griffiths Meeting, 2010 – p.32/32

				
DOCUMENT INFO