# Iterative linear solvers for PDE-constrained optimization

Document Sample

```					         Iterative linear solvers for
PDE-constrained optimization involving
ﬂuid ﬂow

Andy Wathen
Oxford University

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

1 T
minx   2
x Ax   − xT b

equivalent to solving
Ax = b

Dave Grifﬁths Meeting, 2010 – p.2/32

1 T
minx   2
x Ax   − xT b

subject to Bx = c

Dave Grifﬁths Meeting, 2010 – p.3/32

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 Grifﬁths 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 Grifﬁths 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 ﬁnite elements/MAC ﬁnite difference approximation:

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

K : discrete Laplacian

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

Preconditioning?

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

Preconditioning?

Block diagonal/triangular preconditioners

Dave Grifﬁths 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 Grifﬁths 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 Grifﬁths Meeting, 2010 – p.7/32
S for Stokes problem:
Babuska-Brezzi stability ⇒ Schur Complement spectrally
equivalent to the ﬁnite element identity matrix, the mass
matrix, M ie. the Gram matrix of the ﬁnite element basis
functions {φj , j = 1, . . . , n} in L2 (Ω), Ω ⊂ Rd :

BA−1 B T ≈ ∇ · (∇2 )−1 ∇ ≈ Id → M
√
Speciﬁcally ∃γ > 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), Grifﬁths (’95))
Dave Grifﬁths 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) ﬁnite 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 :
nonlinear preconditioner
diagonally scaled Chebyshev (semi-)iteration: is linear
and we have precise eigenvalue inclusion intervals
(W & Rees (2008))                                   Dave Grifﬁths 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 Grifﬁths 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 ﬁnite 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 Grifﬁths 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 Grifﬁths 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 Grifﬁths 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 Grifﬁths 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: ﬁnite 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 Grifﬁths Meeting, 2010 – p.15/32
so, Lagrangian:
1    T        T    β
y My + y b +       uT M u + λT (Ky − M u − d)
2                  2
                       
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

A BT
B 0
Dave Grifﬁths 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 Grifﬁths 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 Grifﬁths 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 Grifﬁths 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 Grifﬁths 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 satisﬁed in this case
(Braess & Peisker (1986))

Dave Grifﬁths 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 Grifﬁths 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) ﬁnite 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 Grifﬁths 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 Grifﬁths 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 Grifﬁths 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 Grifﬁths 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 ﬁnite elements for (forward) Stokes problem:

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

Dave Grifﬁths 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 Grifﬁths 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 Grifﬁths 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 Grifﬁths Meeting, 2010 – p.30/32
Here: Q2-Q1 mixed ﬁnite elements for cavity ﬂow
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 Grifﬁths 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 indeﬁnite 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 Grifﬁths Meeting, 2010 – p.32/32

```
DOCUMENT INFO
Shared By:
Categories:
Stats:
 views: 13 posted: 4/4/2011 language: English pages: 35