Efficient generalized Golub-Kahan based methods for dynamic inverse problems

TL;DR

This paper introduces efficient, matrix-free iterative methods based on the generalized Golub-Kahan bidiagonalization for solving large-scale dynamic inverse problems, enabling uncertainty estimation and regularization without explicit prior covariance matrix computations.

Contribution

The work develops novel, flexible, and scalable algorithms for dynamic inverse problems that incorporate prior information efficiently without explicit covariance matrix operations.

Findings

Methods successfully applied to photoacoustic tomography, deblurring, and seismic tomography.

Achieved rapid solutions for large-scale problems with millions of unknowns.

Demonstrated ability to process real and synthetic data efficiently.

Abstract

We consider efficient methods for computing solutions to and estimating uncertainties in dynamic inverse problems, where the parameters of interest may change during the measurement procedure. Compared to static inverse problems, incorporating prior information in both space and time in a Bayesian framework can become computationally intensive, in part, due to the large number of unknown parameters. In these problems, explicit computation of the square root and/or inverse of the prior covariance matrix is not possible. In this work, we develop efficient, iterative, matrix-free methods based on the generalized Golub-Kahan bidiagonalization that allow automatic regularization parameter and variance estimation. We demonstrate that these methods can be more flexible than standard methods and develop efficient implementations that can exploit structure in the prior, as well as possible structure in the forward model. Numerical examples from photoacoustic tomography, deblurring, and passive seismic tomography demonstrate the range of applicability and effectiveness of the described approaches. Specifically, in passive seismic tomography, we demonstrate our approach on both synthetic and real data. To demonstrate the scalability of our algorithm, we solve a dynamic inverse problem with approximately $43,000$ measurements and $7.8$ million unknowns in under $40$ seconds on a standard desktop.

Figures

Tables

Table 1. Table 1: Regularization parameters computed using WGCV and relative errors for PAT reconstructions corresponding to reorthogonalization. CPU time (in seconds) to obtain dynamic PAT reconstructions without and with reorthogonalization of gen-GK vectors.

	${\lambda }_{\mathrm{wgcv}}$	relative error	without reorth (sec)	with reorth (sec)
Matérn	40.32	2.9923e-01	84.5	273.7
${𝐐}_t=𝐈$	6.24	6.4575e-01	36.5	107.9
${𝐐}_t\ne 𝐈$	35.93	2.3406e-01	36.7	105.4

Table 2. Table 2: Regularization parameters and stopping iteration in decoupled λ wgcv ( i ) superscript subscript 𝜆 wgcv 𝑖 \lambda_{\rm wgcv}^{(i)} . WGCV regularization parameter (along with stopping iteration for the simultaneous approach) and the optimal regularization parameter are provided for comparison.

$i$	1	2	3	4	5	6	7	8	9
${\lambda }_{\mathrm{wgcv}}^{\left(i\right)}$	5.4689	5.1159	4.3018	4.0963	2.3008	4.1293	1.5994	0.6247	0.3141
$k_{\mathrm{stop}}^{\left(i\right)}$	61	82	55	24	17	5	5	4	5
${\lambda }_{\mathrm{wgcv}}$	8.5838 ($k_{\mathrm{stop}}=22$)
${\lambda }_{\mathrm{opt}}$	2.2812

Table 3. Table 3: Dates and number of rays for 8 time periods studied

Begin Date	End Date	Number of rays
Jan 14	Jan. 24	$m_1=13604$
Jan 25	Jan 31	$m_2=\mathrm{\hspace{0.17em}\hspace{0.17em}8512}$
Feb 1	Feb 17	$m_3=\mathrm{\hspace{0.17em}\hspace{0.17em}9255}$
Feb 8	Feb 14	$m_4=27121$
Feb 15	Feb 21	$m_5=41282$
Feb 22	Feb 28	$m_6=13755$
Mar 1	Mar 17	$m_7=11774$
Mar 8	Mar 16	$m_8=66553$

Equations

$$\mathbf{s}=\left[\begin{array}{c}{\mathbf{s}}_1\\ ⋮\\ {\mathbf{s}}_{n_t}\end{array}\right],\text{mbox}and\mathbf{d}=\left[\begin{array}{c}{\mathbf{d}}_1\\ ⋮\\ {\mathbf{d}}_{n_t}\end{array}\right],$$

$$\mathbf{d}=\mathbf{A}\mathbf{s}+\mathbf{\epsilon },$$

$$\pi (\mathbf{s}\mid \mathbf{d})\propto$$

$$\propto$$

$$\mathbf{s}(\lambda )=\underset{\mathbf{s}\in {\mathbb{R}}^{n_s{n}_t}}{\mathrm{argmin}}\frac{1}{2}\Vert \mathbf{A}\mathbf{s}-\mathbf{d}{\Vert }_{{\mathbf{R}}^{-1}}^2+\frac{{\lambda }^2}{2}\Vert \mathbf{s}-\mathbf{\mu }{\Vert }_{{\mathbf{Q}}^{-1}}^2,$$

$$({\mathbf{A}}^{\top }{\mathbf{R}}^{-1}\mathbf{A}+{\lambda }^2{\mathbf{Q}}^{-1})\mathbf{s}={\mathbf{A}}^{\top }{\mathbf{R}}^{-1}\mathbf{d}+{\lambda }^2{\mathbf{Q}}^{-1}\mathbf{\mu }.$$

$$\underset{\mathbf{s}}{\min }\frac{1}{2}\Vert {\mathbf{L}}_{\mathbf{R}}(\mathbf{A}\mathbf{s}-\mathbf{d}){\Vert }_2^2+\frac{{\lambda }^2}{2}\Vert {\mathbf{L}}_{\mathbf{Q}}(\mathbf{s}-\mathbf{\mu }){\Vert }_2^2,$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^{n_s{n}_t}}{\min }\frac{1}{2}{\Vert {\mathbf{L}}_{\mathbf{R}}(\mathbf{A}{\mathbf{L}}_{\mathbf{Q}}^{-1}\mathbf{x}-\mathbf{b})\Vert }_2^2+\frac{{\lambda }^2}{2}{\Vert \mathbf{x}\Vert }_2^2,$$

$$\mathbf{x}\leftarrow {\mathbf{Q}}^{-1}(\mathbf{s}-\mathbf{\mu }),\mathbf{b}\leftarrow \mathbf{d}-\mathbf{A}\mathbf{\mu },$$

$$({\mathbf{A}}^{\top }{\mathbf{R}}^{-1}\mathbf{A}\mathbf{Q}+{\lambda }^2\mathbf{I})\mathbf{x}={\mathbf{A}}^{\top }{\mathbf{R}}^{-1}\mathbf{b}.$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^{n_s{n}_t}}{\min }\frac{1}{2}\Vert \mathbf{A}\mathbf{Q}\mathbf{x}-\mathbf{b}{\Vert }_{{\mathbf{R}}^{-1}}^2+\frac{{\lambda }^2}{2}\Vert \mathbf{x}{\Vert }_{\mathbf{Q}}^2.$$

$$(\mathbf{A}\mathbf{Q}{\mathbf{A}}^{\top }+{\lambda }^2\mathbf{R})\mathbf{\xi }(\lambda )=\mathbf{d}-\mathbf{A}\mathbf{\mu }.$$

$$\mathbf{Q}{\mathbf{A}}^{\top }=$$

$$\mathbf{A}\mathbf{Q}{\mathbf{A}}^{\top }=$$

$${\beta }_{k+1}{\mathbf{u}}_{k+1}$$

$${\alpha }_{k+1}{\mathbf{v}}_{k+1}$$

$${\mathbf{B}}_k\equiv \left[\begin{array}{cccc}{\alpha }_1& & & \\ {\beta }_2& {\alpha }_2& & \\ & {\beta }_3& \ddots & \\ & & \ddots & {\alpha }_k\\ & & & {\beta }_{k+1}\end{array}\right],{\mathbf{U}}_{k+1}\equiv [{\mathbf{u}}_1,\dots ,{\mathbf{u}}_{k+1}],\text{mbox}and{\mathbf{V}}_k\equiv [{\mathbf{v}}_1,\dots ,{\mathbf{v}}_k],$$

$${\mathbf{U}}_{k+1}{\beta }_1{\mathbf{e}}_1=$$

$$\mathbf{A}\mathbf{Q}{\mathbf{V}}_k=$$

$${\mathbf{A}}^{\top }{\mathbf{R}}^{-1}{\mathbf{U}}_{k+1}=$$

