Block $\omega$-circulant preconditioners for parabolic optimal control problems

In recent years, there has been a growing interest in analyzing and solving optimization problems constrained by parabolic equations. We refer to [25, 14, 34, 6] and the references therein for a comprehensive overview.

In this work, we are interested in solving the distributed optimal control model problem. Namely, the following quadratic cost functional is minimized:

$$\min_{y,u}~{}\mathcal{J}(y,u):=\frac{1}{2}\|y-g\|^{2}_{L^{2}(\Omega\times(0,T))}+\frac{\gamma}{2}\|u\|^{2}_{L^{2}(\Omega\times(0,T))},$$

(1)

subject to a parabolic equation with certain initial and boundary conditions

$$\left\{\begin{array}[]{lc}y_{t}-\mathcal{L}y=f+u,\quad(x,t)\in\Omega\times(0,T],\qquad y=0,\quad(x,t)\in\partial\Omega\times(0,T],\\ y(x,0)=y_{0},\quad x\in\Omega,\end{array}\right.\,$$

(2)

where $u,g\in L^{2}$ are the distributed control and the desired tracking trajectory, respectively, $\gamma>0$ is a regularization parameter, $\mathcal{L}=\nabla\cdot(a({x})\nabla)$, and $f$ and $y_{0}$ are given problem dependent functions. Under appropriate assumptions, theoretical aspects such as the existence, uniqueness, and regularity of the solution were well studied in [25]. The optimal solution of (1) & (2) can be characterized by the following system:

$$\left\{\begin{array}[]{lc}y_{t}-\mathcal{L}y-\frac{1}{\gamma}p=f,\quad(x,t)\in\Omega\times(0,T],\qquad y=0,\quad(x,t)\in\partial\Omega\times(0,T],\\ y(x,0)=y_{0}\quad x\in\Omega,\\ -p_{t}-\mathcal{L}p+y=g,\quad(x,t)\in\Omega\times(0,T],\qquad p=0,\quad(x,t)\in\partial\Omega\times(0,T],\\ p(x,T)=0\quad x\in\Omega,\end{array}\right.\,$$

(3)

where the control variable $u$ has been eliminated.

Following [36, 37], we discretize (3) using the $\theta$-method for time and some space discretization, which gives

	$\displaystyle M_{m}\frac{\mathbf{y}_{m}^{(k+1)}-\mathbf{y}_{m}^{(k)}}{\tau}+K_{m}(\theta\mathbf{y}_{m}^{(k+1)}+(1-\theta)\mathbf{y}_{m}^{(k)})$	$\displaystyle=$	$\displaystyle M_{m}(\theta\mathbf{f}_{m}^{(k+1)}+(1-\theta)\mathbf{f}_{m}^{(k)}+\frac{1}{\gamma}(\theta\mathbf{p}_{m}^{(k)}+(1-\theta)\mathbf{p}_{m}^{(k+1)})),$
	$\displaystyle-M_{m}\frac{\mathbf{p}_{m}^{(k+1)}-\mathbf{p}_{m}^{(k)}}{\tau}+K_{m}(\theta\mathbf{p}_{m}^{(k)}+(1-\theta)\mathbf{p}_{m}^{(k+1)})$	$\displaystyle=$	$\displaystyle M_{m}(\theta\mathbf{g}_{m}^{(k)}+(1-\theta)\mathbf{g}_{m}^{(k+1)}-\theta\mathbf{y}_{m}^{(k+1)}-(1-\theta)\mathbf{y}_{m}^{(k)}).$

The backward Euler method corresponds to $\theta=1$, while the Crank-Nicolson method is adopted when $\theta=1/2$.

Combining the given initial and boundary conditions, one needs to solve the following linear system

$$\widetilde{\mathcal{A}}\begin{bmatrix}\mathbf{y}\\ \mathbf{p}\end{bmatrix}=\begin{bmatrix}\mathbf{\widetilde{g}}\\ \mathbf{\widetilde{f}}\end{bmatrix},$$

(4)

where we have $\mathbf{y}=[\mathbf{y}_{m}^{(1)},\cdots,\mathbf{y}_{m}^{(n)}]^{\top}$, $\mathbf{p}=[\mathbf{p}_{m}^{(0)},\cdots,\mathbf{p}_{m}^{(n-1)}]^{\top}$,

$$\mathbf{\widetilde{f}}=\begin{bmatrix}M_{m}(\theta\tau\mathbf{f}_{m}^{(1)}+(1-\theta)\tau\mathbf{f}_{m}^{(0)})+(M_{m}-(1-\theta)\tau K_{m})\mathbf{y}_{m}^{(0)}\\ M_{m}(\theta\tau\mathbf{f}_{m}^{(2)}+(1-\theta)\tau\mathbf{f}_{m}^{(1)})\\ \vdots\\ M_{m}(\theta\tau\mathbf{f}_{m}^{(n)}+(1-\theta)\tau\mathbf{f}_{m}^{(n-1)})\\ \end{bmatrix},\mathbf{\widetilde{g}}=\tau\begin{bmatrix}M_{m}(\theta\mathbf{g}_{m}^{(0)}+(1-\theta)\mathbf{g}_{m}^{(1)}-(1-\theta)\mathbf{y}_{m}^{(0)})\\ M_{m}(\theta\mathbf{g}_{m}^{(1)}+(1-\theta)\mathbf{g}_{m}^{(2)})\\ \vdots\\ M_{m}(\theta\mathbf{g}_{m}^{(n-1)}+(1-\theta)\mathbf{g}_{m}^{(n)})\\ \end{bmatrix},$$

$\displaystyle\widetilde{\mathcal{A}}$

$\displaystyle=$

$\displaystyle\begin{bmatrix}\tau B_{n}^{(2)}\otimes M_{m}&(B_{n}^{(1)})^{\top}\otimes M_{m}+\tau(B_{n}^{(2)})^{\top}\otimes K_{m}\\ B_{n}^{(1)}\otimes M_{m}+\tau B_{n}^{(2)}\otimes K_{m}&-\frac{\tau}{\gamma}(B_{n}^{(2)})^{\top}\otimes M_{m}\end{bmatrix},$

(5)

and the matrices $B_{n}^{(1)},B_{n}^{(2)}$ are, respectively,

$$B_{n}^{(1)}=\begin{bmatrix}1&&&&\\ -1&1&&&\\ &-1&1&&\\ &&\ddots&\ddots&\\ &&&-1&1\end{bmatrix},\quad B_{n}^{(2)}=\begin{bmatrix}\theta&&&&\\ 1-\theta&\theta&&&\\ &1-\theta&\theta&&\\ &&\ddots&\ddots&\\ &&&1-\theta&\theta\\ \end{bmatrix}.$$

We assume that the matrix $M_{m}$ is symmetric positive definite, and the matrix $K_{m}$ is symmetric positive semi-definite. The matrices $M_{m}$ and $K_{m}$ represent the mass matrix and the stiffness matrix, respectively, if a finite element method is employed. For the finite difference method, the linear system is such that $M_{m}=I_{m}$ and $K_{m}=-L_{m}$, where $-L_{m}$ is the discretization matrix of the negative Laplacian.

Following a similar idea in [21], we can further transform (5) into the following equivalent system

$${\mathcal{A}}\begin{bmatrix}\sqrt{\gamma}\widetilde{\mathbf{y}}\\ \widetilde{\mathbf{p}}\end{bmatrix}=\begin{bmatrix}\mathbf{g}\\ \sqrt{\gamma}\mathbf{f}\end{bmatrix},$$

(6)

where $\mathbf{f}=(I_{n}\otimes M_{m}^{-\frac{1}{2}})\mathbf{\widetilde{f}}$, $\mathbf{g}=(I_{n}\otimes M_{m}^{-\frac{1}{2}})\mathbf{\widetilde{g}}$, $\widetilde{\mathbf{y}}=(B_{n}^{(2)}\otimes M_{m}^{\frac{1}{2}})[\mathbf{y}_{m}^{(1)},\cdots,\mathbf{y}_{m}^{(n)}]^{\top}$, $\mathbf{p}=((B_{n}^{(2)})^{\top}\otimes M_{m}^{\frac{1}{2}})[\mathbf{p}_{m}^{(0)},\cdots,\mathbf{p}_{m}^{(n-1)}]^{\top}$, and

	$\displaystyle{\mathcal{A}}$	$\displaystyle=$	$\displaystyle\begin{bmatrix}\alpha{I}_{n}\otimes I_{m}&B_{n}^{\top}\otimes I_{m}+\tau{I}_{n}\otimes M_{m}^{-\frac{1}{2}}K_{m}M_{m}^{-\frac{1}{2}}\\ B_{n}\otimes I_{m}+\tau{I}_{n}\otimes M_{m}^{-\frac{1}{2}}K_{m}M_{m}^{-\frac{1}{2}}&-\alpha{I}_{n}\otimes I_{m}\end{bmatrix}$
		$\displaystyle=$	$\displaystyle\begin{bmatrix}\alpha{I}_{n}\otimes I_{m}&\mathcal{T}^{\top}\\ \mathcal{T}&-\alpha{I}_{n}\otimes I_{m}\end{bmatrix}.$

Note that $\frac{1}{2}\leq\theta\leq 1$, $\alpha=\frac{\tau}{\sqrt{\gamma}}$, $I_{m}$ is the $m\times m$ identity matrix, and

$$\mathcal{T}=B_{n}\otimes I_{m}+\tau{I}_{n}\otimes M_{m}^{-\frac{1}{2}}K_{m}M_{m}^{-\frac{1}{2}},$$

(8)

where $B_{n}$ is a lower triangular Toeplitz matrix whose entries are known explicitly as

$$B_{n}=\begin{bmatrix}\frac{1}{\theta}&&&&\\ \frac{-1}{\theta^{2}}&\frac{1}{\theta}&&&\\ \frac{-(\theta-1)}{\theta^{3}}&\frac{-1}{\theta^{2}}&\frac{1}{\theta}&&\\ \vdots&\ddots&\ddots&\ddots&\\ \frac{-(\theta-1)^{n-2}}{\theta^{n}}&\cdots&\frac{-(\theta-1)}{\theta^{3}}&\frac{-1}{\theta^{2}}&\frac{1}{\theta}\end{bmatrix}.$$

Incidentally, $B_{n}$ can be expressed as the product of two Toeplitz matrices, i.e., $B_{n}=B_{n}^{(1)}(B_{n}^{(2)})^{-1}=(B_{n}^{(2)})^{-1}B_{n}^{(1)}$. As will be explained in Section 2, the Toeplitz matrices $B_{n}^{(1)}$ and $B_{n}^{(2)}$ are respectively generated by the functions

$$f_{1}(\phi)=1-\exp{(\mathbf{i}\phi})$$

(9)

and

$$f_{2}(\phi)=\theta+(1-\theta)\exp{(\mathbf{i}\phi)}.$$

(10)

In what follows, we focus on using the finite difference method to discretize the system (3), namely, $M_{m}=I_{m}$ and $K_{m}=-L_{m}$ in the linear system (6). However, we point out that our proposed preconditioning methods with minimal modification are still applicable when a finite element method is used. We first develop a preconditioned generalized minimal residual (GMRES) method for a nonsymmetric equivalent system of (6), which is

$$\mathcal{\widehat{A}}\begin{bmatrix}\sqrt{\gamma}\widetilde{\mathbf{y}}\\ \widetilde{\mathbf{p}}\end{bmatrix}=\begin{bmatrix}\sqrt{\gamma}\mathbf{f}\\ \mathbf{g}\end{bmatrix}$$

(11)

where

$\displaystyle\mathcal{\widehat{A}}$

$\displaystyle=$

$\displaystyle\begin{bmatrix}\mathcal{T}&-\alpha{I}_{n}\otimes I_{m}\\ \alpha{I}_{n}\otimes I_{m}&\mathcal{T}^{\top}\end{bmatrix}.$

(12)

For $\mathcal{\widehat{A}}$, we propose the following novel block preconditioner:

$\displaystyle\mathcal{P}_{S}=\begin{bmatrix}\mathcal{S}&-\alpha{I}_{n}\otimes I_{m}\\ \alpha{I}_{n}\otimes I_{m}&\mathcal{S}^{*}\end{bmatrix},$

(13)

where

$$\mathcal{S}=S_{n}\otimes I_{m}+\tau I_{n}\otimes(-L_{m}).$$

(14)

Notice that $S_{n}:=S_{n}^{(1)}(S_{n}^{(2)})^{-1}$, where

$\displaystyle S_{n}^{(1)}=\begin{bmatrix}1&&&&-\omega\\ -1&1&&&\\ &-1&1&&\\ &&\ddots&\ddots&\\ &&&-1&1\end{bmatrix},\quad S_{n}^{(2)}=\begin{bmatrix}\theta&&&&\omega(1-\theta)\\ 1-\theta&\theta&&&\\ &1-\theta&\theta&&\\ &&\ddots&\ddots&\\ &&&1-\theta&\theta\end{bmatrix}$

and $\omega=e^{\textbf{i}\zeta}\in\mathbb{C}$ with $\zeta\in[0,2\pi)$. Clearly, both $S_{n}^{(1)}$ and $S_{n}^{(2)}$ are $\omega$-circulant matrices [3, 2], so they admit the eigendecompositions

$$S_{n}^{(k)}=(\Gamma_{n}\mathbb{F}_{n})\Lambda_{n}^{(j)}(\Gamma_{n}\mathbb{F}_{n})^{*},\quad j=1,2,$$

(15)

where $\Gamma_{n}={\rm diag}(\exp{(\textbf{i}\zeta{(\frac{k-1}{n}}))})_{k=1}^{n}$ and $\mathbb{F}_{n}=\frac{1}{\sqrt{n}}[\theta_{n}^{(i-1)(j-1)}]_{i,j=1}^{n}$ with $\theta_{n}=\exp(\frac{2\pi{\bf i}}{n})$ and $\Lambda_{n}^{(j)}={\rm diag}(f_{1}(\frac{\zeta+2\pi k}{n}))_{k=0}^{n-1}$ and $f_{j}$ is defined by (9) and (10).

When $\omega=1$, $\mathcal{P}_{S}$ becomes the block circulant based preconditioner proposed in [36]. The existing block skew-circulant based preconditioner [4] proposed only with the backward Euler method is also included in our preconditioning strategy when $\omega=-1$ and $\theta=1$. As extensively studied in [8, 9], $\omega$-circulant matrices as preconditioners for Toeplitz systems can substantially outperform the Strang type preconditioners [32] (i.e., when $\omega=1$), especially in the ill-conditioned case. For related studies on the unsatisfactory performance of Strang preconditioners for ill-conditioned nonsymmetric (block) Toeplitz systems, we refer to [19, 18].

Since $\omega$-circulant matrices can be efficiently diagonalized by the fast Fourier transforms (FFTs), which can be parallelizable over different possessors. Hence, our preconditioner $\mathcal{P}_{S}$ is especially advantageous in a high performance computing environment.

In order to support our GMRES solver with $\mathcal{P}_{S}$ as a preconditioner, we will show that the singular values of $\mathcal{P}_{S}^{-1}\mathcal{\widehat{A}}$ are clustered around unity. However, despite of its success which can be seen in the numerical experiments from Section 4, the convergence study of preconditioning strategies for nonsymmetric problems is to a great extent heuristic. As mentioned in [35, Chapter 6], descriptive convergence bounds are usually not available for GMRES or any of the other applicable nonsymmetric Krylov subspace iterative methods.

Therefore, as an alternative solver, we develop a preconditioned minimal residual (MINRES) method for the symmetric system (6), instead of (11). Notice that our proposed MINRES method is in contrast with the aforementioned GMRES solvers, such as [36, 4] where a block (skew-)circulant type preconditioner was proposed and the eigenvalues of the preconditioned matrix were shown clustered around unity. As well explained in [13], the convergence behaviour of GMRES cannot be rigorously analyzed by using only eigenvalues in general. Thus, our MINRES solver can get round these theoretical difficulties of GMRES.

Based on the spectral distribution of $\mathcal{A}$, we first propose the following novel SPD block diagonal preconditioner as an ideal preconditioner for $\mathcal{A}$:

$\displaystyle|\mathcal{A}|:=\sqrt{\mathcal{A}^{2}}=\begin{bmatrix}\sqrt{\mathcal{T}^{\top}\mathcal{T}+\alpha^{2}{I}_{n}\otimes I_{m}}&\\ &\sqrt{\mathcal{T}\mathcal{T}^{\top}+\alpha^{2}{I}_{n}\otimes I_{m}}\end{bmatrix}.$

(16)

Despite its excellent preconditioning effect for $\mathcal{A}$, which will be shown in Section 3.2, the matrix $|\mathcal{A}|$ is computational expensive to invert. Thus, we then propose the following parallel-in-time (PinT) preconditioner, which mimics $|\mathcal{A}|$ and can be fast implemented:

$\displaystyle|\mathcal{P}_{S}|:=\sqrt{\mathcal{P}_{S}^{*}\mathcal{P}_{S}}=\begin{bmatrix}\sqrt{\mathcal{S}^{*}\mathcal{S}+\alpha^{2}{I}_{n}\otimes I_{m}}&\\ &\sqrt{\mathcal{S}\mathcal{S}^{*}+\alpha^{2}{I}_{n}\otimes I_{m}}\end{bmatrix}.$

(17)

However, the preconditioner $|\mathcal{P}_{S}|$ requires fast diagonalizability of $L_{m}$ in order to be efficiently implemented. When such diagonalizability is not available, we further propose the following preconditioner $\mathcal{P}_{MS}$ as a modification of $|\mathcal{P}_{S}|$:

			$\displaystyle\mathcal{P}_{MS}$
		$\displaystyle=$	$\displaystyle\begin{bmatrix}\sqrt{S_{n}^{*}S_{n}+\alpha^{2}I_{n}}\otimes I_{m}+\tau I_{n}\otimes(-L_{m})&\\ &\sqrt{S_{n}S_{n}^{*}+\alpha^{2}I_{n}}\otimes I_{m}+\tau I_{n}\otimes(-L_{m})\end{bmatrix}.$

One of our main contributions in this work is to develop a preconditioned MINRES method with the proposed preconditioners, which has theoretically guaranteed convergence based on eigenvalues.

It is worth noting that our preconditioning approaches are fundamentally different from another kind of related existing work (see, e.g., [31, 23]), which is a typical preconditioning approach in the context of preconditioning for saddle point systems. Its effectiveness is based on the approximation of Schur complements (e.g., [30, 21]) and the classical preconditioning techniques [1, 28]. Yet, for instance, our preconditioning proposal extends the MINRES preconditioning strategy proposed in [17] from optimal control of wave equations to that of parabolic equations, resulting in a clustered spectrum around $\{\pm 1\}$. Moreover, the implementation of our preconditioners based on FFTs are parallel-in-time.

The paper is organized as follows. In Section 2, we review some preliminary results on block Toeplitz matrices. In Section 3, we provide our main results on the spectral analysis for our proposed preconditioners. Numerical examples are given in Section 4 for supporting the performance of our proposed preconditioners.

In this section, we provide some useful background knowledge regarding Toeplitz matrices.

We let $L^{1}([-\pi,\pi])$ be the Banach space of all functions that are Lebesgue integrable over $[-\pi,\pi]$ and periodically extended to the whole real line. The Toeplitz matrix generated by $f\in L^{1}([-\pi,\pi])$ is denoted by $T_{n}[f]$, namely,

$\displaystyle T_{n}[f]=\begin{bmatrix}{}a_{0}&a_{-1}&\cdots&a_{-n+2}&a_{-n+1}\\ a_{1}&a_{0}&a_{-1}&&a_{-n+2}\\ \vdots&a_{1}&a_{0}&\ddots&\vdots\\ a_{n-2}&&\ddots&\ddots&a_{-1}\\ a_{n-1}&a_{n-2}&\cdots&a_{1}&a_{0}\end{bmatrix},$

where

$\displaystyle a_{k}=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(\theta)e^{-\mathbf{i}k\theta}\,d\theta,\quad k=0,\pm 1,\pm 2,\dots$

are the Fourier coefficients of $f$. The function $f$ is called the generating function of $T_{n}[f]$. If $f$ is complex-valued, then $T_{n}[f]$ is non-Hermitian for all sufficiently large $n$. Conversely, if $f$ is real-valued, then $T_{n}[f]$ is Hermitian for all $n$. If $f$ is real-valued and nonnegative, but not identically zero almost everywhere, then $T_{n}[f]$ is Hermitian positive definite for all $n$. If $f$ is real-valued and even, $T_{n}[f]$ is symmetric for all $n$. For thorough discussions on the related properties of block Toeplitz matrices, we refer readers to [11] and references therein; for computational features see [29, 7, 12] and references there reported.

In this section, the main results which support the effectiveness of our proposed preconditioners are provided. Also, the implementation issue is discussed.