A 2-Level Domain Decomposition Preconditioner for KKT Systems with Heat-Equation Constraints¹¹1To appear in the proceedings of The 27th International Conference on Domain Decomposition Methods (DD27)

This paper develops a new domain-decomposition method for solving the KKT system with heat-equation constraints. This effort is driven by the quadratic optimization problem of the form

	$\displaystyle\underset{z}{\min}$	$\displaystyle\frac{1}{2}\int_{0}^{T}\|u-\tilde{u}\|_{L^{2}(\Omega)}^{2}\,dt+\frac{\omega}{2}\int_{0}^{T}\|z\|_{L^{2}(\Omega)}^{2}\,dt$			(1)
	s.t.	$\displaystyle\partial_{t}u-\nu\nabla\cdot\nabla u=z,$	$\displaystyle x\in\Omega\subset\mathbb{R}^{2},t\in[0,T]$
		$\displaystyle u(x,t)=0,\;\;x\in{\partial\Omega},t\in[0,T],$	$\displaystyle u(x,0)=u_{0}(x),\;\;x\in\Omega$

This quadratic PDE-constrained optimization problem finds a control $z$ such that the solution $u$ to the heat equation matches the target $\tilde{u}$. The spatial domain is $\Omega$, the time interval is $[0,T]$, and the heat conductivity is $\nu$. Uniform homogenous boundary conditions are assumed for all time, and the initial condition is prescribed by $u_{0}$.

Many nonlinear methods use a series of quadratic approximations of the form represented by Eq. 1 to solve PDE-constrained optimization problems (see for instance sequential quadratic programming methods [8, 15, 27]). There have been several studies focused on developing scalable preconditioners for the saddle-point system that arises from the first-order necessary conditions. Often preconditioners for saddle-point systems take the form of approximate factorization block preconditioners [3]. These were explored for KKT systems in [4, 5]. Our work relies heavily on the block preconditioners from the Wathen group [23, 24, 25].

This effort focuses on transient PDE constraints where the size of the system scales with the number of spatial unknowns times the number of time steps, resulting in substantial computational effort. To alleviate this, a number of efforts have proposed accelerating the time solve using adaptive space-time discretizations [17, 16], parareal [20, 27, 12], multigrid approaches [6, 7, 10, 13, 19], block preconditioning [23, 24], and domain decomposition methods [14].

Our approach is also built on block preconditioning ideas. A difference is that our technique exploits an observation that the Schur-complement of the KKT system is elliptic in time (see [11, 18]). This allows us to leverage existing two level domain decomposition approaches for elliptic systems to improve the parallel scalability of the block preconditioner. Good performance is achieved by algorithmic choices that ensure the forward and backward in time integrators can be applied on the fine level.

In this article the PDE in Eq. 1 will be discretized on a $2$D Cartesian grid using first order backward Euler in time, and a second order finite difference stencil in space. A row of the discrete space-time system for the heat equation satisfies:

$$u_{i,j}^{n+1}-u_{i,j}^{n}+\Delta t\nu\left(\frac{-u_{(i+1)j}^{n+1}+2u_{ij}^{n+1}-u_{(i-1)j}^{n+1}}{\Delta x^{2}}+\frac{-u_{i(j+1)}^{n+1}+2u_{ij}^{n+1}-u_{i(j-1)}^{n+1}}{\Delta y^{2}}\right)=\Delta tz_{ij}^{n+1}.$$

(2)

Here $i,j$ are the interior space indices defined over $1\dots n_{x}-1$ and $1\dots n_{y}-1$. The exterior indices are eliminated using the homogenous boundary conditions. The superscript time index $n$ runs from $0\ldots N_{t}$. Each $n$ is referred to below as a time-node. The control variable $z$’s index matches the implicit index on $u$, therefore $z^{n+1}$ is associated with the $n$th time interval. For a single time interval, Eq. 2 rewritten in matrix form is

$$J_{(n+1)(n+1)}u^{n+1}+J_{(n+1)n}u_{n}+L_{(n+1)(n+1)}z^{n+1}=0,$$

(3)

and the global space-time system is

$$Ju+Lz=f.$$

(4)

The right hand side $f$ includes contributions from the initial conditions. The matrix $J$ is block lower triangular and the matrix $L$ is block diagonal.

The linear system whose solution solves the quadratic optimization problem from Eq. 1 is the celebrated KKT system $K\vec{u}=\vec{f}$ where

$$K=\begin{bmatrix}M_{u}&&J^{T}\\ &\omega M_{z}&L^{T}\\ J&L&\end{bmatrix},\;\vec{u}=\begin{bmatrix}u\\ z\\ w\end{bmatrix},\;\vec{f}=\begin{bmatrix}f_{u}\\ f_{z}\\ f\end{bmatrix}.$$

