Multiscale differential Riccati equations for linear quadratic regulator problems

TL;DR

This paper introduces a multiscale approach using localized orthogonal decomposition to efficiently approximate solutions to differential Riccati equations in linear quadratic regulator problems with multiscale operators, achieving high accuracy with low computational cost.

Contribution

It develops a novel multiscale method for Riccati equations that maintains accuracy independent of multiscale variations and provides detailed analysis and validation.

Findings

Second-order convergence in $L^2$ operator norm

First-order convergence in energy norm

Effective low-rank computations for large-scale problems

Abstract

We consider approximations to the solutions of differential Riccati equations in the context of linear quadratic regulator problems, where the state equation is governed by a multiscale operator. Similarly to elliptic and parabolic problems, standard finite element discretizations perform poorly in this setting unless the grid resolves the fine-scale features of the problem. This results in unfeasible amounts of computation and high memory requirements. In this paper, we demonstrate how the localized orthogonal decomposition method may be used to acquire accurate results also for coarse discretizations, at the low cost of solving a series of small, localized elliptic problems. We prove second-order convergence (except for a logarithmic factor) in the $L^2$ operator norm, and first-order convergence in the corresponding energy norm. These results are both independent of the multiscale variations in the state equation. In addition, we provide a detailed derivation of the fully discrete matrix-valued equations, and show how they can be handled in a low-rank setting for large-scale computations. In connection to this, we also show how to efficiently compute the relevant operator-norm errors. Finally, our theoretical results are validated by several numerical experiments.

Equations

$$\dot{x}$$

$$y$$

$$J(u)={\int }_0^T\left(\mathcal{Q}y,y\right)+\left(\mathcal{R}u,u\right)\mathrm{d}t+\left(\mathcal{G}y(T),y(T)\right),$$

$$\begin{array}{rl}\dot{X}(t)& ={\mathcal{A}}^{\ast }X(t)+X(t)\mathcal{A}+{\mathcal{C}}^{\ast }\mathcal{Q}\mathcal{C}-X(t)\mathcal{B}{\mathcal{R}}^{-1}{\mathcal{B}}^{\ast }X(t),\\ X(0)& =\mathcal{G}.\end{array}$$

$$\left(\dot{X}x,y\right)=\left(Xx,\mathcal{A}y\right)+\left(Xy,\mathcal{A}x\right)+{\left(\mathcal{Q}\mathcal{C}x,\mathcal{C}y\right)}_Z-{\left({\mathcal{R}}^{-1}{\mathcal{B}}^{\ast }Xx,{\mathcal{B}}^{\ast }Xy\right)}_U,$$

$$0$$

$$\infty$$

$$\gamma :=\max {\gamma }_K,{\gamma }_K:=\frac{\mathrm{diam}{B}_K}{\mathrm{diam}K},\forall K\in {\mathcal{T}}_H,H>h.$$

$$\left({\dot{X}}_{h}x,y\right)=\left(X_{h}x,{\mathcal{A}}_{h}y\right)+\left(X_{h}y,{\mathcal{A}}_{h}x\right)+{\left(\mathcal{Q}{\mathcal{C}}_{h}x,{\mathcal{C}}_{h}y\right)}_Z-{\left({\mathcal{R}}^{-1}{\mathcal{B}}_h^{\ast }{X}_{h}x,{\mathcal{B}}_h^{\ast }{X}_{h}y\right)}_U.$$

$$\left({\mathcal{A}}_{h}x,y\right)=\left(\mathcal{A}x,y\right),\left({\mathcal{B}}_{h}u,y\right)=\left(\mathcal{B}u,y\right)\text{and}\left({\mathcal{C}}_{h}x,z\right)=\left(\mathcal{C}x,z\right)$$

$${\mathcal{A}}_H=({\mathrm{Id}}_H^h{)}^{\ast }{\mathcal{A}}_h{\mathrm{Id}}_H^h,{\mathcal{B}}_H=({\mathrm{Id}}_H^h{)}^{\ast }{\mathcal{B}}_h\text{and}{\mathcal{C}}_H={\mathcal{C}}_h{\mathrm{Id}}_H^h,$$

$$H^{-1}\Vert v-I_{H}v{\Vert }_{L^2(K)}+\Vert \nabla I_{H}v{\Vert }_{L^2(K)}\le C\Vert \nabla v{\Vert }_{L^2({\omega }_K)},\forall v\in V_h,$$

$$V_{\text{f}}:=\ker I_H=\{v\in V_h:I_{H}v=0\},$$

$$a({\hat{Q}}_{h}v,w)=a(v,w),\forall w\in V_{\text{f}},$$

$${\omega }_0(K)$$

$${\omega }_k(K)$$

$$a({\hat{Q}}_h^{K}v,w)={\int }_K\kappa \nabla v\cdot \nabla w,\forall w\in V_{\text{f}},v\in V_H,K\in {\mathcal{T}}_H.$$

$$a(Q_h^{K}v,w)={\int }_K\kappa \nabla v\cdot \nabla w,\forall w\in V_{\text{f}}^K,v\in V_H,K\in {\mathcal{T}}_H.$$

$$\left({\dot{X}}_h^{\text{ms}}u,v\right)$$

$$+{\left(\mathcal{Q}{\mathcal{C}}_h^{\text{ms}}u,{\mathcal{C}}_h^{\text{ms}}v\right)}_Z-{\left({\mathcal{R}}^{-1}({\mathcal{B}}_h^{\text{ms}}{)}^{\ast }{X}_h^{\text{ms}}u,({\mathcal{B}}_h^{\text{ms}}{)}^{\ast }{X}_h^{\text{ms}}v\right)}_U$$

$$\left({\mathcal{A}}_h^{\text{ms}}v,w\right)=\left(\mathcal{A}v,w\right),{\left({\mathcal{B}}_h^{\text{ms}}u,w\right)}_U={\left(\mathcal{B}u,w\right)}_U\text{and}{\left({\mathcal{C}}_h^{\text{ms}}v,z\right)}_Z={\left(\mathcal{C}u,z\right)}_Z$$

$${\mathcal{A}}_h^{\text{ms}}=({\mathrm{Id}}_{\text{ms}}^h{)}^{\ast }{\mathcal{A}}_h{\mathrm{Id}}_{\text{ms}}^h,{\mathcal{B}}_h^{\text{ms}}=({\mathrm{Id}}_{\text{ms}}^h{)}^{\ast }{\mathcal{B}}_h\text{and}{\mathcal{C}}_h^{\text{ms}}={\mathcal{C}}_h{\mathrm{Id}}_{\text{ms}}^h.$$

$$\left({\dot{X}}_H^{\text{ms}}x,R_h^{\ast }{R}_{h}y\right)$$

$$+{\left(\mathcal{Q}{\mathcal{C}}_h{R}_{h}x,{\mathcal{C}}_h{R}_{h}y\right)}_Z-{\left({\mathcal{R}}^{-1}{\mathcal{B}}_h^{\ast }{R}_h{X}_H^{\text{ms}}x,{\mathcal{B}}_h^{\ast }{R}_h{X}_H^{\text{ms}}y\right)}_U,$$

$${\mathcal{A}}_H^{\text{ms}}=R_h^{\ast }{\mathcal{A}}_h{R}_h,{\mathcal{B}}_H^{\text{ms}}=R_h^{\ast }{\mathcal{B}}_h\text{and}{\mathcal{C}}_H^{\text{ms}}={\mathcal{C}}_h{R}_h.$$

$$E_h(t)={\mathrm{Id}}_h{\mathrm{e}}^{t{\mathcal{A}}_h}{\mathrm{Id}}_h^{\ast }\text{and}{E}_{\text{ms}}(t)={\mathrm{Id}}_h{\mathrm{Id}}_{\text{ms}}^h{\mathrm{e}}^{t{\mathcal{A}}_h^{\text{ms}}}({\mathrm{Id}}_{\text{ms}}^h{)}^{\ast }{\mathrm{Id}}_h^{\ast },$$

$$\Vert E_h(t)-E_{\text{ms}}(t){\Vert }_{\mathcal{L}(L^2)}\le C{H}^2{t}^{-1}.$$

