We consider a saddle point problem of the form

${\displaystyle {\begin{pmatrix}A&B\\B^{*}&\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\end{pmatrix}}={\begin{pmatrix}b_{1}\\b_{2}\end{pmatrix}},}$ 

where ${\displaystyle A}$  is a symmetric positive-definite matrix. Multiplying the first row by ${\displaystyle B^{*}A^{-1}}$  and subtracting from the second row yields the upper-triangular system

${\displaystyle {\begin{pmatrix}A&B\\&-S\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\end{pmatrix}}={\begin{pmatrix}b_{1}\\b_{2}-B^{*}A^{-1}b_{1}\end{pmatrix}},}$ 

where ${\displaystyle S:=B^{*}A^{-1}B}$  denotes the Schur complement. Since ${\displaystyle S}$  is symmetric positive-definite, we can apply standard iterative methods like the gradient descent method or the conjugate gradient method to

${\displaystyle Sx_{2}=B^{*}A^{-1}b_{1}-b_{2}}$ 

in order to compute ${\displaystyle x_{2}}$ . The vector ${\displaystyle x_{1}}$  can be reconstructed by solving

${\displaystyle Ax_{1}=b_{1}-Bx_{2}.\,}$ 

It is possible to update ${\displaystyle x_{1}}$  alongside ${\displaystyle x_{2}}$  during the iteration for the Schur complement system and thus obtain an efficient algorithm.

We start the conjugate gradient iteration by computing the residual

${\displaystyle r_{2}:=B^{*}A^{-1}b_{1}-b_{2}-Sx_{2}=B^{*}A^{-1}(b_{1}-Bx_{2})-b_{2}=B^{*}x_{1}-b_{2},}$ 

of the Schur complement system, where

${\displaystyle x_{1}:=A^{-1}(b_{1}-Bx_{2})}$ 

denotes the upper half of the solution vector matching the initial guess ${\displaystyle x_{2}}$  for its lower half. We complete the initialization by choosing the first search direction

${\displaystyle p_{2}:=r_{2}.\,}$ 

In each step, we compute

${\displaystyle a_{2}:=Sp_{2}=B^{*}A^{-1}Bp_{2}=B^{*}p_{1}}$ 

and keep the intermediate result

${\displaystyle p_{1}:=A^{-1}Bp_{2}}$ 

for later. The scaling factor is given by

${\displaystyle \alpha :=p_{2}^{*}a_{2}/p_{2}^{*}r_{2}}$ 

and leads to the updates

${\displaystyle x_{2}:=x_{2}+\alpha p_{2},\quad r_{2}:=r_{2}-\alpha a_{2}.}$ 

Using the intermediate result ${\displaystyle p_{1}}$  saved earlier, we can also update the upper part of the solution vector

${\displaystyle x_{1}:=x_{1}-\alpha p_{1}.\,}$ 

Now we only have to construct the new search direction by the Gram–Schmidt process, i.e.,

${\displaystyle \beta :=r_{2}^{*}a_{2}/p_{2}^{*}a_{2},\quad p_{2}:=r_{2}-\beta p_{2}.}$ 

The iteration terminates if the residual ${\displaystyle r_{2}}$  has become sufficiently small or if the norm of ${\displaystyle p_{2}}$  is significantly smaller than ${\displaystyle r_{2}}$  indicating that the Krylov subspace has been almost exhausted.

If solving the linear system ${\displaystyle Ax=b}$  exactly is not feasible, inexact solvers can be applied.^[2][3][4]

If the Schur complement system is ill-conditioned, preconditioners can be employed to improve the speed of convergence of the underlying gradient method.^[2][5]

Inequality constraints can be incorporated, e.g., in order to handle obstacle problems.^[5]

• Chen, Zhangxin (2006). "Linear System Solution Techniques". Finite Element Methods and Their Applications. Berlin: Springer. pp. 145–154. ISBN 978-3-540-28078-1.