Low-rank computation of posterior covariance matrices in Bayesian inverse problems

TL;DR

This paper introduces a low-rank approach to efficiently approximate the posterior covariance matrix in Bayesian inverse problems, significantly reducing computational complexity and storage requirements for large-scale PDE-based models.

Contribution

The authors develop a novel method combining low-rank in time algorithms with low-rank Hessian techniques to efficiently compute posterior covariance matrices in large-scale Bayesian inverse problems.

Findings

Reduces computational complexity from O(n_x n_t) to O(n_x + n_t)

Demonstrates efficiency through numerical experiments

Enables uncertainty quantification in large-scale inverse problems

Abstract

We consider the problem of estimating the uncertainty in statistical inverse problems using Bayesian inference. When the probability density of the noise and the prior are Gaussian, the solution of such a statistical inverse problem is also Gaussian. Therefore, the underlying solution is characterized by the mean and covariance matrix of the posterior probability density. However, the covariance matrix of the posterior probability density is full and large. Hence, the computation of such a matrix is impossible for large dimensional parameter spaces. It is shown that for many ill-posed problems, the Hessian matrix of the data misfit part has low numerical rank and it is therefore possible to perform a low-rank approach to approximate the posterior covariance matrix. For such a low-rank approximation, one needs to solve a forward partial differential equation (PDE) and the adjoint PDE in both space and time. This in turn gives $\mathcal{O}(n_x n_t)$ complexity for both, computation and storage, where $n_x$ is the dimension of the spatial domain and $n_t$ is the dimension of the time domain. Such computations and storage demand are infeasible for large problems. To overcome this obstacle, we develop a new approach that utilizes a recently developed low-rank in time algorithm together with the low-rank Hessian method. We reduce both the computational complexity and storage requirement from $\mathcal{O}(n_x n_t)$ to $\mathcal{O}(n_x + n_t)$. We use numerical experiments to illustrate the advantages of our approach.

Equations

$$Y=g(U,E),$$

$${\pi }_{\mathrm{post}}:=\pi (u|y_{\mathrm{obs}})=\frac{{\pi }_{\mathrm{prior}}(u)\pi (y_{\mathrm{obs}}|u)}{\pi (y_{\mathrm{obs}})}\propto {\pi }_{\mathrm{prior}}(u)\pi (y_{\mathrm{obs}}|u),$$

$$Y=f(U)+E$$

$${\pi }_{\mathrm{noise}}(e)={\pi }_{\mathrm{noise}}(y_{\mathrm{obs}}-f(u)).$$

$${\pi }_{\mathrm{post}}\propto {\pi }_{\mathrm{prior}}(u){\pi }_{\mathrm{noise}}(y_{\mathrm{obs}}-f(u)).$$

$$\begin{array}{rl}{\pi }_{\mathrm{prior}}(u)\propto & \exp \left(-\frac{1}{2}{\left(u-{\bar{u}}_{\mathrm{prior}}\right)}^T{\Gamma }_{\mathrm{prior}}^{-1}\left(u-{\bar{u}}_{\mathrm{prior}}\right)\right),\\ {\pi }_{\mathrm{noise}}(e)\propto & \exp \left(-\frac{1}{2}{\left(e-\bar{e}\right)}^T{\Gamma }_{\mathrm{noise}}^{-1}\left(e-\bar{e}\right)\right),\end{array}$$

$$\begin{array}{rl}{\pi }_{\mathrm{post}}\propto & \exp \left(-\frac{1}{2}{\left(u-{\bar{u}}_{\mathrm{prior}}\right)}^T{\Gamma }_{\mathrm{prior}}^{-1}\left(u-{\bar{u}}_{\mathrm{prior}}\right)-\frac{1}{2}{\left(e-\bar{e}\right)}^T{\Gamma }_{\mathrm{noise}}^{-1}\left(e-\bar{e}\right)\right)\\ =& \exp \left(-\frac{1}{2}{\Vert u-{\bar{u}}_{\mathrm{prior}}\Vert }_{{\Gamma }_{\mathrm{prior}}^{-1}}^2-\frac{1}{2}{\Vert e-\bar{e}\Vert }_{{\Gamma }_{\mathrm{noise}}^{-1}}^2\right).\end{array}$$

$$\begin{array}{rl}{\pi }_{\mathrm{post}}\propto & \exp \left(-\frac{1}{2}{\Vert u-{\bar{u}}_{\mathrm{prior}}\Vert }_{{\Gamma }_{\mathrm{prior}}^{-1}}^2-\frac{1}{2}{\Vert y_{\mathrm{obs}}-Au-\bar{e}\Vert }_{{\Gamma }_{\mathrm{noise}}^{-1}}^2\right).\end{array}$$

$${\bar{u}}_{\mathrm{post}}={\mathrm{argmax}}_u{\pi }_{\mathrm{post}}(u)$$

$${\bar{u}}_{\mathrm{post}}={\mathrm{argmin}}_u\left(\frac{1}{2}{\Vert u-{\bar{u}}_{\mathrm{prior}}\Vert }_{{\Gamma }_{\mathrm{prior}}^{-1}}^2+\frac{1}{2}{\Vert y_{\mathrm{obs}}-Au-\bar{e}\Vert }_{{\Gamma }_{\mathrm{noise}}^{-1}}^2\right).$$

$${\Gamma }_{\mathrm{post}}={\left(A^T{\Gamma }_{\mathrm{noise}}^{-1}A+{\Gamma }_{\mathrm{prior}}^{-1}\right)}^{-1}.$$

$$\begin{array}{rl}{\Gamma }_{\mathrm{post}}& ={\left(A^T{\Gamma }_{\mathrm{noise}}^{-1}A+{\Gamma }_{\mathrm{prior}}^{-1}\right)}^{-1}\\ & ={\Gamma }_{\mathrm{prior}}^{1/2}{\left({\Gamma }_{\mathrm{prior}}^{1/2}{A}^T{\Gamma }_{\mathrm{noise}}^{-1}A{\Gamma }_{\mathrm{prior}}^{1/2}+I\right)}^{-1}{\Gamma }_{\mathrm{prior}}^{1/2}.\end{array}$$

$${\tilde{\mathcal{H}}}_{\text{mis}}={\Gamma }_{\mathrm{prior}}^{1/2}{A}^T{\Gamma }_{\mathrm{noise}}^{-1}A{\Gamma }_{\mathrm{prior}}^{1/2}$$

$${\tilde{\mathcal{H}}}_{\text{mis}}\approx V\Lambda V^T,$$

$${\left({\Gamma }_{\mathrm{prior}}^{1/2}{A}^T{\Gamma }_{\mathrm{noise}}^{-1}A{\Gamma }_{\mathrm{prior}}^{1/2}+I\right)}^{-1}\approx {\left(V\Lambda V^T+I\right)}^{-1}=I-V\tilde{\Lambda }{V}^T$$