$${\mathbf{U}}_{k+1}^{\top }{\mathbf{R}}^{-1}{\mathbf{U}}_{k+1}={\mathbf{I}}_{k+1}\text{mbox}and{\mathbf{V}}_k^{\top }\mathbf{Q}{\mathbf{V}}_k={\mathbf{I}}_k.$$

$${\mathbf{x}}_k\in \text{Span}\{{\mathbf{V}}_k\}={\mathcal{K}}_k({\mathbf{A}}^{\top }{\mathbf{R}}^{-1}\mathbf{A}\mathbf{Q},{\mathbf{A}}^{\top }{\mathbf{R}}^{-1}\mathbf{b})\equiv {\mathcal{S}}_k.$$

$${\mathbf{r}}_k\equiv \mathbf{A}\mathbf{Q}{\mathbf{x}}_k-\mathbf{b}={\mathbf{U}}_{k+1}\left({\mathbf{B}}_k{\mathbf{z}}_k-{\beta }_1{\mathbf{e}}_1\right).$$

$$\underset{{\mathbf{x}}_k\in {\mathcal{S}}_k}{\min }\frac{1}{2}\Vert {\mathbf{r}}_k{\Vert }_{{\mathbf{R}}^{-1}}^2+\frac{{\lambda }^2}{2}\Vert {\mathbf{x}}_k{\Vert }_{\mathbf{Q}}^2\iff \underset{{\mathbf{z}}_k\in {\mathbf{R}}^k}{\min }\frac{1}{2}{\Vert {\mathbf{B}}_k{\mathbf{z}}_k-{\beta }_1{\mathbf{e}}_1\Vert }_2^2+\frac{{\lambda }^2}{2}{\Vert {\mathbf{z}}_k\Vert }_2^2,$$

$${\mathbf{s}}_k=\mathbf{\mu }+\mathbf{Q}{\mathbf{x}}_k=\mathbf{\mu }+\mathbf{Q}{\mathbf{V}}_k{\mathbf{z}}_k,$$

$$G(\lambda )=\frac{n\Vert \mathbf{A}\mathbf{s}(\lambda )-\mathbf{d}{\Vert }_{{\mathbf{R}}^{-1}}^2}{{\left[\mathsf{trace}({\mathbf{I}}_m-{\mathbf{L}}_{\mathbf{R}}\mathbf{A}{\mathbf{A}}_{\lambda }^{†})\right]}^2},$$

$$G_{\text{proj}}(\lambda )\equiv \frac{k{\Vert (\mathbf{I}-{\mathbf{B}}_k{\mathbf{B}}_{k,\lambda }^{†}){\beta }_1{\mathbf{e}}_1\Vert }_2^2}{{\left[\mathsf{trace}({\mathbf{I}}_{k+1}-{\mathbf{B}}_k{\mathbf{B}}_{k,\lambda }^{†})\right]}^2},$$

$$\text{cov}\{Z({\mathbf{p}}_1,t_1),Z({\mathbf{p}}_2,t_2)\}=C({\mathbf{p}}_1-{\mathbf{p}}_2,t_1-t_2),$$

$$C(\mathbf{p},t)=C_S(\Vert \mathbf{p}\Vert )C_T(|t|)\forall (\mathbf{p},t)\in {\mathbb{R}}^d\times \mathbb{R},$$

$$C(\mathbf{p},t)=a_0{C}_S^0(\Vert \mathbf{p}\Vert )C_T^0(|t|)+a_1{C}_S^1(\Vert \mathbf{p}\Vert )+a_2{C}_T^2(|t|)\forall (\mathbf{p},t)\in {\mathbb{R}}^d\times \mathbb{R}.$$

$$C_{\nu ,\ell }(r)=\frac{1}{2^{\nu -1}\Gamma (\nu )}{\left(\frac{r\sqrt{2\nu }}{\ell }\right)}^{\nu }{K}_{\nu }\left(\frac{r\sqrt{2\nu }}{\ell }\right)$$

$$\kappa (r)=\exp \left(-(r/\ell {)}^{\gamma }\right)0<\gamma \le 2.$$

$${\mathbf{s}}_{i+1}={\mathbf{s}}_i+{\mathbf{ϵ}}_i{\mathbf{ϵ}}_i\sim \mathcal{N}(0,{\mathbf{Q}}_{\mathbf{s}}),$$

Full text

Keywords: dynamic inversion, Bayesian methods, Tikhonov regularization, generalized Golub-Kahan, Matérn covariance kernels, tomographic reconstruction

1 Introduction

The goal of an inverse problem is to use data, that is collected or measured, to estimate unknown parameters given some assumptions about the forward model [51, 22]. In many applications, the problem is assumed to be static, in the sense that the underlying parameters do not change during the measurement process. However, in many realistic scenarios such as in passive seismic tomography [53, 54] or dynamic electrical impedance tomography [44, 45], the underlying parameters of interest change dynamically. Incorporating prior information regarding temporal smoothness in reconstruction algorithms can lead to better reconstructions. However, this presents a significant computational challenge since many large spatial reconstructions may need to be computed, e.g., at each time step and for many time points. For example, in passive seismic tomography, geophones are used to collect measurements from seismic events (e.g., earthquakes) occurring in $1-2$ week intervals over $3-4$ months, and the goal is to obtain 3-dimensional spatial reconstructions of the elastic properties of the sub-surface for each time interval (e.g., to monitor changing stress conditions). As an other application, in medical imaging, during the data acquisition process the reconstruction algorithms need to account for patient motion and this can be modeled as a dynamic inverse problem.

Here we consider a discrete dynamic inverse problem with unknowns in space and time where the goal is to reconstruct parameters ${\bf s}_{i}\in\mathbb{R}^{n_{s}}$ from observations ${\bf d}_{i}\in\mathbb{R}^{m_{i}}$ for $i=1,\ldots,n_{t}$. Here $n_{s}$ refers to the number of spatial grid points, $n_{t}$ the number of time points, and $m=\sum_{i}m_{i}$ is the total number of measurements over all time points. For some problems, the number of measurements may correspond to the number of sensors or spatial measurement locations and thus may be the same for all time points. Let

[TABLE]

then we are interested in the following dynamic inverse problem,

[TABLE]

where ${\bf A}\in\mathbb{R}^{m\times n_{s}n_{t}}$ models the forward process which is assumed linear and ${\boldsymbol{\varepsilon}}$ represents noise or measurement errors in the data. We assume that ${\boldsymbol{\varepsilon}}\sim\mathcal{N}({\bf 0},{\bf R}),$ where ${\bf R}$ is a positive definite matrix whose inverse and square root are inexpensive (e.g., a diagonal matrix with positive diagonal entries). Given ${\bf A}$ and ${\bf d},$ the goal of the inverse problem is to reconstruct ${\bf s}.$ Since these problems are typically ill-posed, regularization is often required to compute a reasonable solution.

To solve the dynamic inverse problem we adopt the Bayesian approach. In this approach, the measured data and the parameters to be recovered (here, the spacetime unknowns) are modeled as random variables. Additionally we assume that the prior distribution for ${\bf s}$ is modeled as a Gaussian distribution with mean ${\boldsymbol{\mu}}=\begin{bmatrix}{\boldsymbol{\mu}}_{1}^{\top}&\cdots&{\boldsymbol{\mu}}_{n_{t}}^{\top}\end{bmatrix}^{\top}$ and positive-definite covariance matrix ${\bf Q}.$ That is, ${\bf s}\in\mathcal{N}({\boldsymbol{\mu}},\lambda^{-2}{\bf Q})$, where $\lambda$ is a (yet to be determined) scaling parameter for the precision matrix. For dynamic problems, covariance matrices ${\bf R}$ and ${\bf Q}$ contain information in both space and time. Then Bayes’ rule is used to combine the likelihood and the prior distribution and the posterior distribution,

[TABLE]

where $\|{\bf x}\|_{{\bf M}}=\sqrt{{\bf x}^{\top}{\bf M}{\bf x}}$ for any symmetric positive definite matrix ${\bf M}$. The maximum a posteriori (MAP) estimate, which is the peak of the posterior distribution, can be obtained by minimizing the negative log likelihood of (2), i.e.,

[TABLE]