(5)

The final row is the discrete form of the PDE constraint, enforced by the Lagrange multiplier $w$. We will also refer to $w$ as the adjoint solution. $M_{u}$ and $M_{z}$ are identity matrices scaled by $\Delta t\Delta x\Delta y$. The matrix $K$ is a saddle point matrix, whose structure is frequently observed in numerical optimization. Many effective block preconditioners have been developed for this class of matrix [2, 3, 4]. We focus on the block preconditioning approach developed by Wathen and collaborators for solving linearized PDE-constrained optimization problems [23, 24, 25].

We write a block LDU factorization of the matrix $K$

$$K=\begin{bmatrix}I&&\\ &I&\\ JM_{u}^{-1}&\omega^{-1}LM_{z}^{-1}&I\end{bmatrix}\begin{bmatrix}M_{u}&&\\ &\omega M_{z}&\\ &&-S\end{bmatrix}\begin{bmatrix}I&&M_{u}^{-1}J^{T}\\ &I&\omega^{-1}M_{z}^{-1}L^{T}\\ &&I\end{bmatrix}$$

(6)

where the Schur-complement is $S=JM_{u}^{-1}J^{T}+\frac{1}{\omega}LM_{z}^{-1}L^{T}$. Following [23], $K$ is preconditioned using the block diagonal operator

$$P=\begin{bmatrix}M_{u}&&\\ &\omega M_{z}&\\ &&\hat{S}\end{bmatrix}\mbox{ where }\hat{S}=\hat{J}M_{u}^{-1}\hat{J}^{T},\hat{J}=J+\omega^{-1/2}L.$$

(7)

This preconditioner leverages the result in [21], and approximately inverts the block diagonal in the LDU factorization. The matrix $\hat{J}$ used in the approximate Schur complement $\hat{S}$ is block lower triangular (similar to $J$), a fact that we will exploit below. The choice of $\hat{S}$ integrates the state Jacobian and the effects of the regularization parameter. In [23, 24] and [26], this approximation is developed and shown to lead to robust performance with respect to $\omega$.

We propose a new domain decomposition approach for approximately inverting $\hat{S}$. This is motivated by the observation that the operator $S$ is elliptic in time (see [18] and [11]). For simplicity, we show this discretely using only the term $JM_{u}^{-1}J^{T}$. Consider the ODE $\partial_{t}y=-y$ discretized over three time steps with forward Euler: $y^{n+1}-y^{n}+\Delta ty^{n}=0$. With $M_{u}=I$, the Schur-complement $\hat{S}$ is

$$\begin{bmatrix}1\\ -1+\Delta t&1&&\\ &-1+\Delta t&1&\\ \end{bmatrix}\begin{bmatrix}1&-1+\Delta t&\\ &1&-1+\Delta t\\ &&1\\ \end{bmatrix}\\ =\begin{bmatrix}1&-(1-\Delta t)&\\ -(1-\Delta t)&2(1-\Delta t)+\Delta t^{2}&-(1-\Delta t)\\ &-(1-\Delta t)&2(1-\Delta t)+\Delta t^{2}\end{bmatrix}.$$

(8)

Examining the second row it is clear the operator has a $1$D Laplacian stencil in time, with a positive perturbation on the diagonal. To take advantage of this ellipticity, we will apply existing domain decomposition approaches to the $\hat{S}$ operator. This ellipticity principle enables scalable performance of a preconditioned Krylov method. We also impose an efficiency constraint that the computational kernels in our preconditioner use the time integration method for the state and adjoint unknowns.

The ellipticity principle is realized by considering a restricted additive Schwarz (RAS) method with $N_{D}$ subdomains (see [9]). Each subdomain contains the spatial unknowns associated with a subset of time steps. For instance, if there are $N_{t}=9$ time steps, then for $N_{D}=3$ the subdomains contain time-nodes $\{1,2,3\}$, $\{3,4,5,6\}$, and $\{6,7,8,9\}$ (the $0^{\text{th}}$ time-node is the excluded initial condition). Notice that the time-nodes are overlapped but the time steps are not. With these subsets, boolean operators $R_{s}$ are defined that restrict a space-time vector to the time-nodes in a subdomain, giving the RAS preconditioner

$$\hat{S}_{\text{RAS}}^{-1}=\sum_{s=1}^{N_{D}}R_{s}^{T}D_{s}\left(R_{s}\hat{J}M_{u}^{-1}\hat{J}^{T}R_{s}^{T}\right)^{-1}R_{s}$$

