Skip to main content

Preconditioners for reduced saddle point systems arising in elliptic PDE-constrained optimization problems

Abstract

In this paper, we propose some preconditioning techniques for reduced saddle point systems arising from linear elliptic distributed optimal control problems. The eigenvalues of preconditioned matrices are analyzed. Moreover, the bounds of these eigenvalues with respect to the mesh size h are also obtained. Some numerical tests are presented to validate the theoretical analysis.

1 Introduction

We consider the preconditioning techniques for solving the saddle point system arising from linear elliptic distributed optimal control problems:

$$\begin{aligned}& \min_{u,f}\frac{1}{2}\| u-\hat{u} \|_{L^{2}(\Omega )}^{2}+\beta\| f \|_{L^{2}(\Omega)}^{2} \\& \quad\mbox{subject to } -\Delta u =f \quad \mbox{in } \Omega, \\& \hphantom{\quad\mbox{subject to }} u =g \quad \mbox{on } \partial\Omega, \end{aligned}$$
(1)

where \(\Omega\subset\mathbb{R}^{2}\) is a simply connected polygonal domain with a connected boundary ∂Ω. The data û is the target function, and the parameter β is the Tikhonov regular parameter. u is the state variable, and f is called the control variable. Such problems are simple models, which were originally introduced by Lions in [1]. Some more complex formulations, which include control constraints or state constraints, can be found in [2–4].

For the elliptic PDE-constrained optimization problem (1), we take the method based on discretize-then-optimize [5] approach. More concretely, we use the \(P_{1}\) conforming finite elements for the approximations of the state variable u and the control variable f, which yields the following corresponding minimization problem [5, 6]:

$$ \textstyle\begin{cases} \min\limits_{\mathbf{u},\mathbf{f}}\frac{1}{2}\mathbf{u}^{T}M\mathbf{u}- \mathbf {u}^{T}\mathbf{b}+\alpha+\beta\mathbf{f}^{T}M\mathbf{f}\\ \quad\mbox{subject to } K\mathbf{u}=M\mathbf{f}+\mathbf{d}, \end{cases} $$
(2)

where \(K\in\mathbb{R}^{m\times m}\) is the stiffness matrix, \(M\in \mathbb{R}^{m\times m}\) is the mass matrix. Both M and K are symmetric and positive definite. \(\mathbf{b}\in\mathbb{R}^{m}\) is the discrete Galerkin projection of the target function û. \(\alpha =\frac{1}{2}\|\hat{u}\|_{L^{2}(\Omega)}^{2}\), and \(\mathbf{d}\in\mathbb {R}^{m}\) contains the terms arising from the boundaries of finite element approximation of the state u. Then, using the Lagrange multiplier technique for the minimization problem (2), we can find that f, u and λ satisfy the following linear saddle point system [5]:

$$ A\mathbf{x}= \begin{pmatrix} 2\beta M &0&-M\\ 0&M&K\\ -M&K&0 \end{pmatrix} \begin{pmatrix} \mathbf{f}\\ \mathbf{u}\\ \lambda\end{pmatrix}= \begin{pmatrix} 0\\ \mathbf{b}\\ \mathbf{d}\end{pmatrix}=\mathbf{g}_{1}, $$
(3)

where λ is a vector of Lagrange multipliers. To obtain the above system, we have used the fact that \(K=K^{T}\). Obviously, from the first equation of (3), we can obtain that

$$ \mathbf{f}=\frac{1}{2\beta}\lambda. $$
(4)

Substituting (4) into (3), we can easily obtain the following reduced system:

$$ \mathcal{A}y= \begin{pmatrix} M &K\\ K&-\frac{1}{2\beta}M\end{pmatrix} \begin{pmatrix} \mathbf{u}\\ \lambda\end{pmatrix}= \begin{pmatrix} \mathbf{b}\\ \mathbf{d}\end{pmatrix}=\mathbf{g}. $$
(5)

This paper is devoted to constructing efficient preconditioners for the reduced linear saddle point system (5). Comparing with (3), (5) is a smaller linear system so that it can help us to reduce computational cost. As we know, previous work has been mainly devoted to the development of efficient solution techniques for the original linear saddle point system (3), such as the block diagonal preconditioner and the constraint preconditioner [5], block triangular diagonal preconditioner [7], block-counter-diagonal and block-counter-tridiagonal preconditioning techniques [6]. Other recent work in this direction can also be found in [8–11]. For the reduced saddle point problem (5), there have existed several preconditioning algorithms, we refer the reader to [12–18] for details. In our work, we extend similar preconditioning techniques, which were used in [5, 6] for system (3), to the saddle point system (5). We also refer the reader to [19] for analogous discussions, where the optimal control of Stokes equations was addressed. The major contribution of our work is to establish the bounds of eigenvalues of the preconditioned matrices with respect to the discretization mesh size h, though we also provide the analysis of eigenvalues of the preconditioned matrices using the usual linear algebraic techniques developed in [5, 6, 19]. In particular, for the case that the Tikhonov parameter β is not very small, the bounds of eigenvalues of preconditioned matrices are proved to remain bounded away from 0 as \(h\rightarrow0\) with a bound independent of h. This suggests that they are optimal preconditioners in the sense that their performances are independent of the mesh size h.

The rest of our paper is organized as follows. In Section 2, we provide preconditioning techniques for the situation that the Tikhonov parameter β is not very small, then we discuss the corresponding spectral analysis. In Section 3, we propose preconditioning techniques for the case that the Tikhonov parameter β is sufficiently small, and then the corresponding spectral analysis is presented. In Section 4, some numerical results are presented to demonstrate our theoretical analysis.

2 Preconditioners for the Tikhonov parameter β not very small

