An optimal parallel-in-time preconditioner for parabolic optimal control problems

Developing diagonalization-based parallel-in-time (PinT) for solving optimization problems constrained by partial differential equations (PDEs) has gained substantial attention in recently years. We refer to [19, 12, 27, 9] and the references therein for a comprehensive overview of these optimization problems. Various efficient PinT preconditioners have been proposed to solve all-at-once linear systems arising from time-dependent PDEs [23, 22, 17, 20, 15] and PDE-constrained optimal control problems [14, 29, 30, 31, 13, 8, 10]; see for a review paper on these diagonalization-based methods [21] and the references therein reported.

In this work, we consider solving the distributed optimal control problem constrained by a parabolic equation. Namely, the following quadratic cost functional is minimized:

constrained by a parabolic equation with certain initial and boundary conditions

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

where the control variable $u$ has been eliminated.

Similar to [31, 30], we discretize (3) using the backward Euler method for time and a typical finite element method for space discretization, which gives

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

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}$,

and the matrix $B_{n}\in\mathbb{R}^{n\times n}$ is as follows.

The matrices $M_{m}$ and $K_{m}$ represent the mass matrix and the stiffness matrix, respectively. For the finite difference method, we have $M_{m}=I_{m}$ and $K_{m}=-L_{m}$, where $-L_{m}$ is the discretization matrix of the negative Laplacian. Throughout, the following conditions on both $M_{m}$ and $K_{m}$ are assumed, which are in line with [17, Section 2]:

1. 1.

   The matrices $M_{m}$ and $K_{m}$ are (real) symmetric positive definite;

2. 2.

   The matrices $M_{m}$ and $K_{m}$ are sparse, having only $\mathcal{O}(m)$ nonzero elements;

3. 3.

   The matrix $M_{m}$ has a condition number uniformly upper bounded with respect to $m$, i.e., there exists a positive finite constant $c_{0}$ independent of $m$ and $n$ such that

   (7)

   $$\sup\limits_{m>0}\kappa_{2}(M_{m})\leq c_{0}<+\infty,$$

   where $\kappa_{2}(M_{m}):=||M_{m}||_{2}||M_{m}^{-1}||_{2}$ denotes the condition number of $M_{m}$.

In fact, these assumptions are easily met by applying typical finite element discretization to the spatial term $\mathcal{L}$.

We then further transform (5) into the following equivalent system

where

Note that $\alpha=\frac{\tau}{\sqrt{\gamma}}$, $I_{k}$ denotes the $k\times k$ identity matrix, and

In what follows, we will develop a preconditioned generalized minimal residual (GMRES) method for the nonsymmetric system (8).

For $\mathcal{A}$, we first consider the following rotated block diagonal (RBD) preconditioner proposed in [4]:

where $\mathcal{T}$ is defined in (10). As will be shown in Theorem 2.4, the matrix $\mathcal{P}$ can achieve an excellent preconditioning effect for $\mathcal{A}$. Yet, $\mathcal{P}$ is in general computationally expensive to invert. Thus, we propose the following PinT preconditioner (or denoted by RBD $\epsilon$-circulant preconditioner) as a modification of $\mathcal{P}$, which can be efficiently implemented:

where

Note that $C_{\epsilon,n}$ is defined as

and $\epsilon\in(0,1]$. Since $C_{\epsilon,n}$ is an $\epsilon$-circulant matrix, it admits the following decomposition [7, 6]

where

and

Several efficient preconditioned iterative solvers have been proposed for (5). One notable class of such methods is the circulant-based preconditioners [31, 8, 10], whose main idea is to replace the Toeplitz matrix $B_{n}$ defined by (6) by a circulant or skew-circulant matrix in order to form an efficient preconditioner. Despite their numerically observed success, theoretically these preconditioners have not shown to be parameter-robust due to the low-rank approximation from the circulant-based approximation.

Another major existing method is preconditioning technique based on Schur complement approximations [30, 14, 18], which is a typical preconditioning approach in the context of preconditioning for saddle point systems (see, e.g., [26, 28]) and their effectiveness is based on approximating the Schur complements by multiplications of easily invertible matrices (e.g., [25, 14, 2, 24]). Recently, $\epsilon$-circulant modifications have been introduced to such approximate-Schur-complement preconditioning techniques for the optimal control problems, which leads to a significant improvement on computational efficiency and parallelism along time dimension (see, e.g, [17, 20, 15, 16, 30, 21]).

Our proposed preconditioning method, while still based on $\epsilon$-circulant matrices, diverges from the approximate-Schur-complement approach. Specifically, the Toeplitz matrix $B_{n}$ is formulated using an $\epsilon$-circulant matrix as detailed in (13). An important advancement of our method is that it can achieve optimal convergence rate when $\epsilon=\mathcal{O}(\tau^{1/2})$. This contrasts with $\epsilon$-circulant modified-approximate-Schur-complement preconditioner, such as described in [18], which requires a more restrictive $\epsilon=\mathcal{O}(\tau^{2})$. Additionally, while the other $\epsilon$-circulant modified-approximate-Schur-complement approach in [30] also allows for $\epsilon=\mathcal{O}(\tau^{1/2})$, it imposes the severe theoretical constraint that $n$ must be an even number, a limitation our method does not share.