Notice that for dynamic inverse problems, computing the MAP estimate requires solving for $n_{s}n_{t}$ unknowns. The main challenge here is that for the applications under consideration, $n_{s}$ is typically $O(10^{5}-10^{6})$ and $n_{t}$ is typically $O(10^{2}-10^{3})$. The resulting prior covariance matrices have $(n_{s}n_{t})^{2}$ or $O(10^{14}-10^{18})$ entries. Storing such covariance matrices is completely infeasible, much less performing computations with them. Clearly there is a need for developing specialized numerical methods for tackling the immense computational challenges arising from dynamic inverse problems. Our strategy is to use a combination of highly structured representations of prior information along with efficient numerical methods that can exploit these representations.

Overview of main contributions.

In this work, we adopt a Bayesian framework for solving dynamic inverse problems. We derive two efficient methods for computing MAP estimates, where the distinguishing features of our approach compared to previous methods are that we can incorporate a wide class of spatiotemporal priors, include time-dependent observation operators, and enable automatic regularization parameter selection. The resulting solvers are highly efficient and scalable to large problem sizes. In addition to the MAP estimate, we develop an efficient representation of the posterior covariance matrix using the generalized Golub-Kahan (gen-GK) bidiagonalization. This low-rank approximation can be used for uncertainty quantification, by estimating the variance of the distribution. Using several real-world imaging applications (including both simulated and real data), we show that our methods are well suited for a wide class of dynamic inverse problems and we demonstrate scalability of our algorithms.

Related work.

The literature on dynamic inverse problems is large, and it is not our intention to provide a detailed overview. We mention a few related approaches that are relevant to our work.

A popular approach for spacetime reconstructions is the use of Kalman filters and smoothers. However, textbook implementations of these methods can be prohibitively expensive. This is because they require the storage and computation of covariance matrices that scale as $\mathcal{O}(n_{s}^{2})$. One approach is to use an efficient representation of the state covariance matrix, as a low-rank perturbation of an appropriately chosen matrix [35, 36]. Efficient computational techniques for the Kalman filter especially tailored to the random-walk forecast model were proposed by [29, 43]. In Section 2.4, we show how this forecast model fits within our framework.

Significant simplifications can be made if we assume that measurement errors are independent in time and reconstruct the parameters of interest only using the data available from the current time step. However, several other authors, see for example Schmitt and collaborators [44, 45], have emphasized the importance of including temporal priors in many practical applications. They considered a total-variation type temporal smoothness prior with a simple spatial prior (the identity matrix) and showed that their approach achieved superior results in faster computational time than other statistical approaches such as Kalman smoothers.

Our approach is more general in that we allow for a variety of spatial priors where the resulting covariance matrices are dense, unwieldy, and only available via matrix-vector multiplication. Furthermore, we consider more general forward models and consider hybrid iterative approaches so that the regularization parameter $\lambda$ can be automatically estimated. In previous studies, the regularization (or precision) parameter decoupled in space and time, and the resulting two parameters were required algorithmic inputs [44].

In Section 2, we describe the problem set-up and address various changes of variables that can be used. We also provide a brief overview of generalized hybrid methods to efficiently solve static inverse problems. Efficient methods for approximating the MAP estimate, i.e., solving (3), in the space-time formulation will be described in Section 3, where special cases of problem structure will be considered for efficiency. Efficient variance estimation methods based on the gen-GK bidiagonalization are described in Section 4. Numerical results are presented in Section 5 for various simulated imaging problems and for real data from passive seismic tomography. Conclusions and discussions are provided in Section 6.

2 Problem set up and background

One main goal in the Bayesian framework is to efficiently compute the MAP estimate, and in this paper, we are mainly interested in cases where ${\bf Q}$ is a very large, dense matrix, so computing ${\bf Q}$ and ${\bf Q}^{-1}$ or their factorizations is not feasible. Such scenarios arise, for example, when working with Gaussian random fields, in which case forming ${\bf Q}$ explicitly may not be feasible, but computing matrix-vector products (MVPs) with ${\bf Q}$ can be done efficiently [42].

First, we describe various problem formulations for computing the MAP estimate, and describe a change of variables so that gen-GK methods can be used. Then, a brief overview of generalized hybrid methods is provided in Section 2.2 for completeness, and a discussion on various choices for modeling temporal priors is provided in Section 2.3. Connections to previous works are addressed in Section 2.4.

2.1 Problem formulations

Notice that the desired MAP estimate ${\bf s}(\lambda)$ is the solution to the system of equations

[TABLE]

For problems where matrix factorizations of ${\bf R}$ and ${\bf Q}$ are possible, a common approach to compute the MAP estimate is to solve the following general-form Tikhonov problem,

[TABLE]

where ${\bf Q}^{-1}={\bf L}_{\bf Q}^{\top}{\bf L}_{\bf Q}$ and ${\bf R}^{-1}={\bf L}_{\bf R}^{\top}{\bf L}_{\bf R}$. By transforming the problem to standard form, we get the priorconditioned problem [3, 4],

[TABLE]

where ${\bf s}={\boldsymbol{\mu}}+{\bf L}_{\bf Q}^{-1}{\bf x}$.

In the applications that we consider, the covariance matrices can be very large and dense, so the storage and computational costs to obtain factorizations and/or inverses of ${\bf Q}$ can be prohibitive. In order to avoid matrix factorizations of ${\bf Q}$ and/or expensive linear solves with ${\bf Q}$, a different change of variables was proposed in [8], where

[TABLE]

so that (4) reduces to the modified system of equations

[TABLE]

In summary, with this change of variables, the MAP estimate is given by ${\bf s}(\lambda)={\boldsymbol{\mu}}+{\bf Q}{\bf x}(\lambda)$, where ${\bf x}(\lambda)$ is the solution to the following optimization problem

[TABLE]

Hybrid iterative methods for approximating ${\bf x}(\lambda)$ were described in [8] for generic ${\bf A},{\bf R}$, and ${\bf Q}$. For completeness, we give a brief description of these methods in Section 2.2 and refer the interested reader to the paper for more details.

In this paper, we choose to focus on iterative methods; however, we briefly mention an alternative formulation. For relatively small problems where the number of overall measurements $m$ is small, i.e., $O(10^{3}-10^{4})$, a simple change of variables along with the Sherman-Morrison formula can be used to avoid ${\bf Q}^{-1}$. The MAP estimate can be computed as ${\bf s}(\lambda)={\boldsymbol{\mu}}+{\bf Q}{\bf A}^{\top}{\boldsymbol{\xi}}(\lambda)$ where

[TABLE]

Notice that in terms of solving linear systems, the number of unknowns has reduced from $n_{s}n_{t}$ to $m$. A direct solver could be used to solve (10), but forming ${\bf Q}{\bf A}^{\top}$ may be computationally prohibitive, costing $\mathcal{O}(n_{s}n_{t}m^{2}$). Iterative methods could be used, but in these cases, it may be difficult to know a good regularization parameter a priori. Further simplifications would be possible for problems where $m_{i}$ is constant for all $i=1,\dots,n_{t}$ and ${\bf A}$ is also a Kronecker product, i.e., ${\bf A}={\bf A}_{t}\otimes{\bf A}_{s}$ so that

[TABLE]

As mentioned earlier, we do not pursue this approach since it is computationally expensive when the number of measurements is large.

2.2 Generalized hybrid methods

Given matrices ${\bf A}$, ${\bf R}$, ${\bf Q}$, and vector ${\bf b},$ with initializations $\beta_{1}=\|{\bf b}\|_{{\bf R}^{-1}},{\bf u}_{1}={\bf b}/\beta_{1}$ and $\alpha_{1}{\bf v}_{1}={\bf A}^{\top}{\bf R}^{-1}{\bf u}_{1}$, the $k$th iteration of the gen-GK bidiagonalization procedure generates vectors ${\bf u}_{k+1}$ and ${\bf v}_{k+1}$ such that

[TABLE]

where scalars $\alpha_{i},\beta_{i}\geq 0$ are chosen such that $\|{\bf u}_{i}\|_{{\bf R}^{-1}}=\|{\bf v}_{i}\|_{{\bf Q}}=1$. At the end of $k$ steps, we have

[TABLE]

where the following relations hold up to machine precision,

[TABLE]

Furthermore, in exact arithmetic, matrices ${\bf U}_{k+1}$ and ${\bf V}_{k}$ satisfy the following orthogonality conditions

[TABLE]

An algorithm for the gen-GK bidiagonalization process is provided in Algorithm 1. In addition to MVPs with ${\bf A}$ and ${\bf A}^{\top}$ that are required for the standard GK bidiagonalization [17], each iteration of gen-GK bidiagonalization requires two MVPs with ${\bf Q}$ and two solves with ${\bf R}$ (which are assumed to be cheap); in particular, we emphasize that Algorithm 1 avoids ${\bf Q}^{-1}$ and ${\bf L}_{\bf Q}$, due to the change of variables in (7).