When the Tikhonov parameter β is not very small, we may provide block-counter diagonal preconditioner \(P_{\mathrm{bcd}}\) and block-counter-triangular diagonal preconditioner \(P_{\mathrm{bctd}}\) as the following form:

$$ P_{\mathrm{bcd}}= \begin{pmatrix} 0 & K\\ K &0\end{pmatrix} $$
(6)

and

$$ P_{\mathrm{bctd}}= \begin{pmatrix} M & K\\ K &0\end{pmatrix}. $$
(7)

For the block-counter diagonal preconditioner \(P_{\mathrm{bcd}}\), we have the following results.

Theorem 2.1

Let \(\mathcal{A} \in\mathbb{R}^{2m \times 2m}\) be the coefficient matrix of the linear system (5) and \(P_{\mathrm{bcd}} \in\mathbb{R}^{2m \times2m}\) be defined in (6). Denote by \(\nu_{l}\) the eigenvalue of \(K^{-1}M\in \mathbb{R}^{m \times m}\), where \(\nu_{l}>0\), \(l=1,2,\ldots,m\), and by \(x^{(l)} \in\mathbb{C}^{m}\) the corresponding eigenvector for \(l=1,2,\ldots,m\). Then

  1. (i)

    the eigenvalues of the preconditioned matrix \(P_{\mathrm {bcd}}^{-1}\mathcal{A}\) are

    $$ \lambda_{k}^{(l)}=1-\sqrt{\frac{1}{2\beta}} \nu_{l} e^{\frac{(2k+1) \pi i}{2}},\quad k=0,1, l=1,2, \ldots, m, $$

    where i denotes the imaginary unit;

  2. (ii)

    the eigenvectors of the preconditioned matrix \(P_{\mathrm {bcd}}^{-1}\mathcal{A}\) are of the form

    $$ \begin{pmatrix} {x}^{(l)}\\ -\frac{\sqrt{2\beta}}{\mu_{l} }e^{\frac{(2k+1) \pi i}{2}}K^{-1}M x^{(l)}\end{pmatrix},\quad l=1,2,\ldots,m ,k=0,1. $$

Proof

(i) By a direct calculation, we get

$$ P_{2}^{-1}\mathcal{A}= \begin{pmatrix} 0 &K^{-1}\\ K^{-1}&0\end{pmatrix} \begin{pmatrix} M &K\\ K&-\frac{1}{2\beta}M\end{pmatrix}= \begin{pmatrix} I &-\frac{1}{2\beta}K^{-1}M\\ K^{-1}M&I\end{pmatrix}. $$
(8)

To proceed, let

$$ E=I-P_{2}^{-1}\mathcal{A}= \begin{pmatrix} 0 &\frac{1}{2\beta}K^{-1}M\\ -K^{-1}M&0\end{pmatrix}= \begin{pmatrix} 0 &\frac{1}{2\beta}D\\ -D&0\end{pmatrix}, $$
(9)

where \(D=K^{-1}M\). Therefore,

$$ E^{2}= \begin{pmatrix} -\frac{1}{2\beta}D^{2} &0\\ 0&-\frac{1}{2\beta}D^{2}\end{pmatrix}. $$
(10)

Then from (10) we can know that if \(\nu_{l}\) is the eigenvalue of \(K^{-1}M\), then the eigenvalues of E are \(\sqrt{\frac{1}{2\beta}} \nu_{l}e^{\frac{(2k+1) \pi i}{2}}\), \(k=0,1\). Thus, the eigenvalues of the preconditioned matrix \(P_{\mathrm{bcd}}^{-1}\mathcal{A}\) are \(1-\sqrt{\frac{1}{2\beta}}\nu _{l}e^{\frac{(2k+1) \pi i}{2}}\), \(k=0,1\). This yields (i).

(ii) By similar techniques used in the proof of (ii) in Theorem 3.1, we can obtain the corresponding conclusions. □

For the block-counter-triangular diagonal preconditioner \(P_{\mathrm {bctd}}\), we have the following results.

Theorem 2.2

Let \(\mathcal{A} \in\mathbb{R}^{2m \times2m}\) be the coefficient matrix of the linear system (5) and \(P_{\mathrm{bctd}} \in\mathbb{R}^{2m \times2m}\) be defined in (7). Denote by \(\nu_{l}\) the eigenvalue of \(K^{-1}M\in\mathbb{R}^{m \times m}\), where \(\nu_{l}>0\), \(l=1,2,\ldots,m\), and by \(y^{(l)} \in\mathbb{C}^{m}\) the eigenvector of \((K^{-1}M)^{2}\) for \(l=1,2,\ldots,m\). Then

  1. (i)

    the eigenvalues of the preconditioned matrix \(P_{\mathrm {bctd}}^{-1}\mathcal{A}\) are 1 with algebraic multiplicity m, and the remaining eigenvalues equal \(\frac{1}{2\beta}\nu_{l}^{2}+1\), \(l=1,2,\ldots,m\);

  2. (ii)

    the eigenvectors of the preconditioned matrix \(P_{\mathrm {bctd}}^{-1}\mathcal{A}\) are of the form

    $$ \begin{pmatrix}x\\ 0\end{pmatrix}, \quad\forall x\in{C}^{m}\backslash\{0\},\quad \textit{and}\quad \begin{pmatrix} -\frac{1}{\nu}K^{-1}My^{(l)}\\ y^{(l)}\end{pmatrix},\quad l=1,2,\ldots,m. $$

Proof

(i) A straightforward calculation shows that

$$ P_{\mathrm{bctd}}^{-1}\mathcal{A}= \begin{pmatrix} 0 &K^{-1}\\ K^{-1}&-K^{-1}MK^{-1}\end{pmatrix} \begin{pmatrix} M &K\\ K&-\frac{1}{2\beta}M\end{pmatrix}= \begin{pmatrix} I &-\frac{1}{2\beta}K^{-1}M\\ 0& I+\frac{1}{2\beta}K^{-1}MK^{-1}M\end{pmatrix}. $$
(11)

