92. Structure of Saddle-point Problems#
92.1. Constrained minimization problems#
Let \(A \in {\mathbf R}^{n\times n}\) be spd, and \(f \in {\mathbf R}^n\). Then the minimization problem
is equivalent to the linear system
Since \(A\) is spd, the function to be minimized is convex.
Now let \(B \in {\mathbf R}^{m\times n}\) of full rank, with \(m \leq n\) and, and \(g \in {\mathbf R}^m\). We consider the constrained minimization problem
The minimizers are found among the critical points of the Lagrange function
where \(p \in {\mathbb R}^m\) is the Lagrange parameter. Differentiation with respect to \(u\) and \(p\), and zeroing leads to the Karush-Kuhn-Tucker (KKT) equations
The solution is not a minimizer, but a saddle-point of the Lagrange function \(L(.,.)\).
Theorem: The Hessian matrix of \(L(.,.)\),
\[\begin{split} \left( \begin{array}{cc} A & B \\ B & \end{array} \right) \in {\mathbf R}^{(m+n) \times (m+n)} \end{split}\]has \(n\) positive and \(m\) negative eigenvalues.
Proof: We define the sub-spaces
and evaluate for \(u \neq 0\)
and for \(p \neq 0\)
The matrix \(S := B A^{-1} B^T\) is called the (negative) Schur complement. For an spd matrix \(A\) and for the rank-\(m\) matrix \(B\) it is positive definit.
Since we found an \(n\)-dimensional sub-space where the quadratic form is positive, and an \(m\)-dimensional sub-space where it is negative, the \(\min-\max\)-theorem proves that we have \(n\) positive and \(m\) negative eigenvalues.
92.2. Examples#
Stokes Equation
The vector-field \(u \in [H^1_{0,D}]^2\) is the velocity of an incompressible flow. It minimizes the \(H^1\)-norm, under the incompressibility constraint \(\operatorname{div} u = 0\):
The variational formulation of the KKT system is
Discretization leads to the matrix equation
where \(A\) is the discretization of the \(-\Delta\)-operator for both components of the velocity, and \(B\) is the discretization of the \(\operatorname{div}\)-operator. It’s transposed matrix \(B^T\) is the discretization of the negative \(\nabla\)-operator. The physical meaning of the Lagrange parameters \(p\) is pressure.
Dirichlet boundary conditions
We can interpret a Dirichlet problem as a constrained minimization problem
The KKT system is to find \(u \in H^1(\Omega)\) and \(\lambda\) on the boundary such that
By integration by parts in the first equation, we obtain
We see that the meaning of the Lagrange parameter \(\lambda\) is the negative normal derivative \(- \tfrac{\partial u}{\partial n}\)
92.3. Schur complement iteration#
We solve the block-system
Since \(A\) is spd and thus regular, we can solve for \(u\) in the first block-row
Inserting this \(u\) into the second block-row gives
or equivalently
The Schur complement matrix \(S = B A^{-1} B^T\) is positive definite. The Uzawa algorithm (and there are many variants of it) is now to use a Richardson iteration, or conjugate gradients, to solve the Schur complement equation
with \(\tilde g = B A^{-1} f - g\). After computing the Lagrange parameter \(p\), the first component \(u\) is found from \(u = A^{-1} (f - B^T p)\).
Of course, the drawback of this method is that it requires the inverse \(A^{-1}\), which is usually not cheaply available. If the Uzawa iteration does not take many iteration steps, one can also use an (inner) iterative solver to apply \(A^{-1}\) in every step.
To design preconditioners for the Uzawa iteration, it is important to understand the norm induced by the Schur-complement \(S\):
If \(A\) and \(B\) come from finite-element discretization of bilinear-forms \(a(.,.)\) and \(b(.,.)\), then the last term can be recast as
If we can prove continuity of \(b(.,.)\)
and the famous LBB-condition
we immediately get
If we have a preconditioner for the norm-matrix of the \(Q\)-norm, we have a preconditioner for the Schur complement \(S\).
from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fesu = H1(mesh, order=1)
feslam = H1(mesh, order=1, definedon=mesh.Boundaries(".*"))
u,v = fesu.TnT()
lam,mu = feslam.TnT()
a = BilinearForm(grad(u)*grad(v)*dx + u*v*dx).Assemble()
b = BilinearForm(u*mu*ds).Assemble()
f = LinearForm(10*v*dx).Assemble()
g = LinearForm((x+y)*mu*ds).Assemble()
ainv = a.mat.Inverse()
gtilde = (b.mat @ ainv * f.vec - g.vec).Evaluate()
S = b.mat @ ainv @ b.mat.T
preS = Projector(feslam.FreeDofs(), True)
from ngsolve.krylovspace import CGSolver
Sinv = CGSolver(S, preS, printrates=True, maxiter=200)
gflam = GridFunction(feslam)
gfu = GridFunction(fesu)
gflam.vec.data = Sinv * gtilde
gfu.vec.data = ainv * (f.vec - b.mat.T * gflam.vec)
Draw (gfu);
CG iteration 1, residual = 5.70374341015521
CG iteration 2, residual = 0.3225362849898911
CG iteration 3, residual = 0.08339010831277839
CG iteration 4, residual = 0.01117390734935611
CG iteration 5, residual = 0.0047927265813528675
CG iteration 6, residual = 0.004039693809759069
CG iteration 7, residual = 0.003477548920994166
CG iteration 8, residual = 0.004676900022939997
CG iteration 9, residual = 0.0028443097964922547
CG iteration 10, residual = 0.0010427358439432045
CG iteration 11, residual = 0.0015208669764192647
CG iteration 12, residual = 0.0010007715432088899
CG iteration 13, residual = 0.0008083236800528286
CG iteration 14, residual = 0.00038646050624300133
CG iteration 15, residual = 0.00023345852649526882
CG iteration 16, residual = 0.00019121655931019407
CG iteration 17, residual = 0.0002507057465019636
CG iteration 18, residual = 7.906795223588332e-05
CG iteration 19, residual = 6.20599823174278e-05
CG iteration 20, residual = 5.296784849054979e-05
CG iteration 21, residual = 3.7249857424010985e-05
CG iteration 22, residual = 2.233427661186067e-05
CG iteration 23, residual = 2.5759653496097197e-05
CG iteration 24, residual = 9.522437201035332e-06
CG iteration 25, residual = 7.208353059797299e-06
CG iteration 26, residual = 3.650418161152721e-06
CG iteration 27, residual = 6.150077558234375e-06
CG iteration 28, residual = 4.63861041586411e-06
CG iteration 29, residual = 9.558709031922435e-07
CG iteration 30, residual = 7.351143807631565e-07
CG iteration 31, residual = 4.395662608619207e-06
CG iteration 32, residual = 1.728559224122473e-07
CG iteration 33, residual = 1.9409072280009905e-07
CG iteration 34, residual = 5.965516822381169e-08
CG iteration 35, residual = 3.763709263634211e-08
CG iteration 36, residual = 1.98928958686514e-08
CG iteration 37, residual = 2.6624750020974526e-08
CG iteration 38, residual = 7.27342755478678e-08
CG iteration 39, residual = 3.348927688529326e-09
CG iteration 40, residual = 2.943042808256458e-09
CG iteration 41, residual = 1.6193400996968746e-09
CG iteration 42, residual = 1.9205838369725923e-09
CG iteration 43, residual = 3.6576873786252613e-09
CG iteration 44, residual = 5.15316554116041e-10
CG iteration 45, residual = 4.938952741759179e-10
CG iteration 46, residual = 8.122531987166569e-10
CG iteration 47, residual = 3.306527104434466e-10
CG iteration 48, residual = 2.4383600045391037e-10
CG iteration 49, residual = 1.9216564883158577e-11
CG iteration 50, residual = 4.198500568785393e-11
CG iteration 51, residual = 8.843758072882099e-12
CG iteration 52, residual = 1.1584849958238028e-11
CG iteration 53, residual = 9.66978184735845e-13
If we experiment with the mesh-size, we observe that the iteration number doubles, if we refine the mesh by a factor of 4, which numerically indicate that
We will come back to this, later.
92.4. Block-preconditioning#
Now we want to solve directly the saddle-point system, without forming the Schur complement explicitly. We want to use preconditioners for both fields separately.
First, we show that the spectrum of
consists just of a few, small eigenvalues, and thus this matrix is well conditioned.
Theorem The spectrum \(\sigma\) of the generalized eigenvalue problem
\[\begin{split} \left( \begin{array} {cc} A & B^T \\ B & 0 \end{array} \right) \left( \begin{array}{cc} u \\ p \end{array} \right) = \lambda \left( \begin{array} {cc} A & 0 \\ 0 & S \end{array} \right) \left( \begin{array}{cc} u \\ p \end{array} \right) \end{split}\]satisfies \(\sigma \subset \{ 1, (1+\sqrt{5})/2, (1-\sqrt{5})/2 \}\).
Proof: The eigenvalue problem reads as
with \((u,p) \neq 0\). The first row reads as
or $\( B^T p = (\lambda-1) A u. \)\( This implies \)\( B A^{-1} B^T p = (\lambda-1) B u \)\( or \)\( S p = (\lambda-1) B u \)$
Inserting the second equation \(Bu = \lambda S p\), we obtain
If \(p\) is non-zero, the eigenvalue \(\lambda\) must satisfy the quadratic equation
i.e. \(\lambda = \frac{1}{2} (1 \pm \sqrt{5})\). If \(p = 0\), but \(u \neq 0\), the equation \(B^T p = (\lambda-1) A u\) can be satisfied only if \(\lambda = 1\).
Since
is symmetric, and
provides a norm, the fraction $\( \frac{\max | \sigma |}{\min | \sigma |} = \frac{1+\sqrt{5}}{| 1-\sqrt{5} |} = 2.618... \)$ bounds the norm of the preconditioned matrix. It also bounds the number of iterations of a Krylov-space method. However, since eigenvalues are positive and negative, we cannot use conjugate gradients. One can apply MinRes, or GMRes.
from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fesu = H1(mesh, order=1, autoupdate=True)
feslam = H1(mesh, order=1, definedon=mesh.Boundaries(".*"), autoupdate=True)
u,v = fesu.TnT()
lam,mu = feslam.TnT()
a = BilinearForm(grad(u)*grad(v)*dx + u*v*dx)
b = BilinearForm(u*mu*ds)
pre = Preconditioner(a, type="multigrid")
a.Assemble()
b.Assemble()
for i in range(2):
mesh.Refine()
a.Assemble()
b.Assemble()
f = LinearForm(10*v*dx).Assemble()
g = LinearForm((x+y)*mu*ds).Assemble()
K = BlockMatrix( [[a.mat, b.mat.T], [b.mat, None]] )
F = BlockVector( [f.vec, g.vec])
U = K.CreateRowVector()
prelam = Projector(feslam.FreeDofs(), True)
PRE = BlockMatrix( [[pre.mat, None], [None, prelam]] )
from ngsolve.krylovspace import MinRes
U.data = MinRes(K, pre=PRE, rhs=F, maxsteps=400)
LinearSolver iteration 1, residual = 10.006045236381954
LinearSolver iteration 2, residual = 2.6923699593826087
LinearSolver iteration 3, residual = 1.5965294423392078
LinearSolver iteration 4, residual = 0.17922278593711807
LinearSolver iteration 5, residual = 0.16553184661052575
LinearSolver iteration 6, residual = 0.12886928954088694
LinearSolver iteration 7, residual = 0.041929015261395405
LinearSolver iteration 8, residual = 0.040923325335253814
LinearSolver iteration 9, residual = 0.038049450876329084
LinearSolver iteration 10, residual = 0.0076180294093657385
LinearSolver iteration 11, residual = 0.005474471005680139
LinearSolver iteration 12, residual = 0.005471921389164226
LinearSolver iteration 13, residual = 0.004824883063804365
LinearSolver iteration 14, residual = 0.0019283965718213832
LinearSolver iteration 15, residual = 0.001817186787755464
LinearSolver iteration 16, residual = 0.0018122634640929303
LinearSolver iteration 17, residual = 0.0013064909694559504
LinearSolver iteration 18, residual = 0.0011031206144915993
LinearSolver iteration 19, residual = 0.0011026808151522044
LinearSolver iteration 20, residual = 0.0008766395383082775
LinearSolver iteration 21, residual = 0.0007898291411921694
LinearSolver iteration 22, residual = 0.0007895774460415798
LinearSolver iteration 23, residual = 0.0006500874307245643
LinearSolver iteration 24, residual = 0.0005612430740709964
LinearSolver iteration 25, residual = 0.0005610940467488651
LinearSolver iteration 26, residual = 0.0004940130723187382
LinearSolver iteration 27, residual = 0.00047398999148945165
LinearSolver iteration 28, residual = 0.0004724790675118711
LinearSolver iteration 29, residual = 0.0003340640559150125
LinearSolver iteration 30, residual = 0.0002954863505778272
LinearSolver iteration 31, residual = 0.00029407389742718666
LinearSolver iteration 32, residual = 0.0002858951744502473
LinearSolver iteration 33, residual = 0.0002858943477020491
LinearSolver iteration 34, residual = 0.0002573396840026226
LinearSolver iteration 35, residual = 0.00016632874467311973
LinearSolver iteration 36, residual = 0.00016535344657900392
LinearSolver iteration 37, residual = 0.00016497728585449368
LinearSolver iteration 38, residual = 0.0001325172888976091
LinearSolver iteration 39, residual = 0.00010587170753164257
LinearSolver iteration 40, residual = 0.00010540546822764906
LinearSolver iteration 41, residual = 0.00010523225336641586
LinearSolver iteration 42, residual = 8.45124800338088e-05
LinearSolver iteration 43, residual = 7.218280631649127e-05
LinearSolver iteration 44, residual = 7.212874993302292e-05
LinearSolver iteration 45, residual = 7.164508973011206e-05
LinearSolver iteration 46, residual = 5.951829471843949e-05
LinearSolver iteration 47, residual = 5.812467620620077e-05
LinearSolver iteration 48, residual = 5.809920489294807e-05
LinearSolver iteration 49, residual = 5.256468005904873e-05
LinearSolver iteration 50, residual = 4.827573573321181e-05
LinearSolver iteration 51, residual = 4.826355173665311e-05
LinearSolver iteration 52, residual = 4.738513471387175e-05
LinearSolver iteration 53, residual = 4.3641327621771723e-05
LinearSolver iteration 54, residual = 4.359705537916172e-05
LinearSolver iteration 55, residual = 4.329599026830253e-05
LinearSolver iteration 56, residual = 3.522289251977826e-05
LinearSolver iteration 57, residual = 3.456642835152211e-05
LinearSolver iteration 58, residual = 3.443060784566924e-05
LinearSolver iteration 59, residual = 3.424224823412509e-05
LinearSolver iteration 60, residual = 3.4221850132279816e-05
LinearSolver iteration 61, residual = 3.1073627145448555e-05
LinearSolver iteration 62, residual = 2.9707856164109158e-05
LinearSolver iteration 63, residual = 2.970781611564814e-05
LinearSolver iteration 64, residual = 2.8655448637711687e-05
LinearSolver iteration 65, residual = 2.6921414327920445e-05
LinearSolver iteration 66, residual = 2.691802293153496e-05
LinearSolver iteration 67, residual = 2.6154696719059098e-05
LinearSolver iteration 68, residual = 2.4617985894611936e-05
LinearSolver iteration 69, residual = 2.461457965092768e-05
LinearSolver iteration 70, residual = 2.253928498791778e-05
LinearSolver iteration 71, residual = 2.2111034696236826e-05
LinearSolver iteration 72, residual = 2.2110850172488514e-05
LinearSolver iteration 73, residual = 2.1271994859472737e-05
LinearSolver iteration 74, residual = 1.798256167052536e-05
LinearSolver iteration 75, residual = 1.7930075586512082e-05
LinearSolver iteration 76, residual = 1.792353524357799e-05
LinearSolver iteration 77, residual = 1.7176708938530337e-05
LinearSolver iteration 78, residual = 1.6596824476560285e-05
LinearSolver iteration 79, residual = 1.6595556579047295e-05
LinearSolver iteration 80, residual = 1.6355511573647056e-05
LinearSolver iteration 81, residual = 1.4540774217018864e-05
LinearSolver iteration 82, residual = 1.4511697025156621e-05
LinearSolver iteration 83, residual = 1.3951502118854358e-05
LinearSolver iteration 84, residual = 1.326173089146925e-05
LinearSolver iteration 85, residual = 1.3260842185718308e-05
LinearSolver iteration 86, residual = 1.3205006315114814e-05
LinearSolver iteration 87, residual = 1.3203779858822644e-05
LinearSolver iteration 88, residual = 1.3043058345204111e-05
LinearSolver iteration 89, residual = 1.2981188844192642e-05
LinearSolver iteration 90, residual = 1.2967281626708148e-05
LinearSolver iteration 91, residual = 1.1069502043680783e-05
LinearSolver iteration 92, residual = 1.0462092790837482e-05
LinearSolver iteration 93, residual = 1.0461465581338192e-05
LinearSolver iteration 94, residual = 1.0288591659652555e-05
LinearSolver iteration 95, residual = 8.926724817000584e-06
LinearSolver iteration 96, residual = 8.866807855706634e-06
LinearSolver iteration 97, residual = 8.865276870745839e-06
LinearSolver iteration 98, residual = 8.363092579360034e-06
LinearSolver iteration 99, residual = 7.3464697438074095e-06
LinearSolver iteration 100, residual = 7.340704948792019e-06
LinearSolver iteration 101, residual = 7.29521483369701e-06
LinearSolver iteration 102, residual = 6.184210630157364e-06
LinearSolver iteration 103, residual = 6.124755278070719e-06
LinearSolver iteration 104, residual = 6.123961802106049e-06
LinearSolver iteration 105, residual = 5.703302836149492e-06
LinearSolver iteration 106, residual = 5.04018831668199e-06
LinearSolver iteration 107, residual = 5.037025305243992e-06
LinearSolver iteration 108, residual = 5.01070297777827e-06
LinearSolver iteration 109, residual = 4.929891724395826e-06
LinearSolver iteration 110, residual = 4.862010175644464e-06
LinearSolver iteration 111, residual = 4.047449269515496e-06
LinearSolver iteration 112, residual = 4.000959130828542e-06
LinearSolver iteration 113, residual = 4.000948242141481e-06
LinearSolver iteration 114, residual = 3.914586985583179e-06
LinearSolver iteration 115, residual = 3.2589358798309174e-06
LinearSolver iteration 116, residual = 3.2471057090298423e-06
LinearSolver iteration 117, residual = 3.221296164560807e-06
LinearSolver iteration 118, residual = 2.651391827660408e-06
LinearSolver iteration 119, residual = 2.6033751449924753e-06
LinearSolver iteration 120, residual = 2.6033737951279255e-06
LinearSolver iteration 121, residual = 2.5674527138563625e-06
LinearSolver iteration 122, residual = 2.532953755240753e-06
LinearSolver iteration 123, residual = 2.5300820511137353e-06
LinearSolver iteration 124, residual = 2.0500025237608006e-06
LinearSolver iteration 125, residual = 1.8809856817251082e-06
LinearSolver iteration 126, residual = 1.8809291038613847e-06
LinearSolver iteration 127, residual = 1.829374299982996e-06
LinearSolver iteration 128, residual = 1.7423315064199907e-06
LinearSolver iteration 129, residual = 1.7422606950747696e-06
LinearSolver iteration 130, residual = 1.6827006714301568e-06
LinearSolver iteration 131, residual = 1.5650250702445693e-06
LinearSolver iteration 132, residual = 1.5644977355604277e-06
LinearSolver iteration 133, residual = 1.5505690933947208e-06
LinearSolver iteration 134, residual = 1.5071629470089156e-06
LinearSolver iteration 135, residual = 1.4995223152704887e-06
LinearSolver iteration 136, residual = 1.4315558039131422e-06
LinearSolver iteration 137, residual = 1.4309202702438207e-06
LinearSolver iteration 138, residual = 1.3640081970409472e-06
LinearSolver iteration 139, residual = 1.2545477664918666e-06
LinearSolver iteration 140, residual = 1.2542239242445278e-06
LinearSolver iteration 141, residual = 1.2357693872083088e-06
LinearSolver iteration 142, residual = 1.1704864814003705e-06
LinearSolver iteration 143, residual = 1.1704603209357195e-06
LinearSolver iteration 144, residual = 1.0684054021517988e-06
LinearSolver iteration 145, residual = 9.96581894237824e-07
gfu = GridFunction(fesu)
gfu.vec.data = U[0]
Draw (gfu);