We seek solutions of the form ${\bf x}_{k}={\bf V}_{k}{\bf z}_{k}$, so that

[TABLE]

Define the residual at step $k$ as ${\bf r}_{k}\equiv\>{\bf A}{\bf Q}{\bf x}_{k}-{\bf b}$. It follows from Equations (11)-(13) that

[TABLE]

To obtain coefficients ${\bf z}_{k}$, we take ${\bf z}_{k}$ that minimizes the genLSQR problem,

[TABLE]

where the gen-GK relations were used to obtain the equivalence. Variants of this formulation, e.g., for LSMR, could be used as well [8]. After computing a solution to the projected problem, an approximate MAP estimate can be recovered by undoing the change of variables,

[TABLE]

where, now, ${\bf s}_{k}\in{\boldsymbol{\mu}}+{\bf Q}\mathcal{S}_{k}$. For fixed $\lambda$ and in exact arithmetic, iterates of the genLSQR approach are mathematically equivalent to some pre-existing solvers (e.g., filtered GSVD solutions and priorconditioned solutions) [8]. However, if $\lambda$ is not known a priori, hybrid methods can take advantage of the shift-invariance property of Krylov subspaces and select $\lambda$ adaptively and automatically by utilizing well-known SVD-based regularization parameter selection schemes [22, 20, 51] for the projected problem (15), since ${\bf B}_{k}$ is only of size $(k+1)\times k$. Henceforth, we refer to this approach as genHyBR. Previous work on parameter selection within hybrid methods include [5, 7, 25, 13, 38, 39].

Although a wide range of regularization parameter selection methods can be used in our framework, in this paper we consider a variant of the generalized cross validation (GCV) approach. The GCV parameter is selected to minimize the GCV function [16] corresponding to the general-form Tikhonov problem (5),

[TABLE]

where ${\bf A}_{\lambda}^{\dagger}=({\bf A}^{\top}{\bf R}^{-1}{\bf A}+\lambda^{2}{\bf Q}^{-1})^{-1}{\bf A}^{\top}{\bf L}_{\bf R}^{\top}$. We have assumed ${\boldsymbol{\mu}}={\bf 0}$ for simplicity. At the $k$th iteration, the GCV parameter corresponding to the projected problem (15) minimizes,

[TABLE]

where ${\bf B}_{k,\lambda}^{\dagger}=({\bf B}_{k}^{\top}{\bf B}_{k}+\lambda^{2}{\bf I})^{-1}{\bf B}_{k}^{\top}$. A weighted-GCV (WGCV) approach [5] has been suggested for use within hybrid methods, where a weighting parameter is introduced in the denominator of (18). We denote $\lambda_{\rm wgcv}$ to be the regularization parameter computed using WGCV. As a benchmark for simulated experiments, we also consider the optimal regularization parameter $\lambda_{\rm opt}$, which minimizes the $2$-norm of the error between the reconstruction and the truth.

2.3 Modeling prior covariances

Following the geostatistical approach, we model the unknown field $s({{\bf p}},t)$ as a realization of a spatio-temporal random function $Z({\bf p},t)$ for ${\bf p}\in\mathbb{R}^{d}$ and $t\in\mathbb{R}$. We assume that the covariance function is stationary in space and stationary in time; for simplicity here, we also assume that the mean is zero. In other words, we assume that the covariance function satisfies

[TABLE]

where $C:\mathbb{R}^{d}\times\mathbb{R}\rightarrow\mathbb{R}_{+}$ is a positive definite covariance function.

We briefly review various formulations for spacetime covariance kernels, before delving into the specific choices of kernels we make. Perhaps the most convenient representation can be obtained if we make the assumption that the covariance function is separable in space and time, and isotropic in these variables, then $C$ takes the form

[TABLE]

where $C_{S}(\|{\bf p}\|)$ and $C_{T}(|t|)$ are isotropic, purely spatial and purely temporal covariance functions, respectively. It can be readily seen that the resulting matrices ${\bf Q}$ have the Kronecker product structure. The separability assumption is common in the statistics literature [27, 14], and tests for separability can be found in [12]. The Kronecker product structure has computational advantages which we will exploit in Section 3 to develop efficient algorithms. While mathematically and computationally convenient, it is important to recognize the potential shortcomings of the separability assumption. The key issue is that the prior models do not allow for interactions in variability of space and time; for a detailed discussion, see [27]. However, even if the covariance kernel is not separable, one may approximate it using a separable covariance kernel [14], where the resulting covariance matrix approximation can be represented as a Kronecker product or a sum of Kronecker products. Such approximations were studied in [49] and have been shown to be successful in the context of image deblurring, e.g., [23, 33, 9].

Many nonseparable covariance kernels have been proposed that model space-time interactions of variability. One such approach uses $C({\bf p},t)=\varphi(\sqrt{c_{1}\|{\bf p}\|^{2}+c_{2}|t|^{2}})$, where $c_{1},c_{2}$ are weights that control the correlation of the space and time variables, and $\varphi(\cdot)$ is an appropriate covariance kernel. Another approach is to use a product-sum model

[TABLE]

where $a_{0}$, $a_{1}$ and $a_{2}$ are nonnegative coefficients and $C_{S}^{0}$, $C_{S}^{1}$ and $C_{T}^{0}$, $C_{T}^{2}$ are isotropic, purely spatial and purely temporal covariance functions, respectively. A review of these covariance kernels is provided in [15].

A wide variety of choices for ${\bf Q}$ can be included in our framework and thus incorporated in the methods described below. Next we give a few examples of temporal and spatial priors that are well-suited for our problems.

Specific choices of covariance kernels.

A common approach is to use Gaussian random fields where the entries of the covariance matrix are computed directly as $({\bf Q}_{t})_{ij}=\kappa(|t_{i}-t_{j}|)$, where $\{t_{i}\}_{i=1}^{n_{t}}$ are the time points. A popular choice for $\kappa(\cdot)$ is from the Matérn family of covariance kernels [37], which form an isotropic, stationary, positive-definite class of covariance kernels. We define the covariance kernel in the Matérn class as

[TABLE]

where $\Gamma$ is the Gamma function, $K_{\nu}(\cdot)$ is the modified Bessel function of the second kind of order $\nu$, and $\ell$ is a scaling factor. The choice of parameter $\nu$ in equation (19) defines a special form for the covariance. For example, when $\nu=1/2$, $C_{\nu,\ell}$ corresponds to the exponential covariance function, and if $\nu=1/2+p$ where $p$ is a non-negative integer, $C_{\nu,\ell}$ is the product of an exponential covariance and a polynomial of order $p$. Also, in the limit as $\nu\rightarrow\infty$, $C_{\nu,\ell}$ converges to the Gaussian covariance kernel, for an appropriate scaling of $\ell$. Another related family of covariance kernels is the $\gamma$-exponential function [37],

[TABLE]

Just as $C_{T}(\cdot)$ can be a $\gamma$-exponential function, or chosen from the Matérn class, $C_{S}(\cdot)$ can also be chosen in the same way. Therefore, ${\bf Q}_{s}$ has entries $({\bf Q}_{s})_{ij}=\kappa(\|{\bf p}_{i}-{\bf p}_{j}\|)$, where $\{{\bf p}_{i}\}_{i=1}^{n_{s}}$ are the spatial locations. In the applications of interest, the number of spatial locations $n_{s}$ is much larger than the number of time points $n_{t}$; thus the storage of ${\bf Q}_{s}$ is challenging, as is employing it in iterative methods, since the cost of an MVP is $\mathcal{O}(n_{s}^{2})$. Both the storage and computational cost can be reduced to $\mathcal{O}(n_{s}\log n_{s})$ using the FFT based approach or $\mathcal{H}$-matrix approach. This has been reviewed in [41, 40].

2.4 Other examples that fit our framework

As mentioned in the introduction, the random-walk forecast model [50, 26, 47, 34] was previously considered for its computational advantages [29, 43]. We show how this model also fits within our framework. Assume that the state ${\bf s}_{i}$ undergoes the following dynamics for $i=1,\dots,n_{t}-1$

[TABLE]

with initial conditions ${\bf s}_{1}\sim\mathcal{N}(\bf 0,{\bf Q}_{s})$. We can then express the distribution of the state ${\bf s}$ as

[TABLE]