It follows from the above equality that the eigenvalues of the preconditioned matrix \(P_{\mathrm{bctd}}^{-1}\mathcal{A}\) are 1 with algebraic multiplicity m, and m eigenvalues with the form \(1+\frac {1}{2\beta}\nu^{2}\), where ν is any eigenvalue of the matrix \(K^{-1}M\). This demonstrates the conclusions of (i).

(ii) Let \(\{\lambda, (x,y)^{T}\}\) be an eigenpair for the matrix \(P_{\mathrm {bctd}}^{-1}\mathcal{A}\), then from the expression of \(P_{\mathrm {bctd}}^{-1}\mathcal{A}\) in (11), we have

$$ \begin{pmatrix} I&-\frac{1}{2\beta}K^{-1}M\\ 0&I+\frac{1}{2\beta}K^{-1}MK^{-1}M\end{pmatrix} \begin{pmatrix} x\\ y\end{pmatrix}= \begin{pmatrix} x\\ y\end{pmatrix}=\lambda\mathbf{e}, $$

or the above equation can be rewritten as

$$ \textstyle\begin{cases} (x-\frac{1}{2\beta}K^{-1}My)=\lambda x,\\ (I+\frac{1}{2\beta}K^{-1}MK^{-1}M)y=\lambda y. \end{cases} $$

It is easy to see that if \(\lambda=1\), then \(y=0\), thus the corresponding eigenvectors take the form

$$ \mathbf{e}= \begin{pmatrix}x\\ 0\end{pmatrix}, \quad\mbox{with } \forall x\in\mathbb{C}^{m} \backslash\{0\}. $$
(12)

If \(\lambda\neq1\), we have

$$ x=\frac{1}{2\beta(1-\lambda)}K^{-1}My $$
(13)

and

$$ \bigl(I+2\beta M^{-1}KM^{-1}K\bigr)y=\lambda y. $$
(14)

Thus, \(\lambda=\frac{1}{2\beta}\nu^{2}+1\), with ν an eigenvalue of the matrix \(K^{-1}M\). Direct computations lead to

$$ x=-\frac{1}{\nu}K^{-1}My. $$
(15)

Obviously, \(y\in\mathbb{C}^{m}\) is an eigenvector of the matrix \((K^{-1}M)^{2}\). This completes the proof. □

Furthermore, using the theory of finite element methods, we can use the results stated in the above two theorems to obtain more concrete bounds that involved the discretization mesh size h. Since in our numerical tests we use the \(P_{1}\) conforming finite elements to discretize problem (1), we recall the following lemma from [5, 20] which will be used frequently in our theory analysis.

Lemma 2.3

Let \(\mathcal{T}_{h}=\{T_{i}: i=1,2,\ldots,E_{h}\} \) be a family of quasi-uniform and shape regularity triangulations of the domain \(\Omega\subset\mathbb{R}^{2}\). Denote by \(h_{T_{i}}\) the diameter of the element \(T_{i}\) and \(h=\max_{T_{i}\in\mathcal {T}_{h}}h_{T_{i}}\). For the \(P_{1}\) conforming finite element, the mass matrix M and the stiffness matrix K satisfy

$$\begin{aligned}& c_{1}h^{2} \leq \frac{\mathbf{x}^{T}M\mathbf {x}}{\mathbf{x}^{T}\mathbf{x}}\leq c_{2}h^{2}, \end{aligned}$$
(16)
$$\begin{aligned}& d_{1}h^{2} \leq \frac{\mathbf{x}^{T}K\mathbf {x}}{\mathbf{x}^{T}\mathbf{x}}\leq d_{2}, \end{aligned}$$
(17)

\(\forall\mathbf{x}\neq0\in\mathbb{R}^{m}\). Constants \(c_{1}\), \(c_{2}\), \(d_{1}\) and \(d_{2}\) are positive and independent of h and β.

With the above prepared techniques, we can obtain the following two corollaries.

Corollary 2.4

Let λ be an eigenvalue of \(P_{\mathrm{bcd}}^{-1}\mathcal{A}\), with \(P_{\mathrm{bcd}}\) and \(\mathcal{A}\) being defined in (6) and (5). Then the following bound holds:

$$ \sqrt{1+\frac{c_{1}^{2}h^{4}}{2\beta d_{2}^{2}}} \leq|\lambda| \leq\sqrt {1+\frac {c_{2}^{2}}{2\beta d_{1}^{2}}}. $$
(18)

Proof

From (i) in Theorem 2.1, the eigenvalues of \(P_{\mathrm {bcd}}^{-1}\mathcal{A}\) are

$$ \lambda_{k}^{(l)}=1-\sqrt{\frac{1}{2\beta}} \nu_{l} e^{\frac{(2k+1) \pi i}{2}},\quad k=0,1, l=1,2, \ldots, m. $$

A direct calculation shows that

$$ \bigl|\lambda_{k}^{(l)}\bigr|^{2}=1+ \frac{1}{2\beta}\nu_{l}^{2}-2\sqrt{\frac{1}{2\beta }}\cos \frac{(2k+1)\pi}{2}=1+\frac{1}{2\beta}\nu_{l}^{2}. $$
(19)

In the last step of the above equation we have used the fact that \(\cos\frac{(2k+1)\pi}{2}=0\).

Thus, from the above equality, we know that if we want to bound \(|\lambda_{k}^{(l)}|\), it is sufficient to bound \(\nu_{l}\). Recalling that \(\nu_{l}\) is an eigenvalue of \(K^{-1}M\), we get