(9)

where $D_{s}$ is the (boolean) partition of unity matrix (Defn. 1.11 of [9]). RAS is known to lead to effective preconditioners for elliptic problems and can be extended to a multi-level schemes. While the ellipticity principle is exploited in $\hat{S}_{\text{RAS}}^{-1}$, the explicit formation of the product $\hat{J}M_{u}^{-1}\hat{J}^{T}$ does not satisfy the efficiency constraint.

To satisfy the efficiency constraint note that the range of $\hat{J}^{T}R_{s}^{T}$ is nonzero on time-nodes in the $s^{\text{th}}$ subdomain and one time-node earlier. For instance, if the subdomain contains nodes $\{3,4,5,6\}$ then the range is nonzero on $\{2,3,4,5,6\}$. Let $Q_{s}$ be a new extended restriction operator whose action produces a space-time vector for time-nodes that are nonzero in the range of $\hat{J}^{T}R_{s}^{T}$. Further, choose $Q_{s}$ so that

$$Q_{s}=\begin{bmatrix}W_{s}\\ R_{s}\end{bmatrix}\mbox{ and }Q_{s}Q_{s}^{T}=I.$$

(10)

The operator $W_{s}$ restricts the space-time vector to the time-nodes contained in the earlier time step relative to the $s^{\text{th}}$ subdomain. Because $Q_{s}$ is a restriction operator that represents the nonzero range of $\hat{J}^{T}R_{s}^{T}$, we have that $Q_{s}^{T}Q_{s}\hat{J}^{T}R_{s}^{T}=\hat{J}^{T}R_{s}^{T}$. Recalling that $M_{u}$ is diagonal, we can rewrite the subdomain solve in $\hat{S}_{RAS}^{-1}$ as

$$R_{s}\hat{J}M_{u}^{-1}\hat{J}^{T}R_{s}^{T}=R_{s}\hat{J}Q_{s}^{T}(Q_{s}M_{u}^{-1}Q_{s}^{T})Q_{s}\hat{J}^{T}R_{s}^{T}$$

(11)

Using the constraint in Eq. 10, we have the additional identities $R_{s}=R_{s}Q_{s}^{T}Q_{s}$ and $R_{s}^{T}=Q_{s}^{T}Q_{s}R_{s}^{T}$. This permits a final rewrite of the operator in Eq. 11

$$R_{s}\hat{J}M_{u}^{-1}\hat{J}^{T}R_{s}^{T}=R_{s}Q_{s}^{T}\hat{J}_{s}M_{s}^{-1}\hat{J}_{s}^{T}Q_{s}R_{s}^{T}$$

(12)

where $\hat{J}_{s}=Q_{s}\hat{J}Q_{s}^{T}$ and $M_{s}^{-1}=Q_{s}M_{u}^{-1}Q_{s}^{T}$. The inverse action of $\hat{J}_{s}M_{s}^{-1}\hat{J}_{s}^{T}$ is easily computed in a matrix free way on the time-nodes in the extended subdomain. Motivated by this equivalence, we define a new one-level preconditioner

$$\hat{S}_{\text{RASQ}}^{-1}=\sum_{s=1}^{N_{D}}R_{s}^{T}D_{s}\left(R_{s}Q_{s}^{T}\hat{J}_{s}^{-T}M_{s}\hat{J}_{s}^{-1}Q_{s}R_{s}^{T}\right)R_{s}.$$

(13)

The term in parentheses is different from the term inverted in Eq. 9. The difference is that the inverse computed in Eq. 9 is constrained to have a zero initial condition outside of the subdomain. This revised operator satisfies our efficiency constraint as computing $\hat{J}_{s}^{-1}$ and $\hat{J}_{s}^{-T}$ is done using the time integration method.

To obtain scalability with respect to the number of subdomains a coarse grid correction is required. Again leveraging the elliptic nature of the Schur-complement we consider the Nicolaides coarse space developed for solving the Poisson problem [22]. Following the presentation in [9], the Nicolades coarse space is defined as

$$Z=\begin{bmatrix}D_{1}R_{1}\Phi_{1}&&\ldots&\\ &D_{2}R_{2}\Phi_{2}&\ddots&\\ \vdots&\ddots&\ddots&\\ &\ldots&&D_{N_{D}}R_{N_{D}}\Phi_{N_{D}}\end{bmatrix}\mbox{ where }\Phi_{s}=\begin{bmatrix}1&-1\\ 1&-1+\frac{2}{N_{s}}\\ \vdots&\vdots\\ 1&-1+2\frac{N_{s}}{N_{s}}\end{bmatrix}$$