$$\displaystyle\lVert u_{h}-u_{\text{ms}}\rVert\leq Ct^{-1}\sup_{s\leq t}\Big{\{}s^{2}\lVert\dot{\rho}\rVert+s\lVert\rho\rVert+\lVert\int_{0}^{s}\rho(r)\,\mathrm{d}r\rVert\Big{\}},$$

$$\lVert u_{h}-u_{\text{ms}}\rVert\leq CH^{2}t^{-1}\sup_{s\leq t}\big{(}s^{2}\lVert\ddot{u}_{h}(s)\rVert+s\lVert\dot{u}_{h}(s)\rVert+\lVert u_{h}(s)\rVert+\lVert v\rVert\big{)}.$$

$$\tilde{X}(t)={\mathrm{Id}}_h{X}_h(t){\mathrm{Id}}_h^{\ast }\text{and}\tilde{Y}(t)={\mathrm{Id}}_h{\mathrm{Id}}_{\text{ms}}^h{X}_h^{\text{ms}}({\mathrm{Id}}_{\text{ms}}^h{)}^{\ast }{\mathrm{Id}}_h^{\ast }.$$

$$\Vert \tilde{X}(t){\Vert }_{\mathcal{L}(L^2)}+\Vert \tilde{Y}(t){\Vert }_{\mathcal{L}(L^2)}\le C,$$

$$\lVert\tilde{X}(t)-\tilde{Y}(t)\rVert_{\mathcal{L}(L^{2})}\leq CH^{2}\big{(}\log(H^{-1})+t^{-1}\big{)}.$$

$$X_{h}(t)=\mathrm{e}^{t\mathcal{A}_{h}^{*}}X_{h}(0)\mathrm{e}^{t\mathcal{A}_{h}}+\int_{0}^{t}{\mathrm{e}^{(t-s)\mathcal{A}_{h}^{*}}\Big{(}\mathcal{C}_{h}^{*}\mathcal{Q}\mathcal{C}_{h}-X_{h}(s)\mathcal{B}_{h}\mathcal{R}^{-1}\mathcal{B}_{h}^{*}X_{h}(s)\Big{)}\mathrm{e}^{(t-s)\mathcal{A}_{h}}\,\mathrm{d}s}.$$

Full text

{AMS}

49N10, 65N12, 65N30, 93C20

1 Introduction

In a linear quadratic regulator (LQR) problem, the state $x$ is a model of a system whose evolution can be influenced through the input $u$. The goal is to drive certain measurable quantities of the system, the output $y$, to a given target which is typically zero. The relations between $x$, $u$ and $y$ are given by the state and output equations

[TABLE]

where $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ are given operators. The optimal input function $u^{*}$ is found by minimizing the cost functional

[TABLE]

where $\mathcal{Q}$, $\mathcal{R}$ and $\mathcal{G}$ are given weighting factors. It can be shown (see e.g. [1, 21]) that $u^{*}$ is given in feedback form as $u^{*}(t)=-\mathcal{R}^{-1}\mathcal{B}^{*}X(T-t)x(t)$, where $X$ is the solution to an operator-valued differential Riccati equation (DRE):

[TABLE]

In the case of a nonzero output target, one additional differential equation for the evolution of $u^{*}$ has to be solved.

In this paper, we consider the case when the operator $\mathcal{A}$ exhibits multiscale behaviour. In particular, we consider diffusion problems where the spatial variation of the diffusion coefficient is on a fine scale compared to the computational domain. This e.g. occurs in the modeling of composite materials and flows in porous media. Numerically approximating the solutions to elliptic or parabolic equations given by such operators in the usual way is difficult, because a very fine discretization is necessary to resolve the fine-scale structure. These difficulties are exacerbated when considering DREs such as Eq. 3, as their solution essentially requires solving many parabolic equations.

A by now well established method for multiscale elliptic and parabolic problems is the localized orthogonal decomposition (LOD) [25, 15]. It is a modification of the finite element method (FEM), which incorporates some of the fine-scale structure into a coarse discretization by precomputing a series of localized fine-scale problems. Due to the localization, these are much cheaper to evaluate than the full fine-scale problem and may additionally be solved in parallel.

We note that finite elements were introduced for the approximation of optimal control problems already in the 1970’s, see e.g. [27, 10, 14, 36], and the field has grown much in several different directions since then. When diffusion problems have been considered, the focus has typically been on constant or slowly varying diffusion. Recently, however, also optimal control problems of multiscale type have been considered in e.g. [11, 12, 22]. None of these consider the LOD approach, instead preferring homogenization or asymptotic expansions. Additionally, a common assumption is that the multiscale features are periodic, which is frequently not the case in applications.

The focus in this paper is on the approximation of DREs such as Eq. 3. In contrast to the forward-adjoint approach, which solves a specific optimal control problem, the DRE provides the feedback laws for all problems defined by the operators $\mathcal{A},\mathcal{B},\mathcal{C}$. While more expensive to solve, it can be precomputed and reused in many different situations. We refer to [8, 21] for an overview of Riccati theory, with the latter reference treating very general problems.

Our main result is that LOD-approximations to the solution of Eq. 3 with a mesh size $H$ converge with order $H^{2}\log(H^{-1})$ in the $L^{2}$ operator norm to a given accurate fine-scale FEM approximation. The convergence in the corresponding operator energy norm is shown to be of order $H$. We note that $H^{2}\log(H^{-1})$-convergence of FEM approximations to the exact solution of Eq. 3 has previously been shown in [18], and similar results for algebraic Riccati equations can be found in [21]. (See also [30, 6] for convergence results without orders in related settings.) However, the error constants in these results depend on the multiscale variations of $\mathcal{A}$, and thus such convergence is not observed in practice. This is not the case for our present results.

For practical computations, also a temporal discretization is necessary; for this we consider a low-rank splitting scheme as introduced in [32]. Such methods decompose the DRE into its affine and nonlinear parts and approximate these separately, thereby greatly reducing the computational cost. The affine problem requires the approximation of several parabolic equations involving $\mathcal{A}$ in each time step. As the computational efficiency gain for LOD increases with the number of times the modified basis may be reused, splitting schemes are thus particularly well suited to be combined with the LOD method.

We demonstrate how to transform the FEM and LOD discretizations into matrix-valued equations, and how to implement the fully discrete methods. Even if LOD reduces the need for very fine discretizations, large 2D or 3D-problems may still yield large matrices. We therefore consider the low-rank approach, which greatly reduces the necessary amount of computations. As a side effect, this also allows us to compute errors in the operator norms very efficiently.

A brief outline of the paper is as follows: We formalize the setting and our basic assumptions in Section 2, and define the different spatial discretizations in Section 3. Convergence of the LOD approximations with the appropriate order is then shown in Section 4. The matrix-valued formulations of the discretized DREs and related questions are discussed in Section 5, while Section 6 is devoted to the temporal discretization and low-rank setting. Finally, we present several numerical experiments and their results in Section 7.

2 Setting

Let $\Omega\in\mathbb{R}^{d}$, $d\leq 3$, be a bounded polygonal/polyhedral domain. We consider the separable Hilbert spaces $L^{2}(\Omega)$, $V=H^{1}_{0}(\Omega)$, $U$ and $Z$, where $L^{2}(\Omega)$ corresponds to the state space, $U$ is the control space and $Z$ is the observation space. In the following, the specification of $\Omega$ will be omitted. We write $\left(\cdot,\cdot\right)$ and $\lVert\cdot\rVert$ for the inner product and norm on $L^{2}$, and denote the corresponding quantities on $V$, $U$ and $Z$ by subscripts. To define the state evolution operator $\mathcal{A}$, we assume that the inner product $a(u,v)=\int\kappa\nabla u\cdot\nabla v$ on $V\times V$ is given, with assumptions on $\kappa$ given below. Then $\mathcal{A}\colon L^{2}\supset\mathcal{D}(\mathcal{A})\to L^{2}$ is defined by $\left(\mathcal{A}u,v\right)=-a(u,v)$ and $\mathcal{D}(\mathcal{A})=\{u\in V\;|\;\mathcal{A}u\in L^{2}\}$.