$$\begin{aligned}& K^{-1}M\mathbf{x} = \nu_{l}\mathbf{x}, \\& \textit{i.e.},\ M\mathbf{x} = \nu _{l}K\mathbf{x}, \\& \mbox{thus } \mathbf{x}^{T}M\mathbf{x} = \nu _{l} \mathbf{x}^{T}K\mathbf{x}, \\& \nu_{l} = \frac{\mathbf{x} ^{T}M\mathbf{x}}{\mathbf{x}^{T}K\mathbf{x}} . \end{aligned}$$
(20)

From (17) in Lemma 2.3, we obtain

$$ \frac{1}{d_{2}}\leq\frac{\mathbf{x}^{T}\mathbf{x}}{\mathbf{x}^{T}K\mathbf {x}}\leq\frac{1}{d_{1}h^{2}}. $$
(21)

Thus, from (16), (20) and (21) we get

$$ \frac{c_{1}h^{2}}{d_{2}}\leq\nu_{l}\leq\frac{c_{2}}{d_{1}}. $$
(22)

Then the corollary follows by (19) and (22). □

Corollary 2.5

Let λ be an eigenvalue of \(P_{\mathrm{bctd}}^{-1}\mathcal{A}\), with \(P_{\mathrm{bctd}}\) and \(\mathcal{A}\) defined in (7) and (5). Then the eigenvalues of the preconditioned matrix \(P_{\mathrm{bctd}}^{-1}\mathcal{A}\) are 1 with algebraic multiplicity m, the remaining m eigenvalues can be bounded as follows:

$$ \sqrt{1+\frac{c_{1}^{2}h^{4}}{2\beta d_{2}^{2}}} \leq\lambda\leq\sqrt {1+\frac {c_{2}^{2}}{2\beta d_{1}^{2}}}. $$
(23)

Proof

From (i) in Theorem 2.2, using the same proof procedure as in Corollary 2.4, we can obtain the corresponding conclusion. □

Remark 2.6

From (18) and (23) in Corollaries 2.4 and 2.5, we can see that the eigenvalues of the preconditioned matrices \(P_{\mathrm{bcd}}^{-1}\mathcal{A}\) and \(P_{\mathrm {bctd}}^{-1}\mathcal{A}\) satisfy

$$ 1 \leq|\lambda| \leq\sqrt{1+\frac{c_{2}^{2}}{2\beta d_{1}^{2}}}\quad \mbox{as } h\rightarrow0. $$
(24)

Furthermore, if β is not too small, then from (24) we know that the eigenvalues of the preconditioned matrices \(P_{\mathrm {bcd}}^{-1}\mathcal{A}\) and \(P_{\mathrm{bctd}}^{-1}\mathcal{A}\) are clustered. When the parameter β is not too small, (24) also means that \(P_{\mathrm{bcd}}\) and \(P_{\mathrm{btcd}}\) may perform as optimal preconditioners in the sense that their performances are independent of the mesh size h.

Remark 2.7

In practical computations, we may choose good approximations K̃ and M̃ for K and M, respectively. Then the corresponding preconditioners \(P_{\mathrm{bcd}}\) and \(P_{\mathrm{bctd}}\) are replaced by

$$ \tilde{P}_{\mathrm{bcd}}= \begin{pmatrix} 0 &\tilde{K}\\ \tilde{K}&0 \end{pmatrix} \quad\mbox{and}\quad \tilde{P}_{\mathrm{bctd}}= \begin{pmatrix} \tilde{M} &\tilde{K}\\ \tilde{K}&0\end{pmatrix}, $$

respectively.

3 Preconditioners for the Tikhonov parameter β sufficiently small

In this section, we discuss block diagonal and block-triangular diagonal preconditioning techniques for the linear system (5) and analyze the eigenvalues and eigenvectors of the preconditioned matrices. The block diagonal preconditioner is given by

$$ P_{\mathrm{bd}}= \begin{pmatrix} M &0\\ 0&-\frac{1}{2\beta}M\end{pmatrix}, $$
(25)

and the block-triangular diagonal preconditioner is of the form

$$ P_{\mathrm{btd}}= \begin{pmatrix} M &K\\ 0&-\frac{1}{2\beta}M\end{pmatrix}. $$
(26)

Moreover, we find that both the preconditioned matrices \(P_{\mathrm {bd}}^{-1}\mathcal{A}\) and \(P_{\mathrm{btd}}^{-1}\mathcal{A}\) have eigenvalues around 1 when the Tikhonov parameter β is sufficiently small (i.e., \(\beta\ll1\)). These properties are stated in the following two theorems. The spectral analysis for preconditioned matrices is largely based on the references [5, 19]. First, we consider the block diagonal preconditioner.

Theorem 3.1

Let \(\mathcal{A} \in\mathbb{R}^{2m \times2m}\) be the coefficient matrix of the linear system (5) and \(P_{\mathrm{bd}} \in\mathbb{R}^{2m \times2m}\) be defined in (25). Set \(\mu_{l}\) the eigenvalue \(M^{-1}K\in\mathbb {R}^{m \times m}\), where \(\mu_{l}>0\), \(l=1,2,\ldots,m\), and \(x^{(l)} \in \mathbb{C}^{m}\) the corresponding eigenvector for \(l=1,2,\ldots,m\). Then

  1. (i)

    the eigenvalues of the preconditioned matrix \(P_{\mathrm {bd}}^{-1}\mathcal{A}\) are

    $$ \lambda_{k}^{(l)}=1-\sqrt{2\beta}\mu_{l} e^{\frac{(2k+1) \pi i}{2}},\quad k=0,1, l=1,2, \ldots, m, $$

    where i denotes the imaginary unit;

  2. (ii)

    the eigenvectors of the preconditioned matrix \(P_{\mathrm {bd}}^{-1}\mathcal{A}\) are of the form

    $$ \begin{pmatrix} -\frac{1}{\sqrt{2\beta}\mu_{l} e^{\frac{(2k+1) \pi i}{2}}}M^{-1}K x^{(l)}\\ x^{(l)}\end{pmatrix},\quad l=1,2,\ldots,m, k=0,1. $$