$$\begin{array}{rl}{\left(A^T{\Gamma }_{\mathrm{noise}}^{-1}A+{\Gamma }_{\mathrm{prior}}^{-1}\right)}^{-1}& \approx {\left(V\Lambda V^T+{\Gamma }_{\mathrm{prior}}\right)}^{-1}\\ & ={\Gamma }_{\mathrm{prior}}-{\Gamma }_{\mathrm{prior}}V{\left({\Lambda }^{-1}+V{\Gamma }_{\mathrm{prior}}{V}^T\right)}^{-1}\tilde{\Lambda }{V}^T{\Gamma }_{\mathrm{prior}},\end{array}$$

$$\underset{y,u}{\min }\frac{1}{2}{\Vert y-y_{\mathrm{obs}}\Vert }_{\mathcal{Q}}^2+\frac{\beta }{2}{\Vert u\Vert }_{\mathcal{P}}^2,$$

$$\begin{array}{rl}{y}_t-{\nabla }^{2}y=& u,\text{in }\Omega ,\\ y=& f,\text{on }\partial \Omega .\end{array}$$

$$\begin{split}\underbrace{\left[\begin{array}[]{ccc}{D}_{1}\otimes\tau{M}_{1}&0&-\left({I}_{n_{t}}\otimes{L}+{C}^{T}\otimes{M}\right)\\ 0&{D}_{2}\otimes\beta\tau{M}_{2}&{D}_{3}\otimes\tau{N}^{T}\\ -\left({I}_{n_{t}}\otimes{L}+{C}\otimes{M}\right)&{D}_{3}\otimes\tau{N}&0\\ \end{array}\right]}_{\mathcal{A}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left[\begin{array}[]{c}\mathbf{{y}}\\ \mathbf{{u}}\\ \mathbf{{p}}\\ \end{array}\right]}&\\ =\left[\begin{array}[]{c}{D}_{1}\otimes\tau{M}_{1}\mathbf{y_{\text{obs}}}\\ \mathbf{0}\\ \mathbf{d}\\ \end{array}\right]&,\end{split}$$

$$\mathbf{{y}}=\left[\begin{array}[]{c}\mathbf{{y}}_{1}\\ \vdots\\ \mathbf{{y}}_{n_{t}}\\ \end{array}\right],\mathbf{{u}}=\left[\begin{array}[]{c}\mathbf{{u}}_{1}\\ \vdots\\ \mathbf{{u}}_{n_{t}}\\ \end{array}\right],\textnormal{ and }\mathbf{{p}}=\left[\begin{array}[]{c}\mathbf{{p}}_{1}\\ \vdots\\ \mathbf{{p}}_{n_{t}}\\ \end{array}\right].$$

$$W\otimes V=\left[\begin{array}[]{ccc}w_{11}V&\ldots&w_{1m}V\\ \vdots&\ddots&\vdots\\ w_{n1}V&\ldots&w_{nm}V\\ \end{array}\right]$$

$$\mathop{\mathrm{vec}}\left(W\right)=\left[\begin{array}[]{c}w_{11}\\ \vdots\\ w_{n1}\\ \vdots\\ w_{nm}\\ \end{array}\right]$$

$$\left(W^T\otimes V\right)\mathrm{vec}\left(Y\right)=\mathrm{vec}\left(VYW\right).$$