The paper is organized as follows. In Section 2, we provide our main results on the spectral analysis for our proposed preconditioners. Numerical examples are given in Section 3 for supporting the performance of our proposed preconditioner. As last, conclusions are given in Section 4.

Before presenting our main preconditioning result, we introduce some useful notation in what follows.

Let

and

Then, it is straightforward to verify that

To check the invertibility of $\mathcal{P}_{\epsilon}$, it suffices to verify the invertibility of $\mathcal{G}$ and $\mathcal{H}_{\epsilon}$ respectively. Clearly, $\mathcal{G}$ is invertible. By (14), it is straightforward to see that $\mathcal{C}_{\epsilon}+\alpha I_{n}\otimes M_{m}$ is similar to a block diagonal matrix with each diagonal block having the form $\alpha M_{m}+z_{1}I_{m}+z_{2}{\bf i}I_{m}$ ($z_{1}\geq 0$, $z_{2}\in\mathbb{R}$). Consequently, $\mathcal{C}_{\epsilon}+\alpha I_{n}\otimes M_{m}$ is invertible. Moreover, $\mathcal{C}_{\epsilon}^{\top}+\alpha{I}_{n}\otimes M_{m}$ as the transpose of $\mathcal{C}_{\epsilon}+\alpha{I}_{n}\otimes M_{m}$ is also invertible. Thus, $\mathcal{H}_{\epsilon}$, as a block diagonal matrix with $\mathcal{C}_{\epsilon}^{\top}+\alpha{I}_{n}\otimes M_{m}$ and $\mathcal{C}_{\epsilon}+\alpha{I}_{n}\otimes M_{m}$ as its diagonal blocks, is also invertible. Therefore, $\mathcal{P}_{\epsilon}$ is invertible.

In this work, we propose the use of the matrix $\mathcal{P}_{\epsilon}$ to precondition the saddle point system described in (8). In other words, we shall apply the GMRES solver to solve the following preconditioned system

We will focus on investigating the convergence rate of the GMRES solver for the preconditioned system mentioned above, by theoretically estimating a suitable choice of $\epsilon$.

Denote

It is then easy to see that

With the equations above, it is straightforward to verify that (15) can be rewritten as

meaning that (15) can be equivalently converted into the following system

where $\tilde{\bf b}=\widetilde{\mathcal{H}}_{\epsilon}^{-1}\mathcal{W}^{-\frac{1}{2}}{\bf b}$; the unknown $\tilde{\bf x}$ in (16) and the unknown ${\bf x}$ in (15) are related by

We will examine the theoretical convergence rate of the GMRES solver for the linear system (16) and subsequently relate it to the convergence rate of the GMRES solver for (15) using Theorem 2.2.

The convergence behavior of GMRES solver is closely related to the Krylov subspace. For a square matrix $\mathbf{E}\in\mathbb{R}^{k\times k}$ and a vector ${\bf y}\in\mathbb{R}^{k\times 1}$, a Krylov subspace of degree $j\geq 1$ is defined as follows

For a set $\mathcal{S}$ and a point $z$, we denote

We recall the relation between the iterative solution by GMRES and the Krylov subspace in the following lemma.

According to Theorem 2.2, any convergence rate estimation for the GMRES solver applied to (16) is also applicable to (15). Namely, the GMRES convergence for (15) is guaranteed by its convergence for (16). As a result, in the following discussion, we will focus on investigating the convergence of the GMRES solver for (16) by first examining the spectrum of its coefficient matrix, $\widetilde{\mathcal{H}}_{\epsilon}^{-1}\widetilde{\mathcal{B}}$.

Denote

Clearly, $\widetilde{\mathcal{H}}_{\epsilon}$ is an approximation to $\widetilde{\mathcal{H}}$. In the next, we will study the spectrum of $\widetilde{\mathcal{H}}^{-1}\widetilde{\mathcal{B}}$ and then relate the spectrum of $\widetilde{\mathcal{H}}_{\epsilon}^{-1}\widetilde{\mathcal{B}}$ to the spectrum of $\widetilde{\mathcal{H}}^{-1}\widetilde{\mathcal{B}}$ by the fact $\lim\limits_{\epsilon\rightarrow 0^{+}}\widetilde{\mathcal{H}}_{\epsilon}=\widetilde{\mathcal{H}}$.

The following theorem shows that $\widetilde{\mathcal{H}}$ is an ideal preconditioner for $\widetilde{\mathcal{B}}$. We remark that the proof mirrors that of [4, Theorem 4.1], with modifications. However, for self-containment, we decide to provide a complete proof below.

Denote by $\sigma({\bf C})$, the spectrum of a square matrix ${\bf C}$.

It is straightforward to verify that $\widetilde{\mathcal{T}}+\alpha I_{mn}$ can be rewritten as

and that $(\widetilde{\mathcal{T}}+\alpha I_{mn})^{-1}$ has the following expression

Denote