Proof

(i) A simple calculation shows that

$$ P_{\mathrm{bd}}^{-1}\mathcal{A}= \begin{pmatrix} M^{-1} &0\\ 0&-{2\beta}M^{-1}\end{pmatrix} \begin{pmatrix} M &K\\ K&-\frac{1}{2\beta}M\end{pmatrix} = \begin{pmatrix} I &M^{-1}K\\ -2\beta M^{-1}K&I\end{pmatrix}. $$
(27)

To proceed, let

$$ C=I-P_{\mathrm{bd}}^{-1}\mathcal{A}= \begin{pmatrix} 0 &-M^{-1}K\\ 2\beta M^{-1}K&0 \end{pmatrix}= \begin{pmatrix} 0 &-B\\ 2\beta B&0\end{pmatrix}, $$
(28)

where \(B=M^{-1}K\). Therefore,

$$ C^{2}= \begin{pmatrix} -2\beta B^{2} &0\\ 0&-2\beta B^{2}\end{pmatrix}. $$
(29)

From (29) we can know that if \(\mu_{l}\) is the eigenvalue of \(M^{-1}K\), then the eigenvalues of C are \(\sqrt{2\beta}\mu_{l} e^{\frac {(2k+1)\pi i}{2}}\), \(k=0,1\). Thus, the eigenvalues of the preconditioned matrix \(P_{\mathrm {bd}}^{-1}\mathcal{A}\) are \(1-\sqrt{2\beta}\mu_{l} e^{\frac{(2k+1)\pi i}{2}}\), \(k=0,1\). This leads to (i).

(ii) Let \(\mu_{l}\) be the eigenvalue of \(M^{-1}K\in\mathbb{R}^{m\times m}\) and \(x^{(l)}\in\mathbb{C}^{m}\) be the corresponding eigenvector for \(l=1,2,\ldots,m\), then

$$ M^{-1}KM^{-1}Kx^{(l)}=\mu_{l}^{2}x^{(l)}. $$

Multiplying both sides of the above equation by \(-2\beta\) yields

$$ -2\beta M^{-1}KM^{-1}Kx^{(l)}=-2\beta \mu_{l}^{2}x^{(l)}= \bigl(\sqrt{2\beta } \mu_{l} e^{\frac{(2k+1) \pi i}{2}} \bigr)^{2} x^{(l)}, $$

i.e.,

$$ -2\beta M^{-1}K\frac{1}{\sqrt{2\beta}\mu_{l} e^{\frac{(2k+1) \pi i}{2}}} \cdot M^{-1}Kx^{(l)}= \bigl(\sqrt{2\beta}\mu_{l} e^{\frac{(2k+1) \pi i}{2}} \bigr) x^{(l)}. $$
(30)

Let

$$ y=-\frac{1}{\sqrt{2\beta}\mu_{l} e^{\frac{(2k+1) \pi i}{2}}}M^{-1}K x^{(l)}, $$

then

$$ -M^{-1}Kx^{(l)}= \bigl(\sqrt{2\beta} \mu_{l} e^{\frac{(2k+1) \pi i}{2}} \bigr) y. $$
(31)

From the expressions of (30) and (31) we obtain

$$ \textstyle\begin{cases} -Bx^{(l)}= (\sqrt{2\beta}\mu_{l} e^{\frac{(2k+1) \pi i}{2}} ) y,\\ 2\beta By= (\sqrt{2\beta}\mu_{l} e^{\frac{(2k+1) \pi i}{2}} )x^{(l)} . \end{cases} $$

Thus,

$$ C \begin{pmatrix} y\\ x^{(l)}\end{pmatrix}= \begin{pmatrix} 0& -B\\ 2\beta B&0 \end{pmatrix} \begin{pmatrix} y\\ x^{(l)}\end{pmatrix} = \bigl(\sqrt{2\beta}\mu_{l} e^{\frac{(2k+1) \pi i}{2}} \bigr) \begin{pmatrix} y\\ x^{(l)}\end{pmatrix}. $$

Therefore, if \(x^{(l)}\in\mathbb{C}^{m}\) is an eigenvector of the matrix B, then \({ y\choose x^{(l)} } \) is the eigenvector of the matrix \(C\in\mathbb{R}^{{2m}\times{2m}}\). This finishes the proof. □

Then we turn our attention to the block-triangular diagonal preconditioner.

Theorem 3.2

Let \(\mathcal{A} \in\mathbb{R}^{2m \times2m}\) be the coefficient matrix of the linear system (5) and \(P_{\mathrm{btd}} \in\mathbb{R}^{2m \times2m}\) be defined in (26). Set \(\mu_{l}\) the eigenvalue of \(M^{-1}K\in \mathbb{R}^{m \times m}\), where \(\mu_{l}>0\), \(l=1,2,\ldots,m\), and \(x^{(l)} \in\mathbb{C}^{m}\) the eigenvector of the matrix \((M^{-1}K)^{2}\) for \(l=1,2,\ldots,m\). Then

  1. (i)

    the eigenvalues of the preconditioned matrix \(P_{\mathrm {btd}}^{-1}\mathcal{A}\) are 1 with algebraic multiplicity m, and the remaining eigenvalues are equal to \(2\beta\mu_{l}^{2}+1\), \(l=1,2,\ldots,m\);

  2. (ii)

    the eigenvectors of the preconditioned matrix \(P_{\mathrm {btd}}^{-1}\mathcal{A}\) are of the form

    $$ \begin{pmatrix}0\\ y\end{pmatrix}, \quad\forall y\in{C}^{m}\backslash\{0\}, \quad\textit{and} \quad \begin{pmatrix}x^{(l)}\\ -\frac{1}{\mu}M^{-1}Kx^{(l)}\end{pmatrix}, \quad l=1,2,\ldots,m. $$