Thus, ${\bf s}$ is a Gaussian distribution with zero mean and covariance matrix ${\bf Q}_{t}\otimes{\bf Q}_{s}$, where

[TABLE]

The matrix ${\bf Q}_{t}$ has an explicit representation and is the so-called minij matrix with $(i,j)$-th entry of ${\bf Q}_{t}$ equal to $\min\{i,j\}$. A similar representation is also available for the forecast model ${{\bf s}_{i+1}=\bar{\alpha_{i}}{\bf s}_{i}+\bar{\beta_{i}}{\boldsymbol{\epsilon}}_{i}}$, but will not be considered here.

Another approach assumes that ${\bf Q}_{t}=({\bf L}_{t}^{\top}{\bf L}_{t}+\gamma{\bf I})^{-1}$ where ${\bf L}_{t}$ is a sparse discretization of a differential operator and $\gamma$ is a small positive parameter to ensure ${\bf Q}_{t}$ is positive definite. For example, a common choice is to enforce smoothness in time by selecting

[TABLE]

Although not derived within a Bayesian framework, Schmitt and collaborators [44, 45] considered such temporal priors along with a standard Tikhonov term to enforce spatial smoothness. In fact, it is possible to show that their algorithm approximates the MAP estimate, which in our framework corresponds to ${\bf R}={\bf I},{\boldsymbol{\mu}}={\bf 0},$ and

[TABLE]

where ${\bf B}={\bf L}_{t}\otimes{\bf I}.$ Here $\lambda_{s}$ and $\lambda_{t}$ correspond to regularization parameters in space and time respectively, and the MAP estimate minimizes the function,

[TABLE]

In this paper, we focus on Matérn kernels and their covariance matrices for both the spatial and temporal priors.

3 Generalized hybrid methods for dynamic inverse problems

In this section, we describe various approaches based on the gen-GK bidiagonalization for computing MAP estimates for dynamic inverse problems. These iterative approaches are desirable for problems where both ${\bf A}$ and ${\bf Q}$ may be extremely large, or for problems where these matrices are not explicitly stored but can be accessed via function calls to compute MVPs with ${\bf A},{\bf A}^{\top},$ and ${\bf Q}$ efficiently. In Section 3.1, we describe an “all-at-once” generalized hybrid method that requires MVPs with covariance matrix ${\bf Q}$. This requirement is quite general and includes matrices such as ${\bf Q}$ being a Kronecker product, a sum of Kronecker products, or a convolution operator. Then for problems where the number of time points is small and ${\bf A}$ and ${\bf R}$ are also Kronecker products, we describe an efficient decoupled approach in Section 3.2.

3.1 Simultaneous generalized hybrid approach

The first approach we consider for solving dynamic inverse problems is to use genHyBR as summarized in Section 2.2 to solve for all unknown variables (e.g., in space and time) simultaneously. Since the number of unknowns can be quite large in the “all-at-once” approach and the gen-GK vectors ${\bf V}_{k}$ must be stored for hybrid methods, we assume that solutions can be captured in relatively few iterations or that appropriate preconditioning can be used so that $k$ remains small. Next we describe efficiencies can can be gained for problems where ${\bf Q}$ is a Kronecker product, but we reiterate that the simultaneous approach has applicability beyond the cases presented here.

For problems where ${\bf Q}$ is a Kronecker product [28], MVPs with ${\bf Q}$ can be computed efficiently as

[TABLE]

where $\mathsf{vec}$ and $\mathsf{mat}$ are operations such that $\mathsf{vec}$ unfolds a matrix ${\bf X}\in\mathbb{R}^{m\times n}$ into a vector by stacking column-wise and $\mathsf{mat}$ folds it back, i.e.,

[TABLE]

Assuming the cost of an MVP with ${\bf Q}_{s}$ is $\mathcal{O}(n_{s}\log n_{s})$, the cost of ${\bf Q}{\bf x}$ is $\mathcal{O}(n_{t}n_{s}\log n_{s}+n_{s}n_{t}^{2})$, which is significantly smaller than the naive cost of $\mathcal{O}(n_{s}^{2}n_{t}^{2})$. Further reductions in computational cost can be achieved by parallelizing the matrix-matrix multiplications, e.g., MVPs of ${\bf Q}_{s}$ with the columns of ${\bf X}$.

Separable forward operator.

If, additionally, the number of measurements at each timestep is the same, which we denote by $\bar{m}$, and ${\bf A}={\bf A}_{t}\otimes{\bf A}_{s}$ and ${\bf R}={\bf R}_{t}\otimes{\bf R}_{s}$, where ${\bf R}_{t}\in\mathbb{R}^{n_{t}\times n_{t}}$ and ${\bf R}_{s}\in\mathbb{R}^{\bar{m}\times\bar{m}}$ are the temporal and spatial noise covariance matrices respectively, then MVPs required for the gen-GK bidiagonalization algorithm can be computed efficiently. That is, for vectors ${\bf x}\in\mathbb{R}^{n_{s}n_{t}}$ and $\boldsymbol{y}\textbullet\in\mathbb{R}^{\bar{m}n_{t}}$

[TABLE]

where ${\bf X}\equiv\mathsf{mat}({\bf X})\in\mathbb{R}^{n_{s}\times n_{t}}$ and ${\bf Y}\equiv\mathsf{mat}({\bf y})\in\mathbb{R}^{\bar{m}\times n_{t}}$.

3.2 Decoupled generalized hybrid approach

For problems where ${\bf A},{\bf R}$, and ${\bf Q}$ are all Kronecker products and $n_{t}$ is relatively small such that ${\bf Q}_{t}$ and its factor ${\bf L}_{t}$ are feasible, we describe a decoupled genHyBR approach. The normal equations corresponding to the weighted least-squares problem in (9) can be written as

[TABLE]

We can exploit the fact that ${\bf Q}={\bf L}_{t}^{\top}{\bf L}_{t}\otimes{\bf Q}_{s}\,,$ and let ${\bf y}=({\bf L}_{t}\otimes{\bf I}){\bf x}$, alternatively ${\bf Y}={\bf X}{\bf L}_{t}^{\top}$. We introduce the variable ${\bf H}={\bf A}^{\top}{\bf R}^{-1}{\bf A}$; similarly, we also define ${\bf H}_{s}={\bf A}^{\top}_{s}{\bf R}^{-1}_{s}{\bf A}_{s}$ and ${\bf H}_{t}={\bf A}_{t}^{\top}{\bf R}_{t}^{-1}{\bf A}_{t}$. The normal equations simplify, when we left-multiply by ${\bf L}_{t}^{-\top}\otimes{\bf I}$ to obtain

[TABLE]

Using the properties of Kronecker products, this equation can alternatively be written as the following generalized-Sylvester equation

[TABLE]

where, for simplicity, we introduce ${\bf C}\equiv{\bf Q}_{s}{\bf A}_{s}^{\top}{\bf R}_{s}^{-1}{\bf B}{\bf R}_{t}^{-1/2}$ and $\widehat{{\bf A}}_{t}\equiv{\bf R}_{t}^{-1/2}{\bf A}_{t}{\bf L}_{t}^{\top}$.

Let $\widehat{{\bf A}}_{t}={\bf U}_{t}{\boldsymbol{\Sigma}}_{t}{\bf V}_{t}^{\top}$ be its singular value decomposition, then ${\bf L}_{t}{\bf H}_{t}{\bf L}_{t}^{\top}={\bf V}_{t}{\boldsymbol{\Sigma}}_{t}^{2}{\bf V}^{\top}_{t}$. We make another change of variables ${\bf Z}\equiv{\bf Y}{\bf V}_{t}={\bf X}{\bf L}_{t}^{\top}{\bf V}_{t},$ and we get

[TABLE]

where we have multiplied on the right by ${\bf V}_{t}$. Expanding the above expression column-wise, the key observation is that all the equations decouple so that

[TABLE]

Notice that for $\sigma_{i}=0$, ${\bf z}_{i}=\bf 0$, whereas for $\sigma_{i}>0,$ the solution can be obtained by solving the least-squares problem,

[TABLE]

which can be done using the genHyBR method. Note that a transformation back of variables must be made. This is summarized in Algorithm 2.

This approach is embarrassingly parallel, which may lead to enormous computational savings. Additionally, solving the sequence of decoupled systems may be done efficiently either by recycling of Krylov subspaces or by effective preconditioning.

4 Estimation of posterior variances using gen-GK

For the problem considered here, the posterior distribution ${\bf s}|{\bf d}$ is Gaussian with