Further, let the input operator $\mathcal{B}\colon U\to L^{2}$ and the output operator $\mathcal{C}\colon L^{2}\to Z$ be given. We also consider the output and input weighting operators $\mathcal{Q}\colon Z\to Z$ and $\mathcal{R}\colon U\to U$ (which could be included in $\mathcal{C}$ and $\mathcal{B}$ but are typically not) and the final state weighting operator $\mathcal{G}\colon L^{2}\to L^{2}$. By ∗, we denote Hilbert-adjoint operators with respect to $L^{2}$, so that e.g. $\mathcal{B}^{*}\colon L^{2}\to U$ satisfies $\left(\mathcal{B}x,y\right)=\left(x,\mathcal{B}^{*}y\right)$ for all $x\in U$ and $y\in L^{2}$. Finally, we denote the linear bounded operators from one generic Hilbert space, $Y$, to another, $W$, by $\mathcal{L}(Y,W)$. When $W=Y$, we abbreviate $\mathcal{L}(Y)=\mathcal{L}(Y,Y)$.

In this notation, the weak form of Eq. 3 is to find $X\in\mathcal{L}(L^{2})$ satisfying

[TABLE]

for all $x,y\in\mathcal{D}(\mathcal{A})$.

{assumption}

The diffusion coefficient $\kappa\in L^{\infty}(\mathbb{R}^{d\times d})$ is symmetric and satisfies

[TABLE]

In addition, $\mathcal{B}\in\mathcal{L}(U,L^{2})$, $\mathcal{C}\in\mathcal{L}(L^{2},Z)$, $\mathcal{Q}\in\mathcal{L}(Z)$, $\mathcal{R}\in\mathcal{L}(U)$ is invertible with $\mathcal{R}^{-1}\in\mathcal{L}(U)$ and $X(0)=\mathcal{G}\in\mathcal{L}(L^{2})$.

The first part of Section 2 shows that $a$ is a bounded and coercive bilinear form, which means that $\mathcal{A}$ is the generator of an analytic semigroup $\mathrm{e}^{t\mathcal{A}}\colon L^{2}\to L^{2}$, see e.g. [34, Theorem 3.6.1]. In conjunction with the boundedness assumptions on $\mathcal{B}$, $\mathcal{C}$, $\mathcal{Q}$ and $\mathcal{R}$, this guarantees the existence and uniqueness of a solution to Eq. 4. In fact, there is even a classical solution to Eq. 3 [8, Part IV, Ch. 3], which means that the $\mathcal{A}^{*}X+X\mathcal{A}$ term can be extended to an operator in $\mathcal{L}(L^{2})$. As a consequence, Equation Eq. 4 holds also for $x,y\in L^{2}$. We note that these conclusions are valid also under various weaker forms of Section 2, which additionally permit the treatment of boundary control and observation [21]. A discussion on an extension of our results to such a setting may be found in Section 8.1.

3 Spatial discretization

We first introduce the FEM approximation of Eq. 4. To this end, we let $\mathcal{T}_{h}$ be a triangulation of $\Omega$ with meshwidth $h$ and $N_{h}$ internal nodes. The subspace $V_{h}\subset V$ denotes the space of continuous and piecewise affine functions on $\mathcal{T}_{h}$, and we denote the corresponding nodal basis functions by $\{\varphi^{h}_{i}\}_{i=1}^{N_{h}}$. This discretization is referred to as the fine, or sometimes also reference, mesh, see further Section 3.1 below.

We also consider a coarse discretization space $V_{H}\subset V_{h}$ for $H>h$, with the corresponding family of triangulations $\{\mathcal{T}_{H}\}_{H>h}$, which is assumed to be quasi-uniform. For these triangulations, we let $B_{K}$ be the largest ball contained in the triangle $K$ and denote by $\gamma>0$ the shape regularity of the mesh, defined by

[TABLE]

Furthermore, we let $\operatorname{Id}_{H}^{h}\colon V_{H}\to V_{h}$ denote the identity operator between these spaces, i.e. $\operatorname{Id}_{H}^{h}u=u$ for all $u\in V_{H}$. Similarly, $\operatorname{Id}_{h}\colon V_{h}\to L^{2}$ is the identity operator mapping into $L^{2}$ and $\operatorname{Id}_{h}^{*}$ is the $L^{2}$-orthogonal projection of $L^{2}$ onto $V_{h}$.

The semi-discretized weak form of Eq. 3 is defined by

[TABLE]

for all $x,y\in V_{h}$ and with $X_{h}\colon V_{h}\to V_{h}$ satisfying $X_{h}(0)=\operatorname{Id}_{h}^{*}X(0)\operatorname{Id}_{h}$. Here, the operators $\mathcal{A}_{h}\colon V_{h}\to V_{h}$, $\mathcal{B}_{h}\colon U\to V_{h}$ and $\mathcal{C}_{h}\colon V_{h}\to Z$ satisfy

[TABLE]

for all $x,y\in V_{h}$, $u\in U$ and $z\in Z$. We note that $X$ can be proven to be self-adjoint, so we additionally require that $X_{h}$ is self-adjoint.

For the coarse discretization, we have the same equation but with $H$ instead of $h$. We observe that the coarse and fine operators are related in the following way:

[TABLE]

and that the natural extension of $X_{H}$ to a map on $V_{h}$ is given by $\operatorname{Id}_{H}^{h}X_{H}(\operatorname{Id}_{H}^{h})^{*}$. Here, $(\operatorname{Id}_{H}^{h})^{*}$ is the $L^{2}$-orthogonal projection of $V_{h}$ onto $V_{H}$.

3.1 Localized orthogonal decomposition

If $\kappa$ is varying on a small scale of size $\epsilon>0$, then the classical FEM approximation of a parabolic problem $\dot{x}=\mathcal{A}x+f$ may yield poor results, unless $h$ is sufficiently small to resolve the fine-scale variations. That is, we typically do not observe $O(h^{2})$-convergence until $h<\epsilon$, which requires infeasible amounts of computation. The same behaviour occurs for the $X_{h}$-discretizations of Eq. 4.

To this end, we assume that $h$ is sufficiently small so that $X_{h}$ is a good approximation of $X$. That is, $h<\epsilon$, and we refer to $X_{h}$ as the reference solution. The aim is now to approximate $X_{h}$ by using a multiscale space $V_{\text{ms}}\subset V_{h}$ of the same dimension as the coarse space $V_{H}$. To obtain such a space, we use the localized orthogonal decomposition (LOD) method introduced in [25], which incorporates fine-scale information in the coarse-scale space. The construction involves the solution of several fine-scale, but localized and parallelizable, problems. We briefly summarize the procedure here and refer to [25, 15], for the details.

To define the multiscale space $V_{\text{ms}}$, we first introduce an interpolation operator $I_{H}\colon V_{h}\rightarrow V_{H}$ that fulfills

[TABLE]

for all triangles $K\in\mathcal{T}_{H}$, where $\omega_{K}:=\cup\{\hat{K}\in\mathcal{T}_{H}:\hat{K}\cap K\neq\emptyset\}$. In this paper we use the weighted Clément interpolant as in [25]. Let $V_{\text{f}}$ denote the kernel of $I_{H}$,

[TABLE]

and note that $V_{h}$ can be decomposed as $V_{h}=V_{H}\oplus V_{\text{f}}$, meaning that every $v_{h}\in V_{h}$ can be written as $v_{h}=v_{H}+v_{\text{f}}$ with $v_{H}\in V_{H},v_{\text{f}}\in V_{\text{f}}$. We now introduce the (global) correction operator $\hat{Q}_{h}\colon V_{H}\rightarrow V_{\text{f}}$ by

[TABLE]

and define the (global) multiscale space as $\hat{V}_{\text{ms}}:=\hat{R}_{h}V_{H}=V_{H}-\hat{Q}_{h}V_{H}$, with $\hat{R}_{h}:=\operatorname{Id}_{H}^{h}-\hat{Q}_{h}$. This leads to the decomposition $V_{h}=\hat{V}_{\text{ms}}\oplus V_{\text{f}}$ with the orthogonality $a(\hat{v}_{\text{ms}},v_{\text{f}})=0$, $\hat{v}_{\text{ms}}\in\hat{V}_{\text{ms}}$, $v_{\text{f}}\in V_{\text{f}}$. Note that $\hat{Q}_{h}$ is the orthogonal projection onto $V_{\text{f}}$ with respect to the inner product $a(\cdot,\cdot)$, i.e. the Ritz projection onto $V_{\text{f}}$, and $\hat{V}_{\text{ms}}$ is the orthogonal complement to $V_{\text{f}}$. From the construction it follows that $\dim\hat{V}_{\text{ms}}=\dim V_{H}$. Indeed, a basis for $\hat{V}_{\text{ms}}$ is given by $\{\varphi^{H}_{i}-\hat{Q}_{h}\varphi^{H}_{i}:i=1,...,N_{H}\}$.