Proof

(i) By a direct calculation we get that

$$ P_{\mathrm{btd}}^{-1}\mathcal{A}= \begin{pmatrix} M^{-1} &2\beta M^{-1}KM^{-1}\\ 0&-{2\beta}M^{-1}\end{pmatrix} \begin{pmatrix} M &K\\ K&-\frac{1}{2\beta}M \end{pmatrix} = \begin{pmatrix} I+2\beta M^{-1}KM^{-1}K &0\\ -2\beta M^{-1}K&I \end{pmatrix} . $$
(32)

It follows from the above equality that the eigenvalues of the preconditioned matrix \(P_{\mathrm{btd}}^{-1}\mathcal{A}\) are 1 with algebraic multiplicity m, and m eigenvalues with the form \(1+2\beta \mu^{2}\), where μ is any eigenvalue of the matrix \(M^{-1}K\). This demonstrates the conclusions of (i).

(ii) Let \(\{\lambda, (x,y)^{T}\}\) be an eigenpair for the matrix \(P_{\mathrm {btd}}^{-1}\mathcal{A}\), then from the expression of \(P_{\mathrm{btd}}^{-1}\mathcal{A}\) in (32), we have

$$ \begin{pmatrix} I+2\beta M^{-1}KM^{-1}K&0\\ -2\beta M^{-1}K&I\end{pmatrix} \begin{pmatrix} x\\ y\end{pmatrix} =\lambda \begin{pmatrix} x\\ y\end{pmatrix} =\lambda \mathbf{e} , $$

or the above equation can be rewritten as

$$ \textstyle\begin{cases} (I+2\beta M^{-1}KM^{-1}K)x=\lambda x,\\ -2\beta M^{-1}Kx+y =\lambda y. \end{cases} $$

It is easy to see that if \(\lambda=1\), then \(x=0\), thus the corresponding eigenvectors take the form

$$ \mathbf{e}= \begin{pmatrix}0\\ y\end{pmatrix} , \quad\mbox{with } \forall y\in \mathbb{C}^{m}\backslash\{0\}. $$
(33)

If \(\lambda\neq1\), we have

$$ y=\frac{2\beta}{1-\lambda}M^{-1}Kx $$
(34)

and

$$ \bigl(I+2\beta M^{-1}KM^{-1}K\bigr)x=\lambda x. $$
(35)

Thus, \(\lambda=2\beta\mu^{2}+1\), with μ an eigenvalue of the matrix \(M^{-1}K\). A direct calculation shows that

$$ y=-\frac{1}{\mu}M^{-1}Kx. $$
(36)

Obviously, \(x\in\mathbb{C}^{m}\) is an eigenvector of the matrix \((M^{-1}K)^{2}\). This completes the proof. □

Furthermore, we provide more concrete bounds which involved the discretization mesh size h. The results are stated in the following corollaries.

Corollary 3.3

Let λ be an eigenvalue of \(P_{\mathrm{bd}}^{-1}\mathcal{A}\), with \(P_{\mathrm{bd}}\) and \(\mathcal{A}\) defined in (25) and (5). Then the following bound holds:

$$ \sqrt{1+\frac{2\beta d_{1}^{2}}{c_{2}^{2}}} \leq|\lambda| \leq\sqrt {1+\frac {2\beta d_{2}^{2}}{c_{1}^{2}h^{4}}}. $$
(37)

Proof

From (i) in Theorem 3.1, the eigenvalues of \(P_{\mathrm {bd}}^{-1}\mathcal{A}\) are

$$ \lambda_{k}^{(l)}=1-\sqrt{{2\beta}}\mu_{l} e^{\frac{(2k+1) \pi i}{2}},\quad k=0,1, l=1,2, \ldots, m. $$

A direct calculation leads to

$$ \bigl|\lambda_{k}^{(l)}\bigr|^{2}=1+{2\beta} \mu_{l}^{2}-2\sqrt{2\beta}\cos\frac {(2k+1)\pi}{2}=1+{2\beta} \mu_{l}^{2}. $$
(38)

In the last step of the above equation we have used the fact that \(\cos\frac{(2k+1)\pi}{2}=0\).

Thus, from the above equation, we know that if we want to bound \(|\lambda_{k}^{(l)}|\), it is enough to bound \(\mu_{l}\). Recalling that \(\mu_{l}\) is an eigenvalue of \(M^{-1}K\), we get

$$\begin{aligned}& M^{-1}K\mathbf{x} = \mu_{l}\mathbf{x} , \\& \textit{i.e.},\ K\mathbf{x} = \mu _{l}M\mathbf{x} , \\& \mbox{thus } \mathbf{x}^{T}K\mathbf{x} = \mu _{l} \mathbf{x}^{T}M\mathbf{x} , \\& \mu_{l} = \frac{\mathbf{x}^{T}K\mathbf{x}}{ \mathbf{x}^{T}M\mathbf{x}}. \end{aligned}$$
(39)

From (16) in Lemma 2.3, we obtain

$$ \frac{1}{c_{2}h^{2}}\leq\frac{\mathbf{x}^{T}\mathbf{x}}{\mathbf{x}^{T}M\mathbf {x}}\leq\frac{1}{c_{1}h^{2}}. $$
(40)

Thus, from (17), (39) and (40) we get

$$ \frac{d_{1}}{c_{2}}\leq\mu_{l}\leq\frac{d_{2}}{c_{1}h^{2}}. $$
(41)