(14)

The columns of $\Phi_{s}$ form a constant and linear basis over the $N_{s}$ subdomain time-nodes. The coarse restriction $R_{0}=Z^{T}$ is used in the definition of the coarse operator $\hat{S}_{0}=R_{0}\hat{J}M_{u}^{-1}\hat{J}^{T}R_{0}^{T}$. The coarse solve is applied in a multiplicative way

$$\hat{S}_{2-\text{level}}^{-1}=\hat{S}_{\text{RASQ}}^{-1}R_{0}^{T}\hat{S}_{0}^{-1}R_{0}.$$

(15)

Due to the structure of $R_{0}$ the coarse operator $\hat{S}_{0}$ can be constructed in parallel. This does represent a violation of the efficiency constraint to be addressed by future work.

To demonstrate this approach we discretize the quadratic optimization problem from Eq. 1 as described in Sec. 2. The $2$D spatial domain is $\Omega=(0,1)\times(0,2)\subset\mathbb{R}^{2}$, and the time domain is $[0,1]$. The initial conditions and target solutions are

$$u_{0}(x,y)=-xy(x-1)(y-2),\quad\tilde{u}(x,y,t)=\sin(2.0\pi t)\sin(2.0\pi x)\sin(2.0\pi y).$$

(16)

The regularization parameter $\omega$ varies over five orders of magnitude. Experiments were run with $9\times 9$, $17\times 17$, and $33\times 33$ mesh points. Qualitatively, variation with number of spatial points was not a factor in the convergence. This is not surprising as the implicit operator in space is inverted with a direct solve. As a result, the computations below are all for the case of $17\times 17$ mesh points. Recall that homogeneous boundary conditions are removed, giving $15\times 15$ unknowns in each time step. The linear system $K\vec{u}=\vec{f}$ (Eq. 5) is solved using right preconditioned GMRES from PyAMG [1] iterated until a relative residual tolerance of $10^{-6}$ is achieved.

Figure 1 presents three weak-scaling studies ranging from $100$ to $3200$ time steps. The number of time steps per subdomain is fixed at $80$ in the left plot, $20$ in the center plot, and $5$ in the right plot. For the case of $80$ steps, the fewest number of time steps is $200$ (the minimum number of time steps evenly divisible by $80$ in the chosen sequence). The plots, show the number of iterations as a function of time step count for GMRES preconditioned using $P$ from Eqn 7 with Schur complement approximations $\hat{S}_{\text{RASQ}}^{-1}$ for the one level case (circle markers), and $\hat{S}_{2-\text{level}}^{-1}$ for the two level case (triangle markers). Different values for the regularization parameter $\omega$ are indicated using solid ($10^{-2}$), dashed ($10^{-3}$) or dotted ($10^{-4}$) lines. These plots demonstrate that the performance of the two level method is independent of both the number of subdomains, and the number of time steps. Further, independence holds regardless of the value of the regularization parameter. As anticipated the one level method has substantial growth with the number of time steps, and variability with the regularization parameter. However, it is worth noting that dependent on the number of subdomains and the size of the regularization parameter the one level method may be faster despite its lack of scalability. For instance, when using $40$ subdomains and a regularization parameter of $10^{-4}$ the one-level method takes the same number of iterations but lacks the synchronization and added cost of the two level method.

The scaling with respect to the regularization parameter is investigated in Figure 2. In these plots the preconditioned iteration counts are plotted as a function of the inverse regularization parameter. Data points are excluded when the number of iterations exceeded the maximum iteration count for GMRES (in this case $420$). Here again the two level method scales well, yielding essentially flat iteration counts as a function of the regularization parameter. The one level method shows strong dependence on $\omega$, though it improves dramatically for smaller values.

In this paper, motivated by results in block preconditioning and the elliptic-in-time nature of the KKT system, we develop a two level domain decomposition preconditioner that facilitates a parallel-in-time solver for the discrete optimality system constrained by the heat equation. While limited in their breadth, initial results for this approach show excellent scalability with respect to the number of time steps, subdomains, and the regularization parameter. Future work will focus on achieving improved scaling by including more levels in the hierarchy, and applying this technique to a broader class of problems and discretizations.

The author is indebted anonymous reviewers whose comments resulted in a significant strengthening of the paper. The author also acknowledges support from the SEA-CROGS project and the Early Career program, both funded by the DOE Office of Science. This article has been authored by an employee of National Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee owns all right, title and interest in and to the article and is solely responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan https://www.energy.gov/downloads/doe-public-access-plan. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. Sandia report number: SAND2023-03175O.