In general, the corrections $\hat{Q}_{h}\varphi^{H}_{i}$ have global support and are expensive to compute, since they are posed in the entire fine scale space $V_{\text{f}}\subseteq V_{h}$. To overcome this, it is observed that the corrections have exponential decay away from the $i$:th node of $\mathcal{T}_{H}$ (see [25, 15]), which motivates a truncation of the corrections. For this purpose, we define patches $\omega_{k}(K)$ of size $k$ around each $K\in\mathcal{T}_{H}$ by the following:

[TABLE]

Further, we define $V^{K}_{\text{f}}:=\{v\in V_{\text{f}}:v(z)=0\text{ on }\overline{\Omega}\setminus\omega_{k}(K)\}$ to be the restriction of $V_{\text{f}}$ to the patch $\omega_{k}(K)$. For brevity, we do not include the dependence on $k$ in the notation. Now note that the correction operator $\hat{Q}_{h}$ can be written as the sum $\hat{Q}_{h}=\sum_{K\in\mathcal{T}_{H}}\hat{Q}_{h}^{K},$ where

[TABLE]

We can now localize these computations by replacing $V_{\text{f}}$ with $V_{\text{f}}^{K}$. Define $Q_{h}^{K}\colon V_{H}\rightarrow V_{\text{f}}^{K}$ such that

[TABLE]

Finally, we can define a local operator $Q_{h}:=\sum_{K\in\mathcal{T}_{H}}Q_{h}^{K}$ and a localized space $V_{\text{ms}}:=R_{h}V_{H}=V_{H}-Q_{h}V_{H}$, with $R_{h}:=\operatorname{Id}_{H}^{h}-Q_{h}$.

The approximation properties (and the required computational effort) of the space $V_{\text{ms}}$ depends on the choice of $k$. In [15] it is proven that convergence of order $H^{2}$ is obtained if $k$ is chosen proportional to $\log H^{-1}$. In this paper we therefore assume that $k\sim\log H^{-1}$ to avoid explicitly stating the dependence on $k$.

To define an LOD-approximation to the solution $X_{h}$ in Eq. 5, we additionally need to introduce the identity operator $\operatorname{Id}_{\text{ms}}^{h}:V_{\text{ms}}\to V_{h}$, $\operatorname{Id}_{\text{ms}}^{h}u=u$. Its $L^{2}$-adjoint is the $L^{2}$-orthogonal projection of $V_{h}$ onto $V_{\text{ms}}$. Replacing the space $V_{h}$ with $V_{\text{ms}}$ then results in the problem to find $X^{\text{ms}}_{h}\colon V_{\text{ms}}\to V_{\text{ms}}$ satisfying

[TABLE]

for all $u,v\in V_{\text{ms}}$ and the initial condition $X^{\text{ms}}_{h}(0)=(\operatorname{Id}_{\text{ms}}^{h})^{*}\operatorname{Id}_{h}^{*}X(0)\operatorname{Id}_{h}\operatorname{Id}_{\text{ms}}^{h}$. Here, the operators $\mathcal{A}^{\text{ms}}_{h}\colon V_{\text{ms}}\to V_{\text{ms}}$, $\mathcal{B}^{\text{ms}}_{h}\colon U\to V_{\text{ms}}$ and $\mathcal{C}^{\text{ms}}_{h}\colon V_{\text{ms}}\to Z$ are given by

[TABLE]

for all $v,w\in V_{\text{ms}}$, $u\in U$ and $z\in Z$. Similar to Eq. 6 we have

[TABLE]

The natural $V_{h}$-extension of $X^{\text{ms}}_{h}$ is given by $\operatorname{Id}_{\text{ms}}^{h}X^{\text{ms}}_{h}(\operatorname{Id}_{\text{ms}}^{h})^{*}$, similar to the $X_{H}$-case.

Since $V_{\text{ms}}$ has the same dimension as $V_{H}$, there is a lower-dimensional representative for $X^{\text{ms}}_{h}$, given by $X^{\text{ms}}_{h}=R_{h}X^{\text{ms}}_{H}R_{h}^{-1}$. By inserting $u=R_{h}x$ and $v=R_{h}y$, with $x,y\in V_{H}$, in Eq. 7 we see that

[TABLE]

and we consequently define the corrected coarse-scale operators

[TABLE]

4 Error analysis

In the following, $C$ denotes a generic constant which may take different values at different occasions. It may depend on the problem data and the size of the domain, but is independent of $h$ and $H$. Moreover, it does not depend on the multiscale variations of $\mathcal{A}$, i.e. any derivatives of $\kappa$. We start by gathering some useful results:

4.1 Preliminaries

Recall that $\operatorname{Id}_{h}\colon V_{h}\to L^{2}$ is the identity mapping, $P_{h}=\operatorname{Id}_{h}^{*}\colon L^{2}\to V_{h}$ denotes the $L^{2}$-orthogonal projection onto $V_{h}$, and $P_{\text{ms}}$ is the $L^{2}$-orthogonal projection onto $V_{\text{ms}}$. We have $P_{\text{ms}}=(\operatorname{Id}_{\text{ms}}^{h})^{*}P_{h}$, i.e. we first project onto $V_{h}$ and then onto $V_{\text{ms}}$. Straightforward calculations show the following:

Lemma 4.1**.**

*Under Section 2, it holds that $\operatorname{Id}_{h}\mathcal{B}_{h}\in\mathcal{L}(U,L^{2})$, $\mathcal{C}_{h}P_{h}\in\mathcal{L}(L^{2},Z)$ and $S_{h}:=\operatorname{Id}_{h}\mathcal{B}_{h}\mathcal{R}^{-1}\mathcal{B}_{h}^{*}\operatorname{Id}_{h}^{*}\in\mathcal{L}(L^{2})$. *

Further, let $\mathrm{e}^{t\mathcal{A}_{h}}$ denote the solution operator to the equation $\dot{u}+\mathcal{A}_{h}u=0$, i.e. the semigroup generated by $\mathcal{A}_{h}$. Similarly, $\mathrm{e}^{t\mathcal{A}^{\text{ms}}_{h}}$ is the semigroup generated by $\mathcal{A}^{\text{ms}}_{h}$. Because $\mathcal{A}$ generates an analytic semigroup on $L^{2}$, these operators are analytic semigroups on $V_{h}$ and $V_{\text{ms}}$, respectively. More specifically, we have

Lemma 4.2**.**

Under Section 2, the operators

[TABLE]

*are both in $\mathcal{L}(L^{2})$ for $t\in[0,T]$, with the uniform bounds $\lVert E_{h}(t)\rVert_{\mathcal{L}(L^{2})}\leq 1$ and $\lVert E_{\text{ms}}(t)\rVert_{\mathcal{L}(L^{2})}\leq 1$. *

By arguing as in [24], but for the (simpler) semi-discrete case, we have (choosing $k\sim\log H)$

Lemma 4.3**.**

For $t\in(0,T]$ it holds that

[TABLE]

*Here, the constant $C$ depends on $T$, $\alpha$, $\beta$, and $\gamma$, but not on the multiscale variations of $\mathcal{A}$. *

Proof 4.4**.**

We only comment briefly on the proof here. Let $u_{h}(t)=e^{t\mathcal{A}_{h}}P_{h}v$ and $u_{\text{ms}}(t)=e^{t\mathcal{A}^{\text{ms}}_{h}}P_{\text{ms}}v$, for $v\in L^{2}(\Omega)$. By introducing the Ritz projection $R_{\text{ms}}\colon V_{h}\rightarrow V_{\text{ms}}$ satisfying $a(R_{\text{ms}}v,w)=a(v,w)$ for all $w\in V_{\text{ms}},v\in V_{h}$ we get, see [35, Chapter 3] and [24],

[TABLE]

where $\rho:=u_{h}-R_{\text{ms}}u_{h}$. From the error bounds of $R_{\text{ms}}$ in [25], see also [24], we get

[TABLE]