Then the corollary follows by (38) and (41). □

Corollary 3.4

Let λ be an eigenvalue of \(P_{\mathrm{btd}}^{-1}\mathcal{A}\), with \(P_{\mathrm{btd}}\) and \(\mathcal{A}\) defined in (26) and (5). Then the eigenvalues of the preconditioned matrix \(P_{\mathrm{btd}}^{-1}\mathcal{A}\) are 1 with algebraic multiplicity m, the remaining m eigenvalues can be bounded as follows:

$$ \sqrt{1+\frac{2\beta d_{1}^{2}}{c_{2}^{2}}} \leq\lambda\leq\sqrt {1+\frac {2\beta d_{2}^{2}}{c_{1}^{2}h^{4}}}. $$
(42)

Proof

From (i) in Theorem 3.2, using the same proof procedure as in Corollary 3.3, we can obtain the corresponding conclusion. □

Remark 3.5

From Theorems 3.1 and 3.2, we can see that the eigenvalues of the preconditioned matrices \(P_{\mathrm{bd}}^{-1}\mathcal{A}\) and \(P_{\mathrm{btd}}^{-1}\mathcal{A}\) are around 1 when the Tikhonov parameter β is sufficiently small.

Remark 3.6

In actual computations, we may choose good approximations K̃ and M̃ for K and M, respectively. Then the corresponding preconditioners \(P_{\mathrm{bd}}\) and \(P_{\mathrm{btd}}\) are replaced by

$$ \tilde{P}_{\mathrm{bd}}= \begin{pmatrix} \tilde{M} &0\\ 0&-\frac{1}{2\beta}\tilde{M} \end{pmatrix} \quad\mbox{and}\quad \tilde{P}_{\mathrm{btd}}= \begin{pmatrix} \tilde{M} &\tilde{K}\\ 0&-\frac{1}{2\beta}\tilde{M} \end{pmatrix}, $$

respectively. Some examples of proper M̃ and K̃ were presented in [5, 7].

4 Numerical experiments

The test problem we consider is the following linear elliptic optimal control problem:

$$\begin{aligned}& \min_{u,f}\frac{1}{2} \| u-\hat{u}\|_{2}^{2}+\beta \| f \|_{2}^{2} \\& \quad\mbox{subject to } -\Delta u=f \quad\mbox{in } \Omega, \\& \hphantom{\quad\mbox{subject to }} u=0 \quad\mbox{on } \partial\Omega, \end{aligned}$$

where \(\Omega=(0,1)\times(0,1)\), \(\hat{u}\in L^{2}(\Omega)\) is given by

$$ \hat{u}= \textstyle\begin{cases} 1 &\mbox{if } (x,y)\in[0,\frac{1}{2}]^{2},\\ 0, &\mbox{otherwise}. \end{cases} $$

Our numerical experiments are performed by MATLAB. We consider four uniformly refined meshes, which are constructed by subsequently splitting each triangle into four triangles by connecting the midpoints of the edges of the triangle. First, by setting the regular parameter \(\beta=10^{-2}\) (not very small), we show the eigenvalues of the preconditioned matrices \(P_{\mathrm{bcd}}^{-1}\mathcal{A}\) and \(P_{\mathrm{bctd}}^{-1}\mathcal {A}\) in Figures 1 and 2, they are both clustered. These results demonstrate the theoretical analysis in Theorems 2.1 and 2.2. Also, by setting the regular parameter \(\beta=10^{-8}\) (very small), we show the eigenvalues of the preconditioned matrices \(P_{\mathrm{bd}}^{-1}\mathcal{A}\) and \(P_{\mathrm{btd}}^{-1}\mathcal {A}\) in Figures 3 and 4, both of them are clustered. These results validate the theoretical analysis in Theorems 3.1 and 3.2.

Figure 1
figure 1

The eigenvalue distribution of the preconditioned matrix \(\pmb{P_{\mathrm{bcd}}^{-1}\mathcal{A}}\) with \(\pmb{\beta=10^{-2}}\) and \(\pmb{h=2^{-3}}\) .

Figure 2
figure 2

The eigenvalue distribution of the preconditioned matrix \(\pmb{P_{\mathrm{bctd}}^{-1}\mathcal{A}}\) with \(\pmb{\beta=10^{-2}}\) and \(\pmb{h=2^{-3}}\) .

Figure 3
figure 3

The eigenvalue distribution of the preconditioned matrix \(\pmb{P_{\mathrm{bd}}^{-1}\mathcal{A}}\) with \(\pmb{\beta=10^{-8}}\) and \(\pmb{h=2^{-3}}\) .

Figure 4
figure 4

The eigenvalue distribution of the preconditioned matrix \(\pmb{P_{\mathrm{btd}}^{-1}\mathcal{A}}\) with \(\pmb{\beta=10^{-8}}\) and \(\pmb{h=2^{-3}}\) .

Furthermore, for comparison, we show the number of iterations for different preconditioners with different regular parameters in Table 1. In all implementations, we use zero as initial guess and stop the iteration when \(\| r^{(k)}(=b-\mathcal {A}x^{(k)})\|_{2}/\| r^{(0)}(=b)\|_{2} \leq 10^{-7}\). The iteration methods we use are preconditioned GMRES methods. From the results stated in Table 1, first we find that when the regular parameter β is not very small, the number of iterations of \(P_{\mathrm{bcd}}\) and \(P_{\mathrm{bctd}}\) preconditioners is smaller than \(P_{\mathrm{bd}}\) and \(P_{\mathrm{btd}}\), thus we should choose \(P_{\mathrm{bcd}}\) or \(P_{\mathrm{bctd}}\) as preconditioners. Moreover, in this case, we observe that the number of iterations of the preconditioned GMRES methods is hardly sensitive to the changes in the mesh size h, this validates the analysis in Remark 2.6. On the other hand, when the regular parameter β is very small, the number of iterations of \(P_{\mathrm{bd}}\) and \(P_{\mathrm{btd}}\) preconditioners is smaller than \(P_{\mathrm{bcd}}\) and \(P_{\mathrm {bctd}}\), thus \(P_{\mathrm{bd}}\) or \(P_{\mathrm{btd}}\) preconditioners are our choices.