[TABLE]

where ${\bf H}={\bf A}^{\top}{\bf R}^{-1}{\bf A}$. The MAP estimate corresponds to the mode of the posterior distribution and provides information regarding the “most likely” estimate. On the other hand, the posterior variance defined as the diagonals of the posterior covariance $\boldsymbol{\Gamma}_{\text{post}}$ provides a measure of the spread of the posterior distribution around the posterior mean. For static inverse problems, estimating the posterior covariance matrix is known to be computationally challenging [41, 42]. For dynamic problems, the problem is further exacerbated since the posterior covariance matrix is of size $n_{s}n_{t}\times n_{s}n_{t}$. Moreover, this matrix is dense and forming it explicitly to obtain the diagonal entries is computationally infeasible. In this work, we use intermediate information from genHyBR for computing a MAP estimate to estimate the posterior variance. Following Section 3, we consider a simultaneous and a decoupled approach.

Before we explain how to compute variances, we make the following remark. In the sequel, we will use Algorithm 1 with one minor modification, namely, at each step we explicitly re-orthogonalize the vectors in ${\bf U}_{k}$ and ${\bf V}_{k}$. Numerical experience suggests that this marginally increases the computational cost by $\mathcal{O}(k^{2}(m+n))$ but considerably improves the accuracy of the variance computations.

4.1 Simultaneous approach

Recall that after $k$ iterations of the gen-GK bidiagonalization process, we have matrices ${\bf B}_{k},{\bf U}_{k}$ and ${\bf V}_{k}$ satisfying relations (11)-(13) and (14). Let ${\bf B}_{k}^{\top}{\bf B}_{k}={\bf W}_{k}{\boldsymbol{\Theta}}_{k}{\bf W}_{k}^{\top}$ be the eigenvalue decomposition with eigenvalues $\theta_{1},\ldots,\theta_{k}$ and let ${\bf Z}_{k}={\bf Q}{\bf V}_{k}{\bf W}_{k}$, then we get the following low-rank approximation

[TABLE]

Using (27) and the Woodbury formula, we obtain the approximation

[TABLE]

where

[TABLE]

Notice that we have an efficient representation of $\widehat{\boldsymbol{\Gamma}}_{\text{post}}$ as a low-rank perturbation of the prior $\lambda^{-2}{\bf Q}$. In summary, diagonal entries of $\widehat{\boldsymbol{\Gamma}}_{\text{post}}$ can provide estimates of diagonal entries of $\boldsymbol{\Gamma}_{\text{post}}$, where the main computational requirement is to obtain the diagonals of ${\bf Q}$ and the diagonals of the rank-$k$ perturbation. Therefore, the only additional computational cost for estimating the posterior variance is $\mathcal{O}(k^{3}+k^{2}n_{s}n_{t})$.

An approximation of this kind was previously explored in [42, 11, 1, 2]; however, the error estimates developed in the above references assume that the exact eigenpairs are available. If the Ritz pairs converge to the exact eigenpairs of the matrix ${\bf Q}{\bf H}{\bf Q}$, then furthermore, the optimality result in [48, Theorem 2.3] applies here as well.

4.2 Decoupled approach

For cases where ${\bf A}={\bf A}_{t}\otimes{\bf A}_{s},$ we develop a similar strategy for estimating the posterior variances by exploiting the decoupled structure described in Section 3.2. In contrast to the simultaneous approach, a different Krylov subspace is constructed for each time step. Using notation defined in Section 3.2, the posterior covariance matrix is given by

[TABLE]

where the matrix in the center is a block-diagonal matrix whose diagonal blocks are $(\lambda^{2}{\bf Q}_{s}^{-1}+\sigma_{i}^{2}{\bf H}_{s})^{-1}$ for $i=1,\ldots,n_{t}\,.$ Analogous to the simultaneous approach, gen-GK approximations for each $i$ denoted by ${\bf B}_{k,i}$ and ${\bf V}_{k,i}$ can be used to get low-rank approximations

[TABLE]

where ${\bf B}_{k,i}^{\top}{\bf B}_{k,i}={\bf W}_{k,i}{\boldsymbol{\Theta}}_{k,i}{\bf W}_{k,i}^{\top}$ is an eigenvalue decomposition,

[TABLE]

Notice that depending on the stopping criteria for the gen-GK process, each time point may have a different rank $k$. We omit this dependence for clarity of presentation, and denote the low-rank blocks ${\bf D}_{i}={\bf Z}_{k,i}{\boldsymbol{\Delta}}_{k,i}{\bf Z}_{k,i}^{\top}$.

In summary, an approximation to the posterior covariance matrix is given by

[TABLE]

where the diagonals can be computed for timestep $i=1,\dots,n_{t}$ as

[TABLE]

where $[{\bf Q}_{t}]_{ii}$ is the $i$th diagonal entry of ${\bf Q}_{t}$ and the operation $\mathsf{diag}\,(\cdot)$ returns a vector containing the diagonals of a matrix.

5 Numerical Results

In this section, we provide three examples from image reconstruction. The first is a model problem from dynamic photoacoustic tomography (PAT) reconstruction under motion that illustrates the significant impact of including a temporal prior. Then we consider an image deblurring problem that can exploit the decoupled framework of Section 3.2, and we finish with a severely ill-posed passive seismic tomography (PST) problem, where we include results from real field measurements. In all of the results, genHyBR solutions correspond to (16) where ${\bf x}_{k}$ is the solution to (15).

5.1 Dynamic photoacoustic tomography (PAT) under motion

PAT is a hybrid imaging modality that combines the rich contrast of optical imaging with the high resolution of ultrasound imaging, thereby producing higher resolution in-vivo images with lower patient risk (e.g., requiring no ionizing radiation and no contrast agents) and lower cost and inconvenience than other imaging modalities. In modern PAT systems, transducers are rotated around an object and data is acquired in time. Most current methods for PAT image reconstruction such as explicit inversion formulas, time reversal, or series solution are not suited for such systems. A significant limitation in extending these methods to dynamic PAT is that accurate tomographic reconstruction relies on expensive motion estimation and artifact removal software, since the object being imaged may move during data acquisition. Chung and Nguyen [6] recently studied a continuous model of PAT reconstruction under motion for specific parameterized motion models and developed specialized algorithms for these scenarios. However, in this work we consider PAT reconstruction for more general motion deformations and seek a sequence of reconstructions, rather than just one reconstruction, by incorporating a temporal prior and more informative image priors. We follow the literature in dynamic tomography (e.g., [24, 18, 19]), and note that since we do not consider parameterized models, our framework can incorporate other realistic scenarios (other than motion) such as non-stationary optical illumination and inaccurate transducer responses [46, 55, 52, 30].

We consider the discrete problem where ${\bf s}_{i}\in\mathbb{R}^{n_{s}}$ is the discretized desired solution111Here we assume the 2D image is vectorized column-wise. at time point $i$ and let $\left\{{\bf z}_{i}\right\},i=1,\ldots,n_{t}$ denote the locations of the transducers. At each transducer location ${\bf z}_{i}$, assume there are $r$ radii and let ${\bf A}_{i}\in\mathbb{R}^{r\times N}$ be the corresponding projection matrix for that location (i.e., ${\bf A}_{i}{\bf s}_{i}$ is the discrete circular Radon transform of ${\bf s}$ on circles centered at $z_{i}$). Thus, the observed spherical projection measurements for all $r$ radii are contained in vector

[TABLE]

where ${\bf e}_{i}\in\mathbb{R}^{r}$ is additive Gaussian noise that is independent and identically distributed. Then the forward model has the form

[TABLE]

where the goal is to estimate the desired images ${\bf s}_{i}$, given the observations ${{\bf d}_{i}.}$

For this example, we consider $120$ true images of size $256\times 256$ that were generated with two Gaussians with fixed width, and rotating counterclockwise. See Figure 1(a) for sample true images. Measurements were taken at $120$ equidistant angles between [math] and $357$ at $3$ degree intervals, and each projection consists of $363$ radii. White noise was added to the observations with ${\bf R}=.0082^{2}{\bf I}$ (this corresponds to a noise level of $0.04$). Results presented here use ${\bf R}$ as above, although noise estimation algorithms could be used [10]. The sinogram of size $363\times 120$ is shown in Figure 1(b), providing a total number of $43,560$ observations.

For the dynamic inverse problem, the number of unknowns is $256*256*120=7,864,320$. We provide, for comparison, a static reconstruction in Figure 1(c), where we solve the following (inaccurate) model problem,