$$\begin{array}{rl}Y=& \left[\begin{array}{c}{\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_{n_t}\end{array}\right]\approx W_Y{V}_Y^T\text{ with }{W}_Y\in {\mathbb{R}}^{n_1,k_1},V_Y\in {\mathbb{R}}^{n_t,k_1},\\ U=& \left[\begin{array}{c}{\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_{n_t}\end{array}\right]\approx W_U{V}_U^T\text{ with }{W}_U\in {\mathbb{R}}^{n_2,k_2},V_U\in {\mathbb{R}}^{n_t,k_2},\\ P=& \left[\begin{array}{c}{\mathbf{p}}_1,{\mathbf{p}}_2,\dots ,{\mathbf{p}}_{n_t}\end{array}\right]\approx W_P{V}_P^T\text{ with }{W}_P\in {\mathbb{R}}^{n_1,k_3},V_P\in {\mathbb{R}}^{n_t,k_3},\end{array}$$

$${\Gamma }_{\mathrm{prior}}={\left(D_2\otimes \beta \tau M_2\right)}^{-1}\text{ and }{\Gamma }_{\mathrm{noise}}={\left(D_1\otimes \tau M_1\right)}^{-1},$$

$$\begin{array}{rl}{\Gamma }_{\mathrm{post}}=& [\left(D_2\otimes \beta \tau M_2\right)+\left(D_3\otimes \tau N^T\right){\left(I_{n_t}\otimes L+C^T\otimes M\right)}^{-1}\\ & {\left(D_1\otimes \tau M_1\right){\left(I_{n_t}\otimes L+C\otimes M\right)}^{-1}\left(D_3\otimes \tau N\right)]}^{-1}\end{array}$$

$$\begin{array}{rl}{\tilde{\mathcal{H}}}_{\text{mis}}=& [{\left(D_2\otimes \beta \tau M_2\right)}^{1/2}\left(D_3\otimes \tau N^T\right){\left(I_{n_t}\otimes L+C^T\otimes M\right)}^{-1}\\ & \left(D_1\otimes \tau M_1\right){\left(I_{n_t}\otimes L+C\otimes M\right)}^{-1}\left(D_3\otimes \tau N\right){\left(D_2\otimes \beta \tau M_2\right)}^{1/2}].\end{array}$$

$$\begin{array}{r}\underset{y,u}{\min }\frac{{\beta }_{\mathrm{noise}}}{2}{\Vert y-y_{\mathrm{obs}}\Vert }_{\mathcal{Q}}^2+\frac{{\beta }_{\mathrm{prior}}}{2}{\Vert u\Vert }_{\mathcal{P}}^2,\end{array}$$

$${\Gamma }_{\mathrm{prior}}={\left(D_2\otimes {\beta }_{\mathrm{prior}}\tau M_2\right)}^{-1}\text{ and }{\Gamma }_{\mathrm{noise}}={\left(D_1\otimes \tau {\beta }_{\mathrm{noise}}{M}_1\right)}^{-1}.$$

$${\Gamma }_{\mathrm{prior}}^{-1}={\beta }_{\mathrm{prior}}\tau h^{d}I\text{ and }{\Gamma }_{\mathrm{noise}}^{-1}={\beta }_{\mathrm{noise}}\tau h^{d}I.$$

Full text

{AMS}

65F15 , 65F10, 65F50, 93C20, 62F15

1 Introduction

Computational mathematicians dealing with simulations of large-scale discretizations describing physical phenomena have made tremendous success over the last decades. This has enabled scientists from various areas of engineering, chemistry, geophysics, et al. to ask more relevant and complex questions. One area that has seen a dramatic increase in the number of published results is the field of statistical inverse problems [36, 20, 8]. In particular, the consideration of partial differential equations (PDEs) as models in statistical inverse problems dramatically increases the problem complexity as a refinement of the model in space and time results in an exponential increase in the problem degrees of freedom. By this we mean that a discretized problem is typically represented by a spatial system matrix $A\in\mathbb{R}^{n,n}$ where the number of degrees of freedom ${n_{x}}$ is typically $O(\frac{1}{h^{d}})$ with $d$ being the spatial dimension, and $h$ is the mesh size. It is easily seen that halving the parameter $h$ means the matrix size will grow by a factor of $2,4,8,\ldots$ depending on the spatial dimension. This complexity is further increased when the temporal dimension is incorporated.

While numerical analysis has provided many techniques that allow the efficient handling of such problems, e.g. Krylov methods [28], multigrid techniques [15], we are faced with an even steeper challenge when uncertainty in the parameters of the model is incorporated. For this we consider the approach of Bayesian inverse problems where the goal is to use prior knowledge to obtain information about the conditional mean or the posterior covariance given a set of measured data. While computing the conditional mean typically corresponds to a problem formulation frequently encountered in PDE-constrained optimization [37, 17, 6], the problem of computing the posterior covariance matrix is much more challenging as this matrix is dense and involves the inverse of high-dimensional discretized PDE problems. In [12] and subsequent works, the authors proposed a low-rank approximation of the posterior covariance matrix. While this already reduces the complexity dramatically, the storage requirements for the resulting vectors still suffer from high dimensionality with respect to $n_{x}$ and $n_{t}$. Our aim in this paper is to introduce a method based on low-rank techniques [35] that allows to reduce the complexity from $O({n_{x}}{n_{t}})$ to $O({n_{x}}+{n_{t}})$.

In order to do this, we will first derive the basic problem following [12]. This will be followed by the presentation of a low-rank technique that we previously introduced for PDE-constrained optimization. We then use this to establish a low-rank eigenvalue method based on the classical Lanczos procedure [24] or Arnoldi procedure [29]. After introducing different choices of covariance matrices, we show that our approach can be theoretically justified. We then illustrate the applicability of our proposed methodology to a diffusion problem and a convection diffusion problem, and present numerical results illustrating the performance of the scheme.

2 Statistical/Bayesian Inverse Problems

We refer to [20, 36] for excellent introductions into the subject of statistical inverse problems. We follow [12, 7] in the derivation of our model setup and start with the parameter-to-observable map $g:\mathbb{R}^{n}\times\mathbb{R}^{k}\rightarrow\mathbb{R}^{m}$ defined as

[TABLE]

where $U,Y,E$ are vectors of random variables. Note that here, $u\in\mathbb{R}^{n}$, our vector of model parameters to be recovered, is a realization of $U,$ the error vector $e\in\mathbb{R}^{k}$ is a realization of $E,$ and $y\in\mathbb{R}^{m}$ is a realization of the vector of observables $Y$. The vector ${{y}_{\operatorname{obs}}}\in\mathbb{R}^{m}$ contains the observed values. As discussed in [7], even when using the ‘true’ model parameters $u,$ the observables ${y}$ will differ from the measurements ${{y}_{\operatorname{obs}}}$ due to measurement noise and the inadequacy of the underlying PDE model.

In a typical application such as the one discussed later, evaluating $g(.)$ requires the solution of a PDE potentially coupled to an observation operator representing a domain of interest.

The Bayes’ theorem, which plays a key role in the Bayesian inference, is written as

[TABLE]

where we used the prior probability density function (PDF) $\pi_{\operatorname{prior}}(x)$, the likelihood function $\pi({{y}_{\operatorname{obs}}}|u)$, and the data ${{y}_{\operatorname{obs}}}$ with $\pi({{y}_{\operatorname{obs}}})>0$. The function $\pi_{\operatorname{post}}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is the posterior probability density function. The likelihood is derived under the assumption of additive noise

[TABLE]

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ and $E$ is the additive noise and given as $E=Y-f(U).$ We once more follow [12], assuming that $U$ and $E$ are statistically independent and we can use

[TABLE]

Therefore, Bayes’ theorem can be written as

[TABLE]

Assuming that both probability density functions for $U$ and $E$ are Gaussian, we can rewrite the PDFs in the form

[TABLE]

where $\bar{u}_{\operatorname{prior}}\in\mathbb{R}^{n}$ is the mean of the model parameter prior PDF and $\bar{e}$ is the mean of the noise PDF. We further have the two covariance matrices $\Gamma_{\operatorname{prior}}\in\mathbb{R}^{n,n}$ for the prior and $\Gamma_{\operatorname{noise}}\in\mathbb{R}^{m,m}$ for the noise. The Gaussian assumption allows us to rewrite Bayes’ theorem further to get

[TABLE]

Let us further assume that the parameter-to-observable map $g(U,Y)$ is given as in (3) with $f(U)=AU$. The matrix $A\in\mathbb{R}^{m,n}$ represents a linear map from the parameters $u$ to the observables $y$. We will later see that often this matrix involves the inverse of a discretized representation of a PDE operator. In [33], the authors incorporate a quantity of interest operator into $A$. Therefore, it will typically be dense and very large. We arrive now at a restated version of Bayes theorem (3)

[TABLE]

From this relation we can express several relevant statistical quantities. For example, we can compute the maximum a posteriori point (MAP), which is defined via

[TABLE]

and to compute it, one can solve the following optimization problem

[TABLE]

Note that this problem is a deterministic inverse problem and resembles the structure one finds in PDE-constrained optimization problems [37, 19]. Many strategies are known how to solve this efficiently and in particular underlying function space considerations help in the development of efficient numerical algorithms. An infinite-dimensional discussion of the above problem is given in [7, 36] and we only refer to the infinite-dimensional setup when needed. Our goal in this paper will not be the solution of the MAP problem. The goal of devising low-rank methods for this case has recently been established in [35] and the techniques there are likely to be applicable as the only difference is the use of the weighting matrices $\Gamma_{\operatorname{noise}}\textrm{ and }{\Gamma_{\operatorname{prior}}},$ which are for the classical PDE-constrained optimization problem mass matrices or matrices involving mass matrices. The more challenging question lies in the approximation of the posterior covariance matrix

[TABLE]

The approximation of $\Gamma_{\operatorname{post}}$ is in general very costly and, without further approximation, intractable. The approach presented in [12, 7] computes a low-rank approximation to this matrix using the following relation

[TABLE]

The authors in [12, 7] then compute a low-rank approximation to the so-called prior-preconditioned Hessian of the data misfit $\tilde{\mathcal{H}}_{\text{mis}}\in\mathbb{R}^{n,n}$

[TABLE]

with the approximation

[TABLE]

where $V$ and $\Lambda$ represent the dominant eigenvectors and eigenvalues, respectively. Using this approximation and the Sherman-Morrison-Woodbury formula [13] one obtains for the prior-preconditioned system

[TABLE]

with $\tilde{\Lambda}=\mathrm{diag}(\frac{\lambda_{i}}{\lambda_{i}+1})$, and $\lambda_{i}$ is the $i$-th diagonal entry of $\Lambda$. A similar approach for

[TABLE]

becomes a feasible alternative if $\left(\Lambda^{-1}+V\Gamma_{\operatorname{prior}}V^{T}\right)^{-1}$ can be evaluated in reasonable time.

Our main goal of this paper will be the derivation of an efficient scheme to approximate the matrix $V$ from the low-rank approximation to the misfit Hessian. Before discussing this problem, we introduce an idea that becomes instrumental in realizing this and is motivated by a PDE-constrained optimization problem.

3 A Low-Rank Technique for PDE-Constrained Optimization

In order to better understand the stochastic inverse problem, we investigate it in relation to a PDE-constrained optimization problem. We start the derivation of the low-rank in time method by considering an often used model problem in PDE-constrained optimization (see [18, 19, 37]), minimization of

[TABLE]

with $\mathcal{P}$ and $\mathcal{Q}$ space-time cylinders. The constraint linking state ${y}$ and control ${u}$ of this problem is given by the heat equation with a distributed control term

[TABLE]

Here $\Omega$ denotes the domain and $\partial\Omega$ corresponds to the boundary of the domain. For a more detailed discussion on the well-posedness, existence of solutions, discretization et al., we refer the interested reader to [18, 19, 37]. The solution of such an optimization problem is obtained using a Lagrangian approach and considering the first order conditions, which for our problem results in a linear system of the form

[TABLE]

where ${D}_{1}={D}_{2}=\mathop{\mathrm{diag}}\left(1,1,\ldots,1,1\right)$ and ${D}_{3}={I}_{n_{t}}$ come from the discretization of the temporal parts of the objective function or the right hand side of the PDE-constraint (cf. [5, 27, 34]). The matrices $M_{1}$ and $M_{2}$ are mass matrices corresponding to observation and control domain. The matrix $N$ is essentially representing the incorporation of the control into the constraint, i.e., $N$ is a mass matrix in the above example. The matrix $C$ represents the all-at-once discretization of the time-derivative in the PDE and $L$ the discretized Laplacian. Here, the state, control, and adjoint state are represented by the following space-time vectors

[TABLE]

We point out again that the Kronecker product is defined as

[TABLE]

and remind the reader of the definition of the $\mathop{\mathrm{vec}}\left(\right)$ operator via

[TABLE]

as well as the relation

[TABLE]

In [35], it was shown that the solution to the PDE-constrained optimization problem can be computed in low-rank form

[TABLE]

where the $k_{i}$ are small in comparison to the spatial and temporal dimensions. The authors in [35] illustrated that the low-rank structure of a right hand side is maintained throughout a Krylov subspace iteration and the above described representation. Low-rank techniques for Krylov-subspace methods have recently received much attention and we refer the reader to [21, 22, 1] and for tensor structured equations [14, 26, 9, 10].

We obtain a significant storage reduction if we can base our approximation of the solution using such low-rank factors. It is easily seen that due to the low-rank nature of the factors, we have to perform fewer multiplications with the submatrices by also maintaining smaller storage requirements.

There are several similarities of the problem (14) and the statistical inverse problem presented earlier. It is clear that with the choice

[TABLE]

the PDE-constrained optimization problem can be interpreted as a statistical inverse problem and the posterior covariance matrix $\Gamma_{\operatorname{post}}$ is given by eliminating both state and adjoint state from the system matrix (14) to obtain a reduced Hessian system. Furthermore, it is clear that for many choices of prior and noise covariance, we can utilize the tensor structure to compute low-rank solutions. For this we state the posterior covariance matrix of the PDE optimization problem

[TABLE]

and the misfit Hessian

[TABLE]

We keep this example in mind when we now discuss a low-rank technique to approximate the eigenvectors of the posterior covariance matrix in low-rank form. For this we propose a low-rank Krylov subspace method to compute the dominating eigenvectors and eigenvalues.

Before discussing the eigenvalue approximation strategy, we want to comment on the scaling of the PDE-constrained optimization problem in relation to the statistical inverse problem discussed in [12]. The authors there consider

[TABLE]

which in the simple case of full observations and control leads to the following rescaling of (16)

[TABLE]

Assuming that $M\approx h^{d}I$ and $D_{i}=I$ we get

[TABLE]

In PDE-constrained optimization one typically reduces $\beta$ in (13) to allow for a more expensive control that drives the state closer to the desired state. This would mean that in the statistical inverse setting $\beta_{\operatorname{noise}}=1$ and decreasing the value of $\beta_{\operatorname{prior}},$ which implies that in (9) the role of the prior covariance gets diminished and most contributions are coming from the noise. We have similar settings and observations for stochastic inverse problems in this manuscript. These will be shown by the analysis in Section 5 and numerical experiments in Section 6. The right choice of parameters $\beta_{\text{prior}}$ and $\beta_{\text{noise}}$ depends on the underlying application and we refer to [8] for a discussion of the roles of the parameters $\beta_{\operatorname{prior}}$ and $\beta_{\operatorname{noise}}$ as regularization parameters.

4 Low-rank Lanczos/Arnoldi Method

We recall that our goal is to find a low-rank approximation of the eigenvectors of the posterior covariance matrix. The goal is to compute an approximation to $\tilde{\mathcal{H}}_{\text{mis}}\approx V\Lambda V^{T}$ with $V=[v_{1},v_{2},\ldots,v_{k}]$ and $k$ much smaller than the dimension of $\tilde{\mathcal{H}}_{\text{mis}}.$ For the PDE-constrained optimization problem $\tilde{\mathcal{H}}_{\text{mis}}\in\mathbb{R}^{{n_{x}}{n_{t}},{n_{x}}{n_{t}}}$.

Our main assumption at this point is that storing each $v_{j}$ and especially a number of such space-time vectors can pose serious problems. Additionally, in order to perform the matrix vector multiplication with $\tilde{\mathcal{H}}_{\text{mis}}$, a large number of PDE-solutions need to be computed. For this we point out that in order to apply the matrix $\tilde{\mathcal{H}}_{\text{mis}}$ in an Arnoldi procedure, we need to solve the spatial system over the whole time-domain. A major advantage of our approach is motived by the fact that

[TABLE]

which we assume is well approximated via

[TABLE]

with $W_{j,1}\in\mathbb{R}^{{n_{x}},r_{j}}$ $W_{j,2}\in\mathbb{R}^{{n_{t}},r_{j}}$ with $r_{j}\ll\min\left\{{n_{x}},{n_{t}}\right\}.$ If the Arnoldi- or Lanczos vectors are of this form, then the application of the matrix $\tilde{\mathcal{H}}_{\text{mis}}$ to such vectors requires fewer PDE solves than in the full case.

Note that we need to compute the dominant eigenvectors of the prior-preconditioned data misfit Hessian $\tilde{\mathcal{H}}_{\text{mis}}$ (11), therefore the Lanczos method can be used. At the $j$-th Lanczos iteration, we need to perform $\tilde{\mathcal{H}}_{\text{mis}}v_{j-1}$ using the low-rank approach to get a form like  (22). Here $v_{j-1}$ is the $(j-1)$-th Lanczos vector. Due to the low-rank approximation, the orthogonality of Lanczos vectors is lost. Reorthogonalization should be used to orthogonalize Lanczos vectors. Meanwhile, the symmetric property of $\tilde{\mathcal{H}}_{\text{mis}}$ cannot be preserved for the low-rank form of the matrix-vector product. Therefore, we make use of the more general Arnoldi method to compute the dominant eigenvectors of $\tilde{\mathcal{H}}_{\text{mis}}$. We also observe that when applying the Arnoldi method with the truncation error appropriately chosen we still get real eigenvalues. This will be shown in Section 6.

We now briefly recall the Arnoldi method, which is the more general procedure. We refer to [13] for details. We recall that the Arnoldi process for a matrix $B$ can be written as

[TABLE]

where $V_{k}$ consists of orthonormal columns and $H_{k+1,k}\in\mathbb{R}^{(k+1)\times k}$ is a Hessenberg matrix. The iterative build-up of the columns of $V$ is captured by the recursion

[TABLE]

where $h_{i,k}=v_{i}^{T}Av_{k}$. The vector $\tilde{v}_{k+1}$ is then normalized using the scalar $h_{k+1,k}.$ While this is well-known our goal here is to illustrate how this method is amenable to the use within a low-rank framework. For the Arnoldi process considered in this manuscript, $B=\tilde{\mathcal{H}}_{\text{mis}}$ and the application of $B$ to $v_{k}$ results in a low-rank matrix, i.e.,

[TABLE]

with small rank. This is because $B$ is related with the inverse of a PDE operator in space and time, which has shown in [35] that applying such an operator to a low-rank vector again gives a low-rank vector. We can then write the right hand side of (23) as

[TABLE]

and write the last expression as

[TABLE]

The size of the matrix

[TABLE]

is increased to $r_{B}+r_{k}+r_{k-1}.$ Using truncation techniques this can typically be controlled. For example, one could achieve the truncation by utilizing skinny QR factorizations [21] or truncated singular value decompositions [35].

In our numerical experiments, we show that maintaining the orthogonality of the Lanczos vectors is crucial. Meanwhile, since we use low-rank methods to approximate the Lanzcos vectors, the orthogonality and the symmetry of the Lanzcos process get lost. We hence opt for the more expensive but stable Arnoldi process, where we orthogonalize with respect to all previous vectors. The full reorthogalization also demands more storage than in the Lanczos case. This is another advantage of our approach. With full reorthogonalization, the storage cost are increasing for both full and low-rank scheme but in the low-rank framework stay significantly below the full scheme. Here, we give the low-rank Arnoldi method in Algorithm 1.

We note that the biggest challenge for the low-rank Arnoldi method is to perform the low-rank matrix vector product in line $3$ of Algorithm 1 since $\tilde{\mathcal{H}}_{\text{mis}}$ is large and dense. We propose the tensor-train (TT) format in Section 6 to perform such computations efficiently. The full orthogonalization procedure in line $4$-$7$ of Algorithm 1 is also performed with the TT format.

We use the standard Arnoldi method for low-rank eigenvector computations in Algorithm 1. This is practical for the problems studied in this manuscript since we just need to compute up to a few hundred Arnoldi vectors. Since the computational complexity of full orthogonalization increases with the number of Arnoldi vectors, if more Arnoldi vectors are needed, the restarted Arnoldi method can be implemented with a low-rank version [13].

5 Analysis of the Eigenfunctions

The eigenfunction analysis for the general case presented above is not straightforward. Our goal in this section is to give a theoretical justification for simple cases. We start with the case of a steady state problem involving the two-dimensional Poisson equation. For this we consider the misfit Hessian

[TABLE]

Assuming for now that $\Gamma_{\operatorname{prior}}=\beta_{\operatorname{prior}}I_{d}$ and $\Gamma_{\operatorname{noise}}=\beta_{\operatorname{noise}}I_{d}$ so that we are left with

[TABLE]

Note that we have assumed for the data misfit Hessian to be defined in function space with the prior and noise covariance operators to be multiples of the identity operator. Our goal is to find eigenfunctions of $A^{T}A$ with $A$ the inverse of the PDE operator and $A^{T}$ the inverse of the adjoint PDE operator. For this we use the eigenfunctions of the two-dimensional Laplacian operator satisfying

[TABLE]

with zero Dirichlet boundary conditions (see [23]). Note that for this problem $A^{T}=A$. We here consider the eigenvalue problem on a rectangular domain $[0,a]\times[0,b]$ and state the eigenfunctions as

[TABLE]

with $x_{1}\in[0,a]$ and $x_{2}\in[0,b]$. The associated eigenvalues are given by

[TABLE]

with $m,n=1,2,\ldots$ (see [23]). Coming back to the misfit Hessian we can write this as

[TABLE]

with $A$ the inverse Laplacian. Assuming zero Dirichlet boundary conditions it holds that

[TABLE]

is satisfied for the eigenfunctions of the Laplacian

[TABLE]

with the eigenvalues given by

[TABLE]

where the decay of $\lambda_{m,n}^{-2}$ is quite rapid. This justifies the approximation of $\tilde{\mathcal{H}}_{\text{mis}}$ by a small number of eigenfunctions. In [12], the authors also justify this for a one-dimensional convection diffusion equation.

The goal of our approach is to add even more low-rankness to the computation. We want to illustrate that the eigenfunctions used to represent the misfit Hessian are of small rank. When studying the eigenfunctions above it is clear that the eigenfunction

[TABLE]

is already separated into $x_{1}$ and $x_{2}$ components with separation rank $1$ (cf. [16]), where it is the sum of two functions with one depending on the first and the other on the second variable. For this case we have established the following lemma

Lemma 5.1**.**

The eigenfunctions of the misfit Hessian

[TABLE]

with $A$ the inverse Laplacian with zero Dirichlet conditions defined on the rectangle $[0,a]\times[0,b]$, are separated and given by

[TABLE]

*and hence are of separation rank $1$. *

It is not so straightforward to establish similar results for more complicated equations.

For the space-time PDE-constrained optimization problem discussed earlier, we note that

[TABLE]

where we used ${M}\approx h^{d}{I}$ and collected all scalars in $\alpha(h,\tau,\beta)$. Our aim is to establish eigenvalue and eigenvector results for

[TABLE]

We note that this fits the well known relation that the singular values of a matrix $\mathcal{A}\in\mathbb{R}^{m,m}$ are the square roots of the eigenvalues of the matrix $\mathcal{A}^{T}\mathcal{A},$ which, assuming full rank of $\mathcal{A}$, is a symmetric and positive definite matrix. Now assuming the SVDs

[TABLE]

we obtain

[TABLE]

From this it follows that

[TABLE]

is the eigendecomposition of $\mathcal{A}^{T}\mathcal{A},$ which has the same eigenvectors as $\mathcal{A}^{-1}\mathcal{A}^{-T}.$ As the eigenvalues of $\mathcal{A}^{-1}\mathcal{A}^{-T}$ quickly decay to zero because of the compactness of the operator, we only need a small number of columns of $V$. Our aim in this paper is to express each column of $V$ further in a low-rank fashion. For this we note that $e_{1}^{({n_{x}}{n_{t}})}=e_{1}^{({n_{x}})}\otimes e_{1}^{({n_{t}})}$ and ignoring superindices we get

[TABLE]

and hence a vector of rank one if the eigenvectors are all real. Complex eigenvectors would further introduce a small rank increase and the consideration of ${M}$ instead of $h^{d}{I}$ can with a simultaneous diagonalization of the pencil $({L},{M})$ lead to small eigenvector ranks of the overall system. This justifies our choice of approximating the eigenvectors in low-rank form.

6 Numerical Results

In this section, we study the performance of the low-rank Lanczos algorithm presented in Section 4. The results presented in this section are based on an implementation of the above described algorithms within MATLAB®. We perform the discretization of the PDE-operators within the IFISS [30] framework using $Q1$ finite elements for the heat equation and the streamline upwind Petrov–Galerkin (SUPG) method for the convection diffusion equation. Our experiments are performed for a final time $T=1$ with a varying number of time-steps. As the domain $\Omega$ we consider the unit square but other domains are of course possible. We specify the boundary conditions for each problem separately. Throughout the results section we fixed the truncation at $10^{-8}$ for which we observed good results. Additionally, we also performed not listed experiments with a tolerance of $10^{-10}$ for which we also observed fast convergence. Larger tolerances should be combined with a deeper analysis of the algorithms and a combination with flexible outer solvers. All results are performed on a standard Ubuntu desktop Linux machine with Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz and 8GB of RAM.

The mathematical model we consider in this section is given by the following instationary PDE,

[TABLE]

where $\mathcal{L}$ is a PDE operator and for the numerical experiments we consider the case of the heat equation $\mathcal{L}=-\Delta$ and the convection-diffusion equation $\mathcal{L}=-\mu\Delta+\overrightarrow{\omega}\cdot\triangledown$. For all the numerical tests, the initial concentration $u$ represents the unknown parameter $\mathbf{u}$, and the observation data $\mathbf{y}_{\text{obs}}$ are collected by sensors, which are distributed in part of the domain $\Omega$.

As stated in Section 2, the statistical inverse problem with Gaussian noise and prior using a Bayesian formulation is related to a weighted least squares problem. Here we use the same functional as in [12], which is given by the following functional

[TABLE]

in which $u$ satisfies the PDE model Eq. 25, $b(x,t)$ is the observation operator, and $u$ is the uncertainty term, which is the initial condition for this numerical example. Here, we study the sparse observation case, where $b(x,t)$ is defined by

[TABLE]

Here $\delta_{j}$ is the regional function of sensors, $\delta_{j}=1$ at the region of the $j$-th sensor and $\delta_{j}=0$ elsewhere.

Discretizing the functional Eq. 26 in turn gives

[TABLE]

where $\mathcal{B}$ is the discretization of $b(x,t)$, the variable ${\bf y}$ is the discrete variant of $y$ in Eq. 25 stacked in time, and it satisfies $\mathcal{K}{\bf y}=\mathcal{C}u$. Here $\mathcal{K}$ and $\mathcal{C}$ come from the discretization of the PDE model Eq. 25.

Discretization of the objective function gives that $\Gamma_{\text{noise}}=\nicefrac{{1}}{{\beta_{\text{noise}}}}I_{n_{t}}\otimes M$, and $\Gamma_{\text{prior}}=\nicefrac{{1}}{{\beta_{\text{prior}}}}M$, where $M$ is the mass matrix, $I_{n_{t}}$ is an $n_{t}\times n_{t}$ identity matrix. Here, $n_{t}$ is the number of time variables. Since we need to compute $\Gamma_{\text{prior}}^{\nicefrac{{1}}{{2}}}$, we use $\Gamma_{\text{prior}}\approx\nicefrac{{1}}{{\beta_{\text{prior}}}}h^{-d}I_{n_{x}}$. Here, $h$ is the mesh size, $d$ is the spatial dimension, $I_{n_{x}}$ is a $n_{x}\times n_{x}$ identity matrix, and $n_{x}$ is the number of spatial variables. Other settings such as smoothing operators used for $\Gamma_{\text{prior}}$ have been studied in [7] and can be used with our approach as well.

As analyzed in Section 2, the posterior covariance matrix $\Gamma_{\text{post}}$ is given by the inverse of the Hessian of Eq. 27. Therefore,

[TABLE]

Note that for different PDE models, $\mathcal{K}$ is also different after discretization. In this section, we use two types of PDE models, i.e., the heat equation and the convection-diffusion equation to show the performance of our low-rank algorithm for the approximation of $\Gamma_{\text{post}}$. We argue that our low-rank algorithm also applies to other time-dependent PDE operators. Here we apply our method to symmetric systems and unsymmetric systems.

We also point out that due to the uncertainty being given as the initial condition, the posterior Hessian is only of spatial dimension. Nevertheless, the Arnoldi process applied to this matrix requires the solution of space-time problems and the low-rank form of our approach results in a much reduced number of spatial solves. The complexity reduction is even more pronounced when the uncertainty is part of the system as a space-time variable such as the right hand side of the PDE.

6.1 Implementation Details

According to (28), the prior-preconditioned data misfit part after discretization of (26) is given by

[TABLE]

To apply the Lanczos iteration to (29), we need to solve the space-time discretized PDE $\mathcal{K}$ and adjoint PDE $\mathcal{K}^{T}$. Here we take $\mathcal{K}$ as an example. This asks us to solve a linear system of the following type,

[TABLE]

or

[TABLE]

Here $X$ and $F$ are matrices/tensors of appropriate sizes. Numerical solutions of (30) have been studied intensively from the matrix equation point of view, c.f. [31, 2, 4] for an overview.

In this manuscript, we solve (30) based on its tensor structure and use the alternative minimal energy (AMEn) approach [10] implemented in the tensor-train (TT) toolbox [25] to solve the tensor equation (30). At each AMEn iteration, either a left Galerkin projection or a right Galerkin projection is applied to the system (30). Therefore, after Galerkin projection, we need to solve a linear system either of the format

[TABLE]

or

[TABLE]

Here $\hat{I}$, $\hat{C}$, $\tilde{L}$, $\tilde{M}$ are matrices of appropriate dimensions after Galerkin projection (c.f. [10] for details).

After Galerkin projection, the size of the system (31) is still relatively large while the size of (32) is quite moderate. Therefore, Krylov solvers such as the generalized minimized residual (GMRES) [28] method or the induced dimension reduction IDR(s) [32] method can be used to solve (31) while a direct method can be used to solve (32).

To accelerate the convergence of the Krylov solver, we use the following preconditioner

[TABLE]

to solve (31). Here $\text{diag}(\cdot)$ is an operation that extracts the diagonal entries of a matrix and form a diagonal matrix. One can use standard techniques such as multigrid methods [38] or incomplete LU factorization (ilu) [13] to approximate the preconditioner (33). Here we use backslash implemented in MATLAB®.

We also want to point out that there are many other methods to efficiently solve (30), such as the low-rank factored alternating directions implicit (ADI) method (cf. [3]).

6.2 The Heat Equation

In this part, we use the 2D time-dependent heat equation in a unit square as an example to study the performance of our low-rank algorithm for the heat equation. Discretizing the equation in space using $Q_{1}$ finite elements and in time using the implicit Euler method gives us an $n_{x}\times n_{t}$ linear system, where $n_{x}$ is the number of spatial variables while $n_{t}$ is the number of time steps. First we study spectral properties of the prior-preconditioned data misfit part $\tilde{\mathcal{H}}_{\text{mis}}$ and the posterior covariance matrix. Using a $64\times 64$ grid to discretize the heat equation and set $n_{t}=30,\ 60,\ 90$, respectively, we plot the $50$ largest eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$ in Fig. 1. Here $\beta_{\text{noise}}=10^{4}\beta_{\text{prior}}$.

As shown in Fig. 1, there are only a few dominant eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$. For most cases, a threshold of $10^{-1}$, or even $10^{0}$ is acceptable to approximate $\tilde{\mathcal{H}}_{\text{mis}}$ and to reduce the uncertainty of the system, which will be shown later. Meanwhile, the number of time steps does not influence the decay rate of $\tilde{\mathcal{H}}_{\text{mis}}$. This makes it possible to compute a fixed number of Arnoldi vectors for even more time steps.

We have seen from Fig. 1 that the eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$ have a sufficiently fast decay. Meanwhile, each eigenvector also exhibits a low-rank property. Here we plot the maximum rank used to compute the low-rank approximation of eigenvectors corresponding to the $50$ largest eigenvalues for $n_{t}=30,\ 60,\ 90$, respectively. In Fig. 2, it is shown that the increase of time steps keeps the rank bounded for the low-rank approximation of the eigenvectors of $\tilde{\mathcal{H}}_{\text{mis}}$. The threshold for the low-rank approximation is set to be $10^{-8}$.

As illustrated in Section 5, the eigenvector also admits a low-rank property. We perform a low-rank approximation on each Arnoldi vector throughout the Arnoldi iteration. Since the low-rank approximation is employed, the orthogonality of the basis of Arnoldi vectors is lost. We just need to compute a few Arnoldi vectors in practice. Here, we use a modified Gram-Schmidt method to perform the full reorthogonalization. Other types of reorthogonalization such as selective orthogonalization [13] or periodic orthogonalization [12] are also possible.

Here we use examples discretized by a $32\times 32$ grid, with $30$ and $60$ time steps to illustrate the effectiveness of the low-rank Arnoldi method. First, we plot the largest $40$ eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$ computed using the low-rank Arnoldi method. Next, we compute $\tilde{\mathcal{H}}_{\text{mis}}$ explicitly and use ‘eigs’ in MATLAB to compute its $40$ largest eigenvalues. These results are shown in Fig. 3. It is clearly shown that the low-rank Arnoldi method can recover the eigenvalues exactly by using reorthogonalization.

As shown in Fig. 1, the increase of time steps does not influence the decay rate of the eigenvectors of $\tilde{\mathcal{H}}_{\text{mis}}$. Next we will show that an increase of the spatial parameters does not influence the decay rate of eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$ either.

We fix the number of time steps $n_{t}$ to $30$, and compute $60$ largest eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$, the results are shown in Fig. 4a. The maximum rank used for the low-rank Arnoldi method with different number of spatial parameters is shown in Fig. 4b.

Figure 4a shows that the increase number of parameters do not influence the decay property of eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$. It is expected that for different number of parameters, the eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$ converge to the eigenvalues of the prior-preconditioned operator as long as the discretization of parameter field is good enough. This is illustrated by the eigenvalues shown in Fig. 4a. Meanwhile, maximum ranks used in the low-rank Arnoldi method are also bounded by a constant with the increase of number of parameters, which is shown in Fig. 4b.

As stated before, a threshold of $10^{-1}$ for the eigenvalue computations of $\tilde{\mathcal{H}}_{\text{mis}}$ is enough to reduce the uncertainty and to approximate the posterior covariance matrix. Next, we plot the diagonal entries of the approximated posterior covariance matrix, i.e., the variance of the points for a $64\times 64$ mesh with a different truncation threshold $\epsilon$ for eigenvalue computations of $\tilde{\mathcal{H}}_{\text{mis}}$. We use a $9$ sensors setting for the sparse observation inverse problem, where $9$ sensors are uniformly distributed inside the domain and the size of each sensor is $\nicefrac{{1}}{{256}}$ of the domain. Here we set $\beta_{\text{noise}}=10^{4}\beta_{\text{prior}}$ and the prior covariance matrix $\Gamma_{\text{prior}}=10I$, where $I$ is an identity matrix with appropriate size.

For the sparse observation case, covariance updates are mostly clustered around the area where data are observed while the rest are dominated by the prior. Uncertainty can only be reduced at areas around the location of sensors. This is clearly shown by Fig. 5a - Fig. 5d, where the dark colored areas are placed at the location of the sensors and have the lowest variance. Decreasing the threshold $\epsilon$, we observe that the variance is further reduced around the location of sensors. For smaller values of $\epsilon$ no more reduction in the variance is achieved as all essential information are already captured. Figure 5 shows that as long as the computations of $\Gamma_{\text{post}}$ is convergent, using a threshold $\epsilon=10^{-1}$ is enough to approximate the posterior covariance matrix and to reduce the uncertainty.

As analyzed in Section 5, the eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$ are related to $\frac{\beta_{\text{noise}}}{\beta_{\text{prior}}}$ and this ratio gives different updates of the posterior covariance matrix. Next, we use different $\frac{\beta_{\text{noise}}}{\beta_{\text{prior}}}$ ratios to plot the variance of the parameters. These results are shown in Figure 6 - Figure 7. The prior covariance matrix is set to be $\Gamma_{\text{prior}}=10I$, where $I$ is an identity matrix with appropriate size.

For the case $\beta_{\text{noise}}=10^{6}\beta_{\text{prior}}$, we need $72$ Arnoldi iterations for $\epsilon=10^{0}$ while $163$ Arnoldi iterations for $\epsilon=10^{-1}$. With $\beta_{\text{noise}}=10^{8}\beta_{\text{prior}}$, we need $347$ Arnoldi iterations for $\epsilon=10^{0}$ while $470$ Arnoldi iterations for $\epsilon=10^{-1}$.

Fig. 6 and Fig. 7 show that with the increase of the ratio between $\beta_{\text{noise}}$ and $\beta_{\text{prior}}$, the diagonal entries of $\Gamma_{\text{post}}$ becomes smaller. This implies that to further reduce the posterior variance, we need bigger ratio between $\beta_{\text{noise}}$ and $\beta_{\text{prior}}$. This can be explained as follows, the weight for the data misfit part in the optimization problem (26) is getting bigger for bigger $\frac{\beta_{\text{noise}}}{\beta_{\text{prior}}}$. This means that the data misfit part is more strictly optimized than for smaller $\frac{\beta_{\text{noise}}}{\beta_{\text{prior}}}$. Therefore the error between the estimation and observed data is getting smaller. However, this yields a more ill-conditioned problem and more Arnoldi iterations are needed. Therefore, a balance between covariance reduction and computational effort is needed with our approach enabling the storage of many Arnoldi vectors due to the complexity reduction of the low-rank approach.

6.3 Convection-Diffusion Equation

In this section, we study our low-rank approach for a stochastic inverse convection-diffusion problem. Here, the convection-diffusion operator $\mathcal{L}$ is given by

[TABLE]

The computational domain is chosen as a square domain given by $[0,\ 1]\times[0,\ 1]$, $\overrightarrow{\omega}=(0,\ 1)$, and the inflow is posed on the down boundary while the outflow is posed on the upper boundary. Boundary conditions are prescribed according to the analytic solution of the convection-diffusion equation, which is described as Example 3.3.1 in [11]. We use the streamline upwind Petrov–Galerkin (SUPG) finite element method to discretize the convection-diffusion equation.

First, we show the eigenvalue decay of $\tilde{\mathcal{H}}_{\text{mis}}$ for different settings of the viscosity parameter $\nu$. Here we set the number of time steps $n_{t}$ to be 30, and $\beta_{\text{noise}}=10^{4}\beta_{\text{prior}}$. We plot the $50$ largest eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$ for different $\nu$ in Fig. 8a.

As shown by Fig. 8a, the eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$ decay rapidly for big $\nu$ while this decay rate slows down when $\nu$ gets smaller. Therefore, more Arnoldi iterations are needed to get a satisfactory approximation of $\tilde{\mathcal{H}}_{\text{mis}}$. For smaller $\nu$, the largest eigenvalue is also bigger than that for bigger $\nu$ as shown in Fig. 8a. The first few eigenvectors form a more dominant subspace than for bigger $\nu$. It is therefore possible to choose a larger truncation threshold for smaller $\nu$, which will be shown later.

The maximum rank for the low-rank Arnoldi iteration with different $\nu$ is shown in Fig. 8b. It shows that the maximum rank does not increase with the decrease of $\nu$ and is bounded by a small constant. Therefore, the complexity for both computations and storage is $\mathcal{O}(n_{x}+n_{t})$ for the low-rank Arnoldi approach.

As shown in Fig. 8a, the eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$ decay slower for smaller $\nu$, and more Arnoldi iterations are needed, which is due to the property of the problem. For such a problem with smaller $\nu$, our low-rank approach is much more superior to the standard Arnoldi method introduced in [12] since we need to compute and store more Arnoldi vectors. Note that this is doable with the approach presented here.

Next, we set $\nu$ to be $10^{-2}$ and compute the $50$ largest eigenvalues of $\tilde{\mathcal{H}}_{\text{mis}}$ with different $n_{t}$, which are shown in Fig. 9a. We are also interested in the relation between the maximum rank at each Arnoldi iteration and the number of time steps $n_{t}$. This is shown in Fig. 9b.

As shown by Fig. 9a, with the increase in the number of time steps, the eigenvalue decay of $\tilde{\mathcal{H}}_{\text{mis}}$ behaves similar. Maximum ranks at each low-rank Arnoldi iteration with different $n_{t}$ are bounded by a moderately small constant, which is independent of $n_{t}$.

Fig. 8b and Fig. 9b show that the maximum rank for each low-rank Arnoldi iteration is almost invariant w.r.t. the number of time steps $n_{t}$ and the viscosity parameter $\nu$. This makes our low-rank Arnoldi method quite appealing for even complicated stochastic convection dominated inverse problems over a long time horizon.

Next, we show the diagonal entries of $\Gamma_{\text{post}}$ for different settings of $\nu$ and the threshold ($\epsilon$) of eigenvalues truncation. Here we set the number of time steps $n_{t}$ to be $90$, use a $32\times 32$ uniform grid to discretize the convection-diffusion equation, and $\beta_{\text{noise}}=10^{4}\beta_{\text{prior}}$. The results are given by Fig. 10 and Fig. 11.

For the case $\nu=10^{-2}$, we need $17$ Arnoldi iterations when we use a threshold $\epsilon=10^{1}$, while we need $121$ Arnoldi iterations for $\epsilon=10^{0}$. For the case $\nu=10^{-3}$, we need $52$ low-rank Arnoldi iterations by setting $\epsilon=10^{1}$ and $131$ low-rank Lanczos iterations when we use $\epsilon=10^{0}$. Note that often a further reduction in the truncation parameter does not yield better results as all the essential information is already captured.

Fig. 10b and Fig. 11b illustrate that the uncertainty (variance of unknowns) is already reduced dramatically even if we choose a relatively large threshold.

As shown in Fig. 8a, the smaller $\nu$ becomes, the bigger the largest eigenvalue of $\tilde{\mathcal{H}}_{\text{mis}}$ is going to be. Therefore, the first few eigenvectors for smaller $\nu$ form a more dominant subspace. This in turn implies that the uncertainty is much more reduced for smaller $\nu$ when we use the same truncation threshold of eigenvalues. We observe this in Fig. 10 and Fig. 11.

We can conclude that for the stochastic convection-diffusion inverse problem, our low-rank Arnoldi approach is very flexible and efficient for different time horizon lengths and viscosity parameters. It is even preferred for convection dominated stochastic inverse problems with long time horizon.

7 Conclusions

In this manuscript, we propose a low-rank Arnoldi method to approximate the posterior covariance matrix that appears in stochastic inverse problems. Compared with the standard Arnoldi approach, our approach exploits the low-rank property of each Arnoldi vector and makes a low-rank approximation of such a vector. This reduces the complexity for both computations and storage demand from $\mathcal{O}(n_{x}n_{t})$ to $\mathcal{O}(n_{x}+n_{t})$. Here $n_{x}$ is the degree of freedom in space and $n_{t}$ is the degree of freedom in time. This makes solving large scale stochastic inverse problems possible.

Our low-rank approach introduced in this manuscript solves linear stochastic inverse problems that can be put into the Bayesian framework. The next step of our work is to extend the low-rank approach introduced in this manuscript to nonlinear stochastic inverse problems, which is still a big challenge.

Acknowledgements

This work was supported by the European Regional Development Fund (ERDF/ EFRE: ZS/2016/04/78156) within the research center dynamic systems (CDS).