*The regularity estimates $\lVert D_{t}^{l}u_{h}(t)\rVert\leq Ct^{-l}\lVert v\rVert$, for $t=0,1,2$, [35, Lemma 2.5], completes the proof. *

Finally, from Lemma 4.1, we get the existence and uniqueness of solutions $X_{h}$ and $X^{\text{ms}}_{h}$ to the discretized DREs Eq. 5 and Eq. 7, respectively. Let us abbreviate

[TABLE]

Then we have

Lemma 4.5**.**

There is a constant $C>0$ which is independent of the multiscale variations of $\mathcal{A}$ but may depend on $\alpha$ and $\beta$, such that

[TABLE]

*for $t\in[0,T]$. *

4.2 Error analysis

We are now ready for the main theorem of this paper:

Theorem 4.6**.**

Suppose that Section 2 is fulfilled. Then for $t\in(0,T]$ it holds that

[TABLE]

*Here, the constant $C$ depends on $T$, $\alpha$, $\beta$, $\gamma$, and $\lVert X(0)\rVert_{\mathcal{L}(L^{2})}$, but not on the multiscale variations of $\mathcal{A}$. *

Proof 4.7**.**

We utilize the integral form of Eq. 5. If $X_{h}$ solves Eq. 5 then it satisfies

[TABLE]

(see e.g. [8, Chapter IV:3, Proposition 2.1]). Recalling that $\operatorname{Id}_{h}^{*}\operatorname{Id}_{h}$ and $(\operatorname{Id}_{\text{ms}}^{h})^{*}\operatorname{Id}_{\text{ms}}^{h}$ are the identity operators on $V_{h}$ and $V_{\text{ms}}$, respectively, and using Eq. 8 therefore shows that

[TABLE]

as well as

[TABLE]

(Note the $\tilde{X}(0)$ in the first term, since we suppose $X^{\text{ms}}_{h}(0)=(\operatorname{Id}_{\text{ms}}^{h})^{*}X_{h}(0)\operatorname{Id}_{\text{ms}}^{h}$.)

Subtracting these expressions yields

[TABLE]

so that

[TABLE]

We observe that for all $G\colon L^{2}\to Y$ (with a generic Hilbert space $Y$) it holds that $\lVert G\rVert_{\mathcal{L}(L^{2},Y)}=\lVert G^{*}\rVert_{\mathcal{L}(Y,L^{2})}$. Thus, using Lemmas 4.2, 4.3 and 4.5 we get

[TABLE]

Additionally using Lemma 4.1 shows that the last two integrands satisfy

[TABLE]

Due to the singularity at $s=t$ in the bound on $\lVert E_{h}(t-s)-E_{\text{ms}}(t-s)\rVert_{\mathcal{L}(L^{2})}$, we split the integrals of the remaining $R_{j}$-terms into two parts. For $R_{3}$, we find

[TABLE]

where we have used $t\leq T$ for the crude estimate $\log t\leq Ct^{-1}$, since a $t^{-1}$-term already appears in the bounds of $R_{1}$ and $R_{2}$. The same bound holds for $R_{4}$, and, by Lemma 4.1 and Lemma 4.5, also for $R_{5}$ and $R_{6}$. In conclusion, we thus have

[TABLE]

*which by Grönwall’s lemma yields the statement of the theorem. *

Remark 4.8**.**

*In the common situation that $X(0)=0$, corresponding to the case of no final state penalization, the $t^{-1}$-singularity disappears. *

Remark 4.9**.**

*We note that a bound of the same form has been shown in [18] for the FEM error. However, the error constant then depends on the variations in $\kappa$, and one does not observe the given convergence order until $H<\epsilon$. *

Similar to the parabolic case, the error bound becomes less singular near $t=0$ if we measure in the $V$-norm. To prove this we need the following, slightly stronger, assumptions on the operators (cf. Section 2):

{assumption}

In addition to Section 2, $\mathcal{B}\in\mathcal{L}(U,V)$, $\mathcal{C}\in\mathcal{L}(V,Z)$, and $X(0)=\mathcal{G}\in\mathcal{L}(V)$. Moreover, we assume that the mesh $\mathcal{T}_{h}$ is of a form such that $P_{h}$ is stable in the $V$-norm.

Remark 4.10**.**

*In particular, quasi-uniform meshes satisfy Section 4.2. We refer to [2] for a discussion on more general permissible meshes. *

Theorem 4.11**.**

Suppose that Section 4.2 is fulfilled. For $t\in(0,T]$ it holds that

[TABLE]

*Here, the constant $C$ depends on $T$, $\alpha$, $\beta$, $\gamma$, and $\lVert X(0)\rVert_{\mathcal{L}(V)}$, but not on the multiscale variations of $\mathcal{A}$. *

Proof 4.12**.**

We start by noting that $\lVert\operatorname{Id}_{h}\rVert_{\mathcal{L}(V_{h},V)}\leq 1$. Furthermore, since $P_{h}$ is stable in the $V$-norm, the following bound holds

[TABLE]

Now, note that if the initial data $v\in V$ then we may instead of Lemma 4.3 prove the following, less singular, error bound

[TABLE]

In addition, parabolic regularity gives the bounds $\lVert E_{h}(t)\rVert_{\mathcal{L}(V)},\lVert E_{\text{ms}}(t)\rVert_{\mathcal{L}(V)}\leq C$.

Note that $\lVert\mathcal{C}_{h}\rVert_{\mathcal{L}(V_{h},Z)}\leq\lVert\mathcal{C}\rVert_{\mathcal{L}(V,Z)}$, so from Section 4.2 it follows that

[TABLE]

Similarly, $\lVert\mathcal{B}_{h}\rVert_{\mathcal{L}(U,V_{h})}\leq\lVert\mathcal{B}\rVert_{\mathcal{L}(U,V)}$, and we have

[TABLE]

As in the proof of Theorem 4.6 we can write the difference $\tilde{X}(t)-\tilde{Y}(t)$ as a sum of eight terms so that

[TABLE]

For $R_{1}$ we have

[TABLE]

and similarly we prove $\lVert R_{2}\rVert_{\mathcal{L}(V)}\leq CHt^{-1/2}$, where we have used that

[TABLE]

which is bounded due to Section 4.2. Using the bounds Eqs. 10 and 11 we get

[TABLE]

and

[TABLE]

*By applying Grönwall’s lemma we obtain the desired error bound. *

5 Matrix-valued formulation

To perform actual computations, we write the finite-dimensional equations on matrix form by expressing the equations in the FEM or LOD bases. To this end, let the function $x\in V_{h}$ and the operator $X_{h}\colon V_{h}\to V_{h}$ have the vector and matrix representations $\bm{x}\in\mathbb{R}^{N_{h}}$ and $\tilde{\bm{X}}^{h}\in\mathbb{R}^{N_{h}\times N_{h}}$, i.e.

[TABLE]

Since exactly the same results hold for $V_{H}$ upon replacing $h$ by $H$, we frequently omit the $h$ sub- and superscripts in the following manipulations. They will be reinstated later when we compare different discretizations. The coordinates satisfy

[TABLE]

where $\bm{M}$ denotes the (symmetric) mass matrix, $\bm{M}_{i,j}=\left(\varphi_{j},\varphi_{i}\right)$. Unfortunately, we will not recover the usual form of the matrix-valued DRE when working in these coordinates. Therefore, we perform the change of variables

[TABLE]

Coincidentally, this means that we actually have

[TABLE]

Equation 5 is equivalent to

[TABLE]

for $1\leq i,j\leq N$, and since $X_{h}\varphi_{i}=\sum_{k=1}^{N}{(\bm{X}\bm{M})_{k,i}\varphi_{k}}$, the first term becomes

[TABLE]

Likewise, with the (negative) stiffness matrix $\bm{A}_{i,j}=\left(\mathcal{A}\varphi_{j},\varphi_{i}\right)$, the second and third terms become

[TABLE]

due to the symmetry of $\bm{M}$ and $\bm{X}$. (Recall that we search for a self-adjoint operator $X_{h}$.) Finally, the last two terms can be written

[TABLE]