[TABLE]

with ${\bf d}=\begin{bmatrix}{\bf d}_{1}^{\top}&\cdots&{\bf d}_{n_{t}}^{\top}\end{bmatrix}^{\top}$ and ${\bf A}=\begin{bmatrix}{\bf A}_{1}^{\top}&\cdots&{\bf A}_{n_{t}}^{\top}\end{bmatrix}^{\top}.$ The least squares problem above is solved using genHyBR where the regularization parameter is picked using the WGCV criterion to obtain one reconstructed image of size $256\times 256$. Here we set ${\bf Q}_{s}$ to be a covariance matrix that is determined from the Matérn covariance function $C_{S}=C_{1,.01}(\cdot)$ as described in Section 2.3. It is evident that the static reconstruction is able to locate the object but can neither distinguish the objects nor provide dynamic information.

Next we consider three cases that use simultaneous genHyBR to solve (28):

• •

genHyBR where ${\bf Q}$ is generated from a Matérn kernel $C_{1,.01}(\sqrt{c_{1}\|{\bf p}\|^{2}+c_{2}|t|^{2}})$ where $c_{1}=1$ and $c_{2}=0.0025$. Here, ${\bf Q}$ can not be represented as a Kronecker product, but MVPs can still be done efficiently.

• •

genHyBR with ${\bf Q}={\bf Q}_{t}\otimes{\bf Q}_{s}$ where ${\bf Q}_{t}={\bf I}$ and ${\bf Q}_{s}$ corresponds to $C_{S}(\cdot)=C_{1,.01}(\cdot)$.

• •

genHyBR with ${\bf Q}={\bf Q}_{t}\otimes{\bf Q}_{s}$ where ${\bf Q}_{t}$ and ${\bf Q}_{s}$ correspond to $C_{T}(\cdot)=C_{\infty,.01}(\cdot)$ and $C_{S}(\cdot)=C_{1,.01}(\cdot)$ respectively.

Image reconstructions for various time points, along with the corresponding true images, are provided in Figure 2, where all of the results use the WGCV parameter after $10$ iterations. Compared to the static reconstruction in Figure 1(c), the Matérn reconstructions in the second row reveal changes over time. However, the more striking comparison occurs when using the separable covariance functions and comparing the results with ${\bf Q}_{t}={\bf I}$ to those with ${\bf Q}_{t}\neq{\bf I}$ (c.f. rows 3 and 4 in Figure 2). Dynamic PAT is a severely underdetermined problem, and this example illustrates that including a temporal prior can be crucial to revealing dynamics of the imaged object. In Table 1 we provide the computed regularization parameters for each approach, along with the relative errors computed as $\|{\bf s}_{k}-{\bf s}_{\mathrm{true}}\|_{2}/\|{\bf s}_{\mathrm{true}}\|_{2}$. These values are consistent with the quality of the reconstructions in Figure 2 and are not significantly improved with reorthogonalization of gen-GK vectors. For reconstructions that assume ${\bf Q}$ is a Kronecker product, genHyBR took around $37$ seconds and, with reorthogonalization, around $106$ seconds222All timings were recorded on a MacPro, OSX Yosemite, 2.7 GHz 12-Core Intel Xeon E5, 64G memory in Matlab 2014b using default computing options.. Partial reorthogonalization could be used but was not investigated here.

Next, we show that variance estimates (i.e., approximations to diagonals of the posterior covariance matrix) can be obtained with minimal additional costs (here, in 15 seconds). In Figure 3 we provide results for Matérn and ${\bf Q}_{t}\neq{\bf I}$ (corresponding to the MAP estimates in the 2nd and 4th rows of Figure 2), where we note that both approaches provide overall variances on the order of $10^{-4}$. We observe that solutions corresponding to earlier and later time points (e.g., $i=1$ and $i=120$) contain higher variances (i.e., greater uncertainty), with smaller variances in the center regions of the images, especially for ${\bf Q}_{t}\neq{\bf I}.$ Variance images for ${\bf Q}_{t}={\bf I}$ were essentially constant with mean value $0.0257$ and standard deviation $2.3\times 10^{-6}$ and thus are omitted.

In summary, we have shown that the gen-GK bidiagonalization can be used for the efficient computation of MAP estimates and variance estimates for dynamic PAT problems where the underlying object is changing slowly relative to the rate of image acquisition. Various choices for the prior covariance matrices could be included in this framework.

5.2 Space-time image deblurring

In dynamic image deblurring, the goal is to reconstruct a sequence of images from a sequence of blurred and noisy images. We consider a simulated problem where $9$ true images of size $50\times 50$ are shown in Figure 4(a) and the corresponding observed images are shown in Figure 4(b). The blur matrix was taken to be ${\bf A}={\bf A}_{t}\otimes{\bf A}_{s}$ where ${\bf A}_{s}$ represents a 2D Gaussian point spread function with spread parameter $\sigma=.07$ and bandwidth $3$ and ${\bf A}_{t}$ represents a 1D Gaussian blur with spread parameter $\sigma=1$ and bandwidth $3$. The noise level is set to be $0.02$ such that ${\bf R}=0.0437^{2}{\bf I}.$ The problem set-up is a modification of the ‘blur’ example from the Regularization Tools toolbox [21].

We compare LSQR, HyBR opt, genLSQR and genHyBR with two different regularization parameter selection techniques: optimal regularization parameter and the WGCV parameter. Here genLSQR means that $\lambda=0$, and LSQR means that $\lambda=0$ and ${\bf Q}={\bf I}$. For genLSQR and genHyBR, we used ${\bf Q}={\bf Q}_{t}\otimes{\bf Q}_{s}$ where ${\bf Q}_{t}$ and ${\bf Q}_{s}$ correspond to $C_{T}(\cdot)=C_{1.5,.3}(\cdot)$ and $C_{S}(\cdot)=C_{.5,.007}(\cdot)$ respectively. Relative errors per iteration provided in Figure 5 reveal similar behavior as that described in [8]. In particular, LSQR and genLSQR are plagued by semiconvergence (i.e., the “U”-shaped error curve that results from noise contamination during inversion), which can be avoided in the hybrid variants with the selection of the optimal regularization parameter. WGCV is able to provide a fairly good regularization parameter, but the process terminated at iteration $24$ due to a flat GCV curve.

Since ${\bf A}$ also has a Kronecker product structure, the decoupled approach applies here. We computed MAP approximations using the decoupled approach, where genHyBR was used to solve each subproblem (26). We denote ‘decoupled $\lambda_{\rm wgcv}$’ to be the solution using the decoupled approach with a fixed regularization parameter $\lambda_{\rm wgcv}$, and ‘decoupled $\lambda_{\rm wgcv}^{(i)}$’ refers to using a different regularization parameter for each subproblem. WGCV-selected regularization parameters $\lambda_{\rm wgcv}^{(i)}$ and corresponding stopping iterations $k_{\rm stop}^{(i)}$ are provided in Table 2, along with regularization parameters $\lambda_{\rm wgcv}$ and $\lambda_{\rm opt}$. We remark that the regularization parameters $\lambda_{\rm wgcv}^{(i)}$ in decoupled approach decrease with increasing index $i$; this can be attributed to the scaling factor from the singular values and the changing right hand sides.

The relative reconstruction error for the decoupled $\lambda_{\rm wgcv}$ reconstruction was $0.2461$, which is slightly smaller than that of (simultaneous) genHyBR WGCV which was $0.2507$ at termination. Furthermore, we observed that allowing different regularization parameters can result in smaller reconstruction error. The relative error for the decoupled $\lambda_{\rm wgcv}^{(i)}$ reconstruction was $0.2198$. Also, allowing a different Krylov subspace for each reconstruction can be beneficial in reducing “ghosting” errors from neighboring images, as evident in the absolute error images provided in Figure 6, where the same color axis is used per $i$ and black correspond to larger errors. HyBR opt and genHyBR opt are only provided for reference, since they require the optimal regularization parameter that is not available in practice.

Variance estimates for genHyBR WGCV are provided in the first row of Figure 7, and variance estimates in the second row illustrate that the decoupled variance estimate approach with fixed regularization parameter (as explained in Section 4.2) can provide a good approximation.

We note that decoupled $\lambda_{\rm wgcv}^{(i)}$ does not directly fit our framework; however, a modification of ${\bf Q}_{t}^{-1}$ may be used to incorporate the changing regularization parameters. In summary, the decoupled approach can be used for both MAP and variance estimation if ${\bf A}$, ${\bf Q},$ and ${\bf R}$ are all Kronecker products.