Table 1 The number of iterations for different preconditioners

References

  1. Lions, JL: Optimal Control of Systems. Springer, New York (1968)

    Google Scholar 

  2. Bergounioux, M, Ito, K, Kunisch, K: Primal-dual strategy for constrained optimal control problems. SIAM J. Control Optim. 37, 1176-1194 (1999)

    Article  MATH  MathSciNet  Google Scholar 

  3. Bergounioux, M, Kunisch, K: Primal-dual strategy for state-constrained optimal control problems. Comput. Optim. Appl. 22, 193-224 (2002)

    Article  MATH  MathSciNet  Google Scholar 

  4. Casas, E: Control of an elliptic problem with pointwise state constraints. SIAM J. Control Optim. 31, 1297-1327 (1993)

    Article  MathSciNet  Google Scholar 

  5. Rees, T, Dollar, HS, Wathen, AJ: Optimal solvers for PDE-constrained optimization. SIAM J. Sci. Comput. 32, 271-298 (2010)

    Article  MATH  MathSciNet  Google Scholar 

  6. Bai, Z: Block preconditioners for elliptic PDE-constrained problems. Computing 91, 379-395 (2011)

    Article  MATH  MathSciNet  Google Scholar 

  7. Rees, T, Stoll, M: Block-triangular preconditioners for PDE-constrained optimization. Numer. Linear Algebra Appl. 17, 977-996 (2010)

    Article  MATH  MathSciNet  Google Scholar 

  8. Herzog, R, Sachs, E: Preconditioned conjugate gradient method for optimal control problems with control and state constraints. SIAM J. Matrix Anal. Appl. 31, 2291-2317 (2010)

    Article  MATH  MathSciNet  Google Scholar 

  9. Pearson, JW, Wathen, A: A new approximation of the Schur complement in preconditioners for PDE-constrained optimization. Numer. Linear Algebra Appl. 19, 816-829 (2012)

    Article  MATH  MathSciNet  Google Scholar 

  10. Schöberl, J, Zulehner, W: Symmetric indefinite preconditioners for saddle point problems with applications to PDE-constrained optimization problems. SIAM J. Matrix Anal. Appl. 29, 752-773 (2007)

    Article  MATH  MathSciNet  Google Scholar 

  11. Zhang, G, Zheng, Z: Block-symmetric and block-lower-triangular preconditioner PDE-constrained optimization problems. J. Comput. Math. 31, 370-381 (2013)

    Article  MATH  MathSciNet  Google Scholar 

  12. Bai, Z, Benzi, M, Chen, F, Wang, Z: Preconditioned MHSS iteration methods for a class of block two-by-two linear systems with applications to distributed control problems. IMA J. Numer. Anal. 33, 343-369 (2013)

    Article  MATH  MathSciNet  Google Scholar 

  13. Borzì, A, Schulz, V: Multigrid methods for PDE optimization. SIAM Rev. 51, 361-395 (2009)

    Article  MATH  MathSciNet  Google Scholar 

  14. Lass, O, Vallejos, M, Borzi, A, Douglas, CC: Implementation and analysis of multigrid schemes with finite elements for elliptic optimal control problems. Computing 84, 27-48 (2009)

    Article  MATH  MathSciNet  Google Scholar 

  15. Schöberl, J, Simon, R, Zulehner, W: A robust multigrid method for elliptic optimal control problems. SIAM J. Numer. Anal. 49, 1482-1503 (2011)

    Article  MATH  MathSciNet  Google Scholar 

  16. Simoncini, V: Reduced order solution of structured linear systems arising in certain PDE-constrained optimization problems. Comput. Optim. Appl. 53, 591-617 (2012)

    Article  MATH  MathSciNet  Google Scholar 

  17. Vallejos, M, Borzì, A: Multigrid optimization methods for linear and bilinear elliptic optimal control problems. Computing 82, 31-52 (2008)

    Article  MATH  MathSciNet  Google Scholar 

  18. Zulehner, W: Nonstandard norms and robust estimates for saddle point problems. SIAM J. Matrix Anal. Appl. 32, 536-560 (2011)

    Article  MATH  MathSciNet  Google Scholar 

  19. Wang, S: A class of preconditioned iterative methods for the optimal control of the Stokes equations. Master thesis, Nanjing normal university (May 2014)

  20. Elman, HC, Silvester, DJ, Wathen, AJ: Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics. Numerical Mathematics and Scientific Computation Series. Oxford University Press, Oxford (2005)

    Google Scholar 

Download references

Acknowledgements

We thank the anonymous referees for their valuable comments and suggestions which led to an improved presentation of this paper. This work was supported by the opening fund of Jiangsu Key Lab for NSLSCS (Grant No. 201402), the Training Program for Outstanding Young Teachers in Guangdong Province (Grant No. 20140202), and the Educational Commission of Guangdong Province (Grant No. 2014KQNCX210).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Yuping Zeng.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

The main idea of this paper was proposed by YZ. All authors contributed equally in writing this article and read and approved the final manuscript.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zeng, Y., Wang, S., Xu, H. et al. Preconditioners for reduced saddle point systems arising in elliptic PDE-constrained optimization problems. J Inequal Appl 2015, 355 (2015). https://doi.org/10.1186/s13660-015-0879-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13660-015-0879-x

MSC

Keywords