where $\bm{B}_{i,j}=\left(\mathcal{B}\varphi^{U}_{j},\varphi_{i}\right)$, $\bm{Q}_{i,j}=\left(\mathcal{Q}\varphi^{Z}_{j},\varphi^{Z}_{i}\right)$, $\bm{R}_{i,j}=\left(\mathcal{R}\varphi^{U}_{j},\varphi^{U}_{i}\right)$, $\bm{C}_{i,j}=\left(\mathcal{C}\varphi_{j},\varphi^{Z}_{i}\right)$ and $\{\varphi^{U}_{i}\}$, $\{\varphi^{Z}_{i}\}$ denote orthonormal bases for $U$ and $Z$, respectively. Summarizing, we can write the equation on matrix form as

[TABLE]

Similar to the relations between the fine and coarse operators Eq. 6, it is easily shown that their matrix representations satisfy

[TABLE]

where $\bm{I}_{H}^{h}\in\mathbb{R}^{N_{h}\times N_{H}}$ is the prolongation matrix that satisfies $\bm{I}_{H}^{h}\bm{x}^{H}=\bm{x}^{h}$ if $x=\sum_{j=1}^{N_{H}}{\bm{x}^{H}_{j}\varphi^{H}_{j}}$ and $\operatorname{Id}_{H}^{h}x=\sum_{j=1}^{N_{h}}{\bm{x}^{h}_{j}\varphi^{h}_{j}}$. By expressing the $\varphi^{H}$ functions in terms of $\varphi^{h}$, it can be seen that $(\bm{I}_{H}^{h})_{i,j}=\varphi^{H}_{j}(z_{i})$, where $z_{i}$ is the $i$:th node of $\mathcal{T}_{h}$. Thus the coarse systems are easily constructed when the fine system is known. Note, however, that the matrix representation of $(\operatorname{Id}_{H}^{h})^{*}$ is not $(\bm{I}_{H}^{h})^{T}$ but $\bm{M}_{H}^{-1}(\bm{I}_{H}^{h})^{T}\bm{M}_{h}$.

For the LOD case, we let $\bm{Q}_{h}$ and $\bm{R}_{h}=\bm{I}_{H}^{h}-\bm{Q}_{h}$ be the matrix representations of $Q_{h}$ and $R_{h}$, respectively. To compute them efficiently, we follow [13]. Then

[TABLE]

where $\bm{X}^{\text{ms}}_{H}$ is symmetric and satisfies

[TABLE]

with the matrices

[TABLE]

Finally, we note that if $u\in V_{h}$, $w\in V_{H}$ and $(\operatorname{Id}_{\text{ms}}^{h})^{*}u=R_{h}w$, then in coordinates we have $\bm{w}=\bm{M}_{\text{ms}}^{-1}\bm{R}_{h}^{T}\bm{M}_{h}\bm{u}$. This means that the matrix representation of $\operatorname{Id}_{\text{ms}}^{h}X^{\text{ms}}_{h}(\operatorname{Id}_{\text{ms}}^{h})^{*}$ is $\bm{R}_{h}\bm{X}^{\text{ms}}_{H}\bm{R}_{h}^{T}\bm{M}_{h}$.

5.1 Error computation

We measure the quality of different approximations as the $\mathcal{L}(L^{2})$-normed distance to a reference approximation at the final time $T$. In order to find a matrix representation for this, we first observe that since $\lVert P_{h}x\rVert\leq\lVert x\rVert$, we have

[TABLE]

But $P_{h}x=x$ for $x\in V_{h}$, so since $V_{h}\subset L^{2}$ we also get

[TABLE]

To compute the $\mathcal{L}(L^{2})$-norm it is thus enough to test with $x=\sum_{i=1}^{N_{h}}{\bm{x}_{i}\varphi^{h}_{i}}\in V_{h}$. Again omitting the $h$ sub- and superscripts, we have that $\left(x,x\right)=\bm{x}^{T}\bm{M}\bm{x}$, and similarly

[TABLE]

Since $\bm{M}$ is symmetric positive definite, we may do a Cholesky factorization $\bm{M}=\bm{L}_{\bm{M}}\bm{L}_{\bm{M}}^{T}$, and the change of variables $\bm{y}=\bm{L}_{\bm{M}}^{T}\bm{x}$ yields

[TABLE]

where $\lVert\cdot\rVert_{\mathbb{R}^{N\times N}}$ denotes the standard spectral matrix norm. Recalling the matrix representation $\bm{I}_{H}^{h}$ of $\operatorname{Id}_{H}^{h}$, we now get that

[TABLE]

The LOD error is completely analogous, using instead $\bm{R}_{h}$ and $\bm{X}^{\text{ms}}_{H}$.

A similar approach also allows us to compute $\mathcal{L}(V)$-errors. Let $\bm{A}=\bm{L}_{\bm{A}}\bm{L}_{\bm{A}}^{T}$ be a Cholesky factorization of $\bm{A}$. Then

[TABLE]

We also get that $\lVert\operatorname{Id}_{h}X_{h}P_{h}-\operatorname{Id}_{H}X_{H}P_{H}\rVert_{\mathcal{L}(V)}$ is bounded from below by $\lVert\bm{L}_{\bm{A}}^{T}(\bm{X}_{h}-\bm{I}_{H}^{h}\bm{X}_{H}(\bm{I}_{H}^{h})^{T})\bm{M}\bm{L}_{\bm{A}}^{-T}\rVert_{\mathbb{R}^{N\times N}}$, i.e. the latter quantity can be thought of as an equivalent norm. Since $\bm{L}_{\bm{A}}$ is triangular, the extra cost required for the computation of $\bm{L}_{\bm{A}}^{-T}$ is negligible. If the low-rank formulation is used (see Section 6.1), only a small number of linear equation systems involving $\bm{L}_{\bm{A}}$ needs to be solved, reducing the cost even further.

6 Temporal discretization

We discretize the matrix-valued DREs in time by means of a low-rank splitting scheme, since the basic operation in such methods is the application of $\mathrm{e}^{t\bm{A}^{T}}$, i.e. essentially solving a parabolic problem. Let $\tau$ denote a fixed time step, and let $t_{j}=j\tau$, $j=0,\ldots,N_{t}$, be the time discretization of the interval $[0,T]$. We split Equation Eq. 15 into two parts, $\dot{\bm{X}}=\mathscr{F}\bm{X}+\mathscr{G}\bm{X}$, where

[TABLE]

Then the Strang splitting approximation at time $t_{j}$ is given by $\bm{X}^{j}$, with $\bm{X}^{0}=\bm{X}(0)$ and

[TABLE]

Here, the solution operators $\mathrm{e}^{t\mathscr{F}}$ and $\mathrm{e}^{t\mathscr{G}}$ satisfy

[TABLE]

where the first equality is apparent from the integral formulation Eq. 9, while the second is easily verified by differentiation.

The low-rank version of the method relies on the assumption that the solution $\bm{X}$ has low rank. This is general true for LQR problems and dramatically reduces the computational cost. In that case, we may factorize $\bm{X}=\bm{L}\bm{D}\bm{L}^{T}$, where $\bm{L}\in\mathbb{R}^{N_{h}\times r}$ and $\bm{D}\in\mathbb{R}^{r\times r}$ with the rank $r\ll N_{h}$. Also $\mathrm{e}^{\tau\mathscr{F}}\bm{X}$ and $\mathrm{e}^{\tau\mathscr{G}}\bm{X}$, and thus the iterates $\bm{X}_{j}$, may then be factorized in such a way. After a reformulation, $\mathrm{e}^{\tau\mathscr{G}}\bm{X}$ is very cheap to compute, and the computation of $\mathrm{e}^{\tau\mathscr{F}}\bm{X}$ reduces to an evaluation of $\mathrm{e}^{\tau\bm{M}^{-T}\bm{A}^{T}}\bm{L}$ (plus preliminary, similar work for the integral term). The latter operation is equivalent to solving $\bm{M}\dot{x}=\bm{A}^{T}x$, $x(0)=\bm{L}$, and the matrix $\bm{M}$ is thus never explicitly inverted. For further details, we refer to [32, 33].

6.1 Low-rank errors