As a final remark, in Section 3.2 we assume that the regularization parameter is fixed; however, in this section, we present results for which the regularization parameter is allowed to be different. This certainly has benefits as demonstrated; however, its statistical meaning is not fully clear and is worth exploring in future work.

5.3 Passive seismic tomography (PST)

Recent advances in PST have enabled the monitoring of mining-induced stress redistribution in coal and hardrock mines [31, 32, 53]. The basic goal of microseismic tomography is to image subsurface properties by using the many low-magnitude seismic events (e.g., microearthquakes) that are recorded by a microseismic monitoring system in a deep mine. Using time-lapse PST tomogram reconstructions, we can better understand the stress redistribution within the rock mass so that trends preceding and following significant seismic events can be analyzed. However, obtaining these 3D spatial reconstructions in real-time is a computationally challenging task that consists of solving a sequence of very large, often nonlinear, inverse problems.

In this work, we consider a simplified, linear PST problem in a dynamic framework and investigate gen-GK methods for computing reconstructions. The basic formulation of the problem is the same as (1) where ${\bf s}$ is a discretization of the velocity model, ${\bf d}$ contains the observed travel times or recorded sinogram, and ${\bf A}$ simulates a ray trace operation. We consider the situation in which measurements are taken in periodic time intervals, and the goal is to generate velocity models over time, from which the changing conditions within the rock mass, inferred to be changing due to induced seismicity, can be obtained.

PST simulated data.

As is commonly done in practice, we begin with a simulated problem where the goal is to reconstruct a “checkerboard” image [53]. We create eight checkerboard volumes of size $50\times 66\times 61$ voxels to represent true images, where the values of the checkerboard structure are generated to be reciprocals of $20,000\pm 10\%$ (i.e., values are $4.545\times 10^{-5}$ and $5.555\times 10^{-5}$) in the region of the volume which is seismically “observable.” Cross-sections from the $4$th such generated structures are provided in the top row of Figure 9. Then we used straight-path ray trace matrices ${\bf A}_{i}$ for $i=1,\ldots,8$ from real mine data to generate sinograms. Each matrix corresponds to seismic events that occurred and were detected in a given time period; see Table 3. Observed data were constructed using (1) where ${\boldsymbol{\varepsilon}}$ represented Gaussian noise with ${\bf R}=.0015^{2}\,{\bf I}.$ In summary, the number of unknowns for this problem is $50*66*61*8=1,610,400$ and the total number of observations is $m=191,856$ where the number of observations per time point is provided in Table 3.

In Figure 8, we provide relative reconstruction errors in the observable regions per iteration for various methods:

• •

genLSQR corresponds to $C_{T}(\cdot)=C_{3.5,0.09}(\cdot)$, $C_{S}=C_{2.5,.025}(\cdot)$ and $\lambda=0$.

• •

genHyBR ${\bf Q}_{t}\neq{\bf I}$ corresponds to $C_{T}(\cdot)=C_{3.5,0.09}(\cdot)$, $C_{S}=C_{2.5,.025}(\cdot)$.

• •

genHyBR ${\bf Q}_{t}={\bf I}$ corresponds to $C_{S}(\cdot)=C_{2.5,.025}(\cdot)$.

Here $C_{\nu,\ell}(\cdot)$ corresponds to a covariance matrix determined from a Matérn kernel with parameters $\nu$ and $\ell$ as defined in Section 2.3, and we assume, for simplicity, that ${\bf p}$ and $t$ are equally spaced on a grid that is normalized to $[0,1].$ Regularization parameter $\lambda_{\rm opt}$ was used for genHyBR.

For LSQR and HyBR (defined in Section 5.2), we follow current practice and take an initial guess ${\bf s}_{0}$ to be a constant image with all entries equal to $5\times 10^{-5}$. For the genHyBR reconstructions, we use a physically-informed prior mean and take ${\boldsymbol{\mu}}$ to be a constant image with all entries equal to $5\times 10^{-5}$. We found these inclusions to be critical for obtaining physically meaningful results.

Oftentimes in PST, the temporal prior is ignored (e.g., ${\bf Q}_{t}={\bf I}$) and reconstructions for each time point are done independently (with ${\bf Q}_{s}={\bf I}$). Comparing HyBR and genHyBR for ${\bf Q}_{t}={\bf I}$, we observe that better reconstructions can be obtained by including a spatial prior ${\bf Q}_{s}$. Furthermore, these results show that incorporating a temporal prior may lead to additional improvements. Cross-sections of the $4$th volume are provided in Figure 9 for HyBR, and genHyBR for ${\bf Q}_{t}={\bf I}$ and ${\bf Q}_{t}\neq{\bf I}$.

PST real data.

As demonstrated in the simulated problem, the small ray path coverage makes dynamic PST a highly ill-posed problem. Next, we consider the real field data measurements and reconstruct a time-lapse of eight volumes using genHyBR. The true volumes are unknown; here we only show isosurfaces and comparisons to currently used algorithms. Along with expert field knowledge, this information can aid in evaluating the reconstructions and the potential for future improvements.

We present three reconstructions using genHyBR. The first approach essentially mimics what is done in practice, which is to compute reconstructions independent of time. Here we used genHyBR for each time period with ${\bf Q}_{s}$ corresponding to Matérn covariance function $C_{S}(\cdot)=C_{10.5,.006}(\cdot)$. In the time independent approach, different regularization parameters and stopping iterations were selected for each reconstruction. Then we used simultaneous genHyBR with ${\bf Q}_{t}={\bf I}$ and ${\bf Q}_{t}$ corresponding to $C_{T}(\cdot)=C_{.3,.3}(\cdot)$. In both cases, we used ${\bf Q}_{s}$ as defined above. For all of these experiments, WGCV was used to select the regularization parameter and automatic stopping criteria was used as described in [5] with a maximum of 10 iterations. Obtaining one dynamic PST reconstruction after 10 simultaneous genHyBR iterations took approximately $11$ seconds.

High velocity iso-values (corresponding to a value of $20,025$ or a slowness of $4.9938\times 10^{-5}$) and contours for different time intervals are shown in Figures 10 and 11. Notice that including the temporal prior can result in better reconstructions, especially for time intervals with very few observations (e.g., Jan 25–31 and Feb 1–17). In addition, genHyBR with ${\bf Q}_{t}\neq{\bf I}$ gives more detailed information and locality of high stresses than the time independent reconstructions, where the isosurfaces are more dispersed.

We compare our reconstructions to a typical field reconstruction using the simultaneous iterative reconstruction technique (SIRT) for the time period February 22–28. In Figure 12 we display cross-sections of genHyBR reconstructions, along with the SIRT reconstruction. These correspond to the $24$th, $40$th, and $27$th slices in the $x$, $y$, and $z$ dimensions respectively. We remark that the genHyBR reconstructions all provide smoother, more localized reconstructions of high-velocity zones. It is worth mentioning that these reconstructions have assumed that the forward model has “straight-ray” paths and a typical approach in mining would be to use this reconstruction as an initial guess for obtaining reconstructions with a more sophisticated, nonlinear “curved-ray” forward model. This is another topic of future work.

6 Conclusions

We consider the problem of dynamic inverse problems using the Bayesian framework, and we developed efficient, iterative, matrix-free methods based on the gen-GK bidiagonalization. A wide range of priors can be incorporated in our framework. We focused on priors that are modeled as Gaussian random fields with special attention to space-time covariance kernels for which MVPs can be computed efficiently. We first focus on computing the MAP estimate. In the simultaneous approach, a solution for the entire unknown in space-time is solved in an “all-at-once” manner. When the observation operator also has Kronecker product structure, a series of variable transformations enables the problem to decouple in time. Both the simultaneous and decoupled approaches leverage the efficient iterative solvers developed in our previous work [8], and has the added benefit that the simultaneous approach allows for automatic selection of regularization parameter selection. In addition to the MAP estimate, we describe methods that reuse intermediate information contained in the iterative solvers to estimate the variance of the posterior distribution. Several examples from image processing, including new applications to PST, demonstrate scalability of our algorithms and illustrate the broad applicability of our work.

7 Acknowledgements

We would like to acknowledge Creighton mines for providing the raw data for the PST application. Furthermore, some of this work was conducted as a part of the SAMSI Program on Optimization 2016-2017. This material was based upon work partially supported by the National Science Foundation under Grant DMS-1127914 to the Statistical and Applied Mathematical Sciences Institute.