Also the error computations outlined in Section 5.1 benefit from being formulated in a low-rank setting. Assume that $\bm{X}_{h}=\bm{L}_{h}\bm{D}_{h}\bm{L}_{h}^{T}$ and $\bm{X}_{H}=\bm{L}_{H}\bm{D}_{H}\bm{L}_{H}^{T}$ with $\bm{L}_{h}\in\mathbb{R}^{N_{h}\times r_{h}}$ and $\bm{L}_{H}\in\mathbb{R}^{N_{H}\times r_{H}}$ with $r_{h},r_{H}\ll N_{h}$, and let $\bm{M}_{h}=\bm{L}_{\bm{M}}\bm{L}_{\bm{M}}^{T}$ be a Cholesky factorization. By setting

[TABLE]

we see that $\bm{V}\in\mathbb{R}^{N_{h}\times(r_{h}+r_{H})}$, $\bm{D}\in\mathbb{R}^{(r_{h}+r_{H})\times(r_{h}+r_{H})}$ and it follows that

[TABLE]

Since $\bm{V}\bm{D}\bm{V}^{T}$ is not necessarily an eigenvalue decomposition, we cannot immediately determine the norm by inspection. However, performing a QR-factorization $\bm{V}=\bm{Q}\bm{R}$ is cheap if the number of columns is low, and $\bm{R}\bm{D}\bm{R}^{T}\in\mathbb{R}^{(r_{h}+r_{H})\times(r_{h}+r_{H})}$ can also be diagonalized cheaply. (This is precisely the $LDL^{T}$ column compression procedure which is applied in each time step.) We acquire $\bm{V}\bm{D}\bm{V}^{T}=(\bm{Q}\bm{W})\tilde{\bm{D}}(\bm{Q}\bm{W})^{T}$, for some $\bm{W}$, where $\lVert\bm{V}\bm{D}\bm{V}^{T}\rVert=|\tilde{\bm{D}}_{1,1}|$.

For errors in the $\mathcal{L}(V)$-norm, we do not get a symmetric matrix as above. But if $\bm{A}=\bm{L}_{\bm{A}}\bm{L}_{\bm{A}}^{T}$ we can still write

[TABLE]

with the same $\bm{D}$, and with

[TABLE]

We can cheaply $QR$-factorize both $\bm{G}_{1}=\bm{U}\bm{R}_{1}$ and $\bm{G}_{2}=\bm{V}\bm{R}_{2}$; this means that

[TABLE]

where $\bm{S}=\bm{R}_{1}\bm{D}\bm{R}_{2}^{T}$ is a small matrix.

7 Numerical experiments

We have performed a number of numerical experiments in order to verify our a priori error bounds for the LOD discretizations, and to demonstrate their efficiency in comparison to the classical FEM.

In all experiments, we compute the relevant matrices for both FEM and LOD by using efficient code written by Fredrik Hellman and Daniel Elfverson111Available on request from Fredrik Hellman, .. These pre-solve computations were run on a Intel® Core™ i5-4690 processor. We note that the localized elliptic fine-scale problems were not solved in parallel. Doing so would further improve the performance of LOD.

For approximating the solutions to the DREs, we employ in all cases the low-rank Strang splitting scheme (as described in Section 6) with $N_{t}=256$ time steps. This ensures that the temporal error is small compared to the spatial error, which is our interest here. Our implementation utilizes the DREsplit222Available on request from Tony Stillfjord, , or from http://www.tonystillfjord.net. library. These computations were performed on resources at Chalmers Centre for Computational Science and Engineering (C3SE) provided by the Swedish National Infrastructure for Computing (SNIC). Each simulation used a single Intel® Xeon® E5-2650 v3 processor.

The multiscale diffusion coefficients $\kappa$ considered in the numerical examples are of two distinct types. In Examples 1, 2, and 4 we consider a piecewise constant coefficient, generated randomly with no spatial correlation, that varies on a fine scale, see Fig. 1. In Examples 3 and 5 $\kappa$ takes two values. One value in the background and one in the thin channels, see Fig. 6. This is a common setup for reinforced (composite) materials. Both these cases are challenging for the finite element method.

7.1 Example 1

In this first example, we consider diffusion on the unit square. More specifically, we take $\Omega=[0,1]^{2}$ and set \mathcal{A}x=\nabla\cdot\big{(}\kappa\nabla x\big{)} with Dirichlet boundary conditions. Here, $\kappa$ is piecewise constant on a square grid of size $2^{-7}$ and taking randomly chosen values in $[10^{-3},1]$; see Fig. 1 for an illustration. We consider $3$ independent inputs and define the input operator $\mathcal{B}$ as the sum

[TABLE]

Thus we can control the system on three small squares. As the output operator we take the mean, i.e. $\mathcal{C}x=\int_{\Omega}{x}$. We choose $\mathcal{Q}$ and $\mathcal{R}$ to be the identity operators and take $\mathcal{G}=X(0)=0$.

For the discretization in space, we start with a coarse mesh containing $8$ triangles, and then refine this $6$ times, giving meshes with $2^{3+2j}$ triangles, for ${j=0,\ldots,6}$. One additional refinement provides the reference grid with $2^{17}=131072$ triangles. This results in matrices $\bm{A}_{j}\in\mathbb{R}^{n\times n}$, $\bm{B}_{j}\in\mathbb{R}^{n\times 3}$ and $\bm{C}_{j}\in\mathbb{R}^{1\times n}$, $j=0,\ldots,7$, with $n=1,9,49,225,961,3969,16129,65025$ (since we only consider the interior nodes).

The approximations are compared only at the final time, in the $\mathcal{L}(L^{2})$- and $\mathcal{L}(V)$-norms as outlined in Section 5.1, and the computed errors are shown in Fig. 2. We see that the classical FEM initially struggles due to not resolving the multiscale coefficient properly, but converges with order $2$ when the mesh becomes fine enough. The LOD approach converges with order $2$ also for the coarse meshes, and additionally results in approximations that are about one order of magnitude more accurate. The plot to the right shows the errors against the actual computation time, including the time spent on constructing the LOD bases. As can be seen, this extra effort is low enough that except for the most inaccurate cases it is always worthwhile to use the LOD approach.

7.2 Example 2

Here, we consider an L-shaped domain $\Omega$, where $[0.5,1]\times[0.5,0.5]$ has been removed from the unit square. The diffusion coefficient $\kappa$ is piecewise constant on a square grid of size $2^{-7}$ and taking random values in $[10^{-3},1]$. We use one control input, given by the characteristic function of the square $[0.65,0.85]^{2}$, and one output, the mean over the square $[0.15,0.35]^{2}$. The meshes are setup as in the previous example, but now with $n=5,33,161,705,2945,12033,48641$ interior nodes ($n=195585$ for the reference solution). The time discretization and other parameters are the same as in the previous example.

The results are shown in Fig. 3. Due to the reentrant corner the errors behave more erratically than in the previous example, but LOD is still clearly first- and second-order convergent in contrast to standard FEM, which performs very poorly. We also observe that LOD is more efficient in all but the coarsest cases.

7.3 Example 3

We again consider the setting of Example 1, but replace the diffusivity constant. Here, $\kappa$ takes the constant value $1$ everywhere, except for in $7$ horizontal stripes where it is $10^{-2}$. The stripes are centered around the heights $j/8$, $j=1,\ldots,7$, and have a width of $2^{-7}$.

The results are shown in Fig. 4. This time, the detrimental effect on the FEM discretization is even more pronounced, with almost no convergence until the thin stripes can be resolved. The LOD approximations are once again more accurate for all $H$. We note that the $\mathcal{L}(L^{2})$-error is not quite $\mathcal{O}(H^{2})$ in this case, but rather close to $\mathcal{O}(H^{2}\log{H^{-1}})$ as predicted by Theorem 4.6. Like in the previous example, computing the LOD bases is cheap enough that the LOD approach is more efficient in all but the least accurate cases.

7.4 Example 4

In this example, we deviate from the basic setting described in Section 4 by considering a boundary control application. All parameters except for the boundary conditions and the input operator are the same as in Example 1. We call the union of the top and bottom edges of the unit square $\Gamma_{D}$ and impose homogeneous Dirichlet boundary conditions there. The left and right edges we denote $\Gamma_{1}$ and $\Gamma_{2}$, respectively, and there we impose nonhomogeneous Neumann boundary conditions. In particular, with the outward-pointing normal denoted by $n$, we consider functions $x$ satisfying

[TABLE]

Here, $u_{1}$ and $u_{2}$ are the two control inputs, and

[TABLE]

is a fixed function. The operator $\mathcal{A}$ now corresponds to x\mapsto\nabla\cdot\big{(}\kappa\nabla x\big{)} on the space $\{x\in H^{1}(\Omega)\;|\;x_{\rvert_{\Gamma_{D}}}=0\}$ with no conditions imposed on $\Gamma_{1}$, $\Gamma_{2}$, while the (unbounded) operator $\mathcal{B}$ implements the Neumann boundary conditions. We refrain from elaborating further on this here, and simply note that the FEM matrix representation becomes

[TABLE]

Fur further details on the proper abstract framework, see e.g. [21] and Section 8.1.

Since Section 2 is no longer satisfied, we may not apply Theorem 4.6. However, the results plotted in Fig. 5 are similar to the results in previous examples. Again, the LOD approximations are more efficient except for the very coarsest meshes. This indicates that our theory could be extended also to the case of unbounded operators $\mathcal{B}$ and $\mathcal{C}$.

7.5 Example 5

As a final experiment, we consider another boundary control application. The domain is formed like a lying U, see Fig. 6. The thickness of each of the “handles” is $1/6$, the total horizontal extent $1$ and the vertical extent $4/6$. Inside the domain are three evenly spaced stripes with a diameter of $0.0052$. As previously, we consider \mathcal{A}u=\nabla\cdot\big{(}\kappa\nabla u\big{)} where $\kappa=10^{-2}$ everywhere except for in the stripes where instead $\kappa=1$. We use homogeneous Neumann boundary conditions over the whole boundary, except for the two vertical sections on the left. On the top-most vertical part, $\Gamma_{1}$, we impose a nonhomogeneous Neumann condition $\kappa\nabla x\cdot n=\Psi u$ with $\Psi$ having the same hat-shaped form as in Example 4. On the bottom vertical part, $\Gamma_{2}$, we impose a homogeneous Dirichlet condition. These correspond to an insulated edge, a controllable heat input and a heat sink, respectively. The operator $\mathcal{B}$ is again given by $u\mapsto u\int_{\Gamma_{1}}{\Psi\varphi}$, and as output we take the mean of the temperature over the domain; $\mathcal{C}x=\int_{\Omega}x$. The meshes in this example have $n=28,84,280,1008,3808,14784$ interior nodes, respectively, while the reference solution uses $n=58240$.

The results are plotted in Fig. 7, where we can once again observe error behaviour consistent with the bounds given in Theorem 4.6.

Remark 7.1**.**

*In all the experiments, we have chosen the fine-scale structure of the multiscale coefficient such that the reference FEM solution can resolve it, since otherwise we can not properly compute the respective errors. Decreasing the size of the fine-scale features even further would mean that the FEM convergence is further delayed, while we may still compute accurate LOD approximations. In such a case, the efficiency of LOD in comparison to FEM is further (greatly) improved. *

Remark 7.2**.**

*We note that the finest discretizations of the demonstrated numerical experiments are representative of large-scale DRE problems. While there is of course no strict limit, the authors would at the time of writing classify large-scale as problems of size $n\geq 10^{4}$. Due to the matrix-valued nature of the equations, this is naively equivalent to solving a vector-valued differential equation of size $10^{8}$. While we do not employ a naive method, the number of unknowns still number in the millions. For readers more familiar with the theory of algebraic Riccati equations (AREs), i.e. the stationary counterparts of DREs, we note that one time step for a DRE solver is roughly equivalent to the solution of one corresponding ARE. The computational effort for solving a DRE is therefore usually at least two orders of magnitude higher, and large-scale in the ARE setting is therefore larger; starting rather at around $n=10^{5}$. We also note that the approximations were here computed on the equivalent of a modern desktop computer. With the increasing availability of parallelization on clusters or GPUs, we expect to see a shift towards even larger problems in the near future. However, for multiscale problems such as these, it is still critical to employ LOD. *

8 Generalizations and future work

In this section we provide some notes on possible extensions of our theory and draw connections to related problems and methods.

8.1 Boundary control

Boundary control applications such as Example 4 occur frequently within the field of optimal control. Then either the input or output operator (or both) acts on the boundary of the computational domain. In order to put such problems into the semigroup framework, one has to allow for unbounded operators $\mathcal{B}$ and $\mathcal{C}$ [21]. Clearly, our convergence analysis is no longer valid in that case, since we can no longer guarantee that $S_{h}\in\mathcal{L}(L^{2})$ or that $\mathcal{C}_{h}P_{h}\in\mathcal{L}(L^{2},Z)$. However, it is typically assumed that $\mathcal{B}$ and $\mathcal{C}$ are not too unbounded. More specifically, if we suppose that $(-\mathcal{A})^{-\beta}\mathcal{B}\in\mathcal{L}(U,L^{2})$ and $\mathcal{C}(-\mathcal{A})^{-\gamma}\in\mathcal{L}(L^{2},Z)$, where $0\leq\beta+\gamma<1$, we cover a large class of applications. Here, $(-\mathcal{A})^{-\alpha}$, denotes fractional powers of $\mathcal{A}$ which exist due to Section 2. They give rise to the spaces $X_{-\alpha}\supset L^{2}$ as the completions of $L^{2}$ in the norm $\lVert x\rVert_{-\alpha}=\lVert(-\mathcal{A})^{-\alpha}x\rVert$. When $\gamma=0$ we then have that $S_{h}\in\mathcal{L}(X_{-\beta})$, and by properly extending also the other involved operators to $X_{-\beta}$ we may follow the line of proof of Theorem 4.6 and show convergence in $\mathcal{L}(X_{-\beta})$.

Obviously, this is a sub-optimal estimation, as $\lVert\cdot\rVert_{-\beta}$ is a weaker norm than $\lVert\cdot\rVert_{L^{2}}$ for $\beta>0$. However, from [21, Theorem 1.2.1.1] we have that $\tilde{X}(t)S_{h}\tilde{X}(t)$ is actually bounded in $L^{2}$, at least away from $t=0$. It therefore seems likely that one could use similar ideas to prove that the same holds for $\tilde{Y}(t)S_{h}\tilde{Y}(t)$, in which case we would have convergence in $\mathcal{L}(L^{2})$. Unfortunately, the theory required for such estimations is rather extensive, and we expect it to be even more so for the LOD approximations. We therefore leave such questions as future work.

8.2 Systems of equations and applications in multiphysics

In this paper we consider problems where the evolution operator $\mathcal{A}$ in the state equation defines an inner product of the form $a(u,v)=\int\kappa\nabla u\cdot\nabla v$. However, many interesting applications requires coupled systems to be modeled appropriately, for instance, multiphysical features such as thermoelasticity [9], which describes temperature and displacement in a material. Another example is the singularly perturbed systems [17, 29], which appear when modeling, for instance, fluid catalytic crackers. These are ill-conditioned problems due to a significantly larger time derivative for one (or more) of the equations.

The LOD method has successfully been applied to thermoelasticity and poroelasticity problems, see [23]. With more complicated models, the computational gain in using a coarse representation of the underlying partial differential equation is even greater. Analysis of such problems should be considered in the future.

8.3 Other time discretizations

It should also be noted that the LOD approach could be used with other time discretizations of the DRE. We have here chosen the Strang splitting scheme due to it being familiar to one of the authors and because an efficient implementation was readily available. However, there are also other types of splitting schemes [33, 28]. Additionally, one might instead consider e.g. BDF and Rosenbrock methods [6, 5, 20], projection-based methods [16] or even peer methods [19]. These depend on solving linear equation systems rather than computing the solutions to parabolic problems, and the error analysis approach would thus differ. However, bounds similar to that given in Lemma 4.3 naturally exist also for stationary problems [26].

8.4 Algebraic Riccati equations

The latter fact is even more relevant if one considers algebraic Riccati equations (AREs). These are the stationary counterpoints to the time-dependent DREs and arise when the final time $T$ in the cost functional goes to infinity. In this case, splitting does not apply, but we may still apply LOD to the equation to reduce its complexity. Then any method for AREs may be applied to solve this smaller problem, such as Newton-Kleinman ADI [4], rational Krylov subspace methods [31], or RADI [3]. See also [7] for a survey. Clearly, for each of these cases one would have to perform an error analysis such as the one provided in this paper.

Acknowledgments

We are grateful to Fredrik Hellman for his assistance with the code for computing the LOD bases.