Randomized Krylov methods for inverse problems

TL;DR

This paper introduces randomized Krylov subspace methods, including a randomized Golub-Kahan approach, for efficiently solving large-scale linear inverse problems with regularization, demonstrated through image deblurring and seismic tomography.

Contribution

It develops new randomized iterative solvers for inverse problems, extending existing methods to rectangular matrices and enabling automatic regularization parameter selection.

Findings

Effective in large-scale inverse problems

Improves computational efficiency with randomized methods

Demonstrated on image deblurring and seismic tomography

Abstract

In this paper we develop randomized Krylov subspace methods for efficiently computing regularized solutions to large-scale linear inverse problems. Building on the recently developed randomized Gram-Schmidt process, where sketched inner products are used to estimate inner products of high-dimensional vectors, we propose a randomized Golub-Kahan approach that works for general rectangular matrices. We describe new iterative solvers based on the randomized Golub-Kahan approach and show how they can be used for solving inverse problems with rectangular matrices, thus extending the capabilities of the recently proposed randomized GMRES method. We also consider hybrid projection methods that combine iterative projection methods, based on both the randomized Arnoldi and randomized Golub-Kahan factorizations, with Tikhonov regularization, where regularization parameters can be selected automatically during the iterative process. Numerical results from image deblurring and seismic tomography show the potential benefits of these approaches.

Tables

Table 1. Table 1 : Summary of the new hybrid randomized methods, as compactly described in ( 44 ).

method	definition	$ℛ\left({𝐕}_k\right)$	${𝐏}_k$	${𝐜}_k$
hybrid rGMRES	(33)	${𝒦}_k(𝐀,𝐛)$	${𝐇}_k^{\top }{𝐇}_k$	${𝐇}_k^{\top }\beta {𝐞}_1$
hybrid rLSQR	(36)-(37)	${𝒦}_k({𝐀}^{\top }𝐀,{𝐀}^{\top }𝐛)$	${𝐌}_k^{\top }{𝐌}_k$	${𝐌}_k^{\top }\beta {𝐞}_1$
hybrid rCGLS	(40)-(41)	${𝒦}_k({𝐀}^{\top }𝐀,{𝐀}^{\top }𝐛)$	${\widetilde{𝐓}}_{k+1}{𝐌}_k$	$\beta t_{1,1}{𝐞}_1$
hybrid rLSMR	(42)-(43)	${𝒦}_k({𝐀}^{\top }𝐀,{𝐀}^{\top }𝐛)$	${\left({𝐓}_{k+1}{𝐌}_k\right)}^{\top }\left({𝐓}_{k+1}{𝐌}_k\right)$	$\beta t_{1,1}{\left({𝐓}_{k+1}{𝐌}_k\right)}^{\top }{𝐞}_1$

Equations

$$\mathbf{b}=\mathbf{A}{\mathbf{x}}_{\text{true}}+\mathbf{e}$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }\Vert \mathbf{A}\mathbf{x}-\mathbf{b}\Vert ,$$

$${\bf x}_{\text{\scriptsize reg}}=\operatorname*{arg\,min}_{{\bf x}\in\mathbb{R}^{n}}\|{\bf A}{\bf x}-{\bf b}\|^{2}+\lambda^{2}\|{\bf x}\|^{2}=\operatorname*{arg\,min}_{{\bf x}\in\mathbb{R}^{n}}\left\|\left[\begin{array}[]{cc}{\bf A}\\ {\lambda}{\bf I}_{n}\end{array}\right]{\bf x}-\left[\begin{array}[]{cc}{\bf b}\\ {\bf 0}\end{array}\right]\right\|^{2},$$

$$⟨\cdot ,\cdot ⟩\approx ⟨{\mathbf{\Theta }}^{(n)}\cdot ,{\mathbf{\Theta }}^{(n)}\cdot ⟩.$$

$$\left|\langle{\bf x},{\bf y}\rangle-\langle{\bm{\Theta}}^{(n)}{\bf x},{\bm{\Theta}}^{(n)}{\bf y}\rangle\right|\leq\varepsilon\|{\bf x}\|\|{\bf y}\|,\quad\mbox{for all ${\bf x},{\bf y}\in\mathcal{V}$}.$$

$$(1+\epsilon {)}^{-1}\Vert {\mathbf{\Theta }}^{(n)}\mathbf{x}{\Vert }^2\le \Vert \mathbf{x}{\Vert }^2\le (1-\epsilon {)}^{-1}\Vert {\mathbf{\Theta }}^{(n)}\mathbf{x}{\Vert }^2.$$

$$\mathbf{A}{\bar{\mathbf{Q}}}_k={\bar{\mathbf{Q}}}_{k+1}{\bar{\mathbf{H}}}_k.$$

$${\bar{\mathbf{x}}}_k={\bar{\mathbf{Q}}}_k{\bar{\mathbf{z}}}_k,\text{mbox}where{\bar{\mathbf{z}}}_k=\underset{\mathbf{z}}{\mathrm{argmin}}\Vert {\bar{\mathbf{H}}}_k\mathbf{z}-\bar{\beta }{\mathbf{e}}_1\Vert \text{mbox}with\bar{\beta }=\Vert \mathbf{b}\Vert .$$

$${\mathbf{r}}_{k-1}=(\underset{=:{\mathbf{S}}_{k-1}}{\underbrace{{\mathbf{\Theta }}^{(n)}{\mathbf{Q}}_{k-1}}}{)}^{\top }({\mathbf{\Theta }}^{(n)}{\mathbf{w}}_k),$$

$${\mathbf{q}}_k=\text{nicefrac}({\mathbf{w}}_k-{\mathbf{Q}}_{k-1}{\mathbf{r}}_{k-1})\Vert {\mathbf{\Theta }}^{(n)}({\mathbf{w}}_k-{\mathbf{Q}}_{k-1}{\mathbf{r}}_{k-1})\Vert .$$

$${\mathbf{r}}_{k-1}=\underset{\mathbf{y}\in {\mathbb{R}}^{k-1}}{\mathrm{argmin}}\Vert {\mathbf{\Theta }}^{(n)}{\mathbf{Q}}_{k-1}\mathbf{y}-{\mathbf{\Theta }}^{(n)}{\mathbf{w}}_k\Vert .$$

$$\mathbf{A}{\mathbf{Q}}_k={\mathbf{Q}}_{k+1}{\mathbf{H}}_k$$

$${\mathbf{x}}_k={\mathbf{Q}}_k{\mathbf{z}}_k,\text{mbox}where\text{mbox}{\mathbf{z}}_k=\underset{\mathbf{z}}{\mathrm{argmin}}\Vert {\mathbf{H}}_k\mathbf{z}-\beta {\mathbf{e}}_1\Vert \text{mbox}and\beta =\Vert {\mathbf{\Theta }}^{(n)}\mathbf{b}\Vert .$$

$$\begin{array}{rl}\Vert {\mathbf{H}}_k{\mathbf{z}}_k-\beta {\mathbf{e}}_1\Vert & =\underset{\mathbf{z}}{\min }\Vert {\mathbf{H}}_k\mathbf{z}-\beta {\mathbf{e}}_1\Vert =\underset{\mathbf{z}}{\min }\Vert {\mathbf{\Theta }}^{(n)}{\mathbf{Q}}_{k+1}({\mathbf{H}}_k\mathbf{z}-\beta {\mathbf{e}}_1)\Vert \\ & =\underset{\mathbf{z}}{\min }\Vert {\mathbf{\Theta }}^{(n)}(\mathbf{A}{\mathbf{Q}}_k\mathbf{z}-\mathbf{b})\Vert =\Vert {\mathbf{\Theta }}^{(n)}(\mathbf{A}{\mathbf{x}}_k-\mathbf{b})\Vert .\end{array}$$

$${\mathbf{x}}_k\in {\mathcal{K}}_k(\mathbf{A},\mathbf{b})\text{mbox}and{\mathbf{r}}_k{\perp }_{{\mathbf{\Theta }}^{(n)}}\mathbf{A}{\mathcal{K}}_k(\mathbf{A},\mathbf{b}),$$

$$\Vert \mathbf{b}-\mathbf{A}{\mathbf{x}}_k{\Vert }^2\le \frac{(1+\epsilon )}{(1-\epsilon )}\Vert \mathbf{b}-\mathbf{A}{\bar{\mathbf{x}}}_k{\Vert }^2,$$

$$(1-\epsilon )\Vert \mathbf{b}-\mathbf{A}{\mathbf{x}}_k{\Vert }^2\le \Vert {\mathbf{\Theta }}^{(n)}(\mathbf{b}-\mathbf{A}{\mathbf{x}}_k){\Vert }^2\le \Vert {\mathbf{\Theta }}^{(n)}(\mathbf{b}-\mathbf{A}{\bar{\mathbf{x}}}_k){\Vert }^2\le (1+\epsilon )\Vert \mathbf{b}-\mathbf{A}{\bar{\mathbf{x}}}_k{\Vert }^2,$$

$${\mathbf{z}}_k=\underset{\mathbf{z}}{\mathrm{argmin}}\Vert {\mathbf{H}}_k\mathbf{z}-\beta {\mathbf{e}}_1\Vert ,\text{mbox}where{\mathbf{H}}_k=[\mathbf{R}{]}_{(1:k+1,2:k+1)}.$$

$$\mathcal{R}({\bar{\mathbf{V}}}_k)={\mathcal{K}}_k({\mathbf{A}}^{\top }\mathbf{A},{\mathbf{A}}^{\top }\mathbf{b})\subset {\mathbb{R}}^n\text{mbox}and\mathcal{R}({\bar{\mathbf{U}}}_{k+1})={\mathcal{K}}_{k+1}(\mathbf{A}{\mathbf{A}}^{\top },\mathbf{b})\subset {\mathbb{R}}^m,$$

$$\mathbf{A}{\bar{\mathbf{V}}}_k={\bar{\mathbf{U}}}_{k+1}{\bar{\mathbf{B}}}_k\text{mbox}and{\mathbf{A}}^{\top }{\bar{\mathbf{U}}}_{k+1}={\bar{\mathbf{V}}}_{k+1}({\bar{\mathbf{B}}}_{k+1,k}{)}^{\top },$$

$$⟨\cdot ,\cdot ⟩\approx ⟨{\mathbf{\Theta }}^{(n)}\cdot ,{\mathbf{\Theta }}^{(n)}\cdot ⟩\text{mbox}and⟨\cdot ,\cdot ⟩\approx ⟨{\mathbf{\Theta }}^{(m)}\cdot ,{\mathbf{\Theta }}^{(m)}\cdot ⟩,$$

$$\bar{\bf q}_{j}^{\top}\widehat{{\bf A}}\bar{\bf q}_{i}=h_{j,i}\mbox{ ($=0$ if $j>i+1$)}\quad\mbox{and}\quad\bar{\bf q}_{i}^{\top}\widehat{{\bf A}}\bar{\bf q}_{j}=h_{i,j}\mbox{ ($=0$ if $i>j+1$)}\,.$$

$${\mathbf{q}}_i^{\top }\hat{\mathbf{A}}{\mathbf{q}}_j\ne {\mathbf{q}}_j^{\top }\hat{\mathbf{A}}{\mathbf{q}}_i={\mathbf{q}}_j^{\top }{\mathbf{Q}}_{k+1}{\mathbf{H}}_k{\mathbf{e}}_i\ne 0\text{for i,j=1,…,k,}$$

$${\mathbf{Q}}_{k+1}^{\top }(({\mathbf{\Theta }}^{(n)}{)}^{\top }){\mathbf{\Theta }}^{(n)}\hat{\mathbf{A}}{\mathbf{Q}}_k={\mathbf{H}}_k,$$

$$\mathbf{A}{\mathbf{V}}_k={\mathbf{U}}_{k+1}{\mathbf{M}}_k,\text{mbox}and{\mathbf{A}}^{\top }{\mathbf{U}}_{k+1}={\mathbf{V}}_{k+1}{\mathbf{T}}_{k+1},$$

$${\mathbf{P}}_k=\left[\begin{array}{cc}{\mathbf{p}}_1& \dots {\mathbf{p}}_k\end{array}\right]\in {\mathbb{R}}^{{\ell }_n\times k}\text{mbox}and{\mathbf{S}}_k=\left[\begin{array}{cc}{\mathbf{s}}_1& \dots {\mathbf{s}}_k\end{array}\right]\in {\mathbb{R}}^{{\ell }_m\times k}.$$

$$[\mathbf{M}{]}_{(1:k,k)}=\underset{\mathbf{y}}{\mathrm{argmin}}\Vert {\mathbf{S}}_k\mathbf{y}-{\tilde{\mathbf{s}}}_{k+1}\Vert \text{mbox}and[\mathbf{T}{]}_{(1:k,k+1)}=\underset{\mathbf{y}}{\mathrm{argmin}}\Vert {\mathbf{P}}_k\mathbf{y}-{\tilde{\mathbf{p}}}_{k+1}\Vert ,$$

$${\tilde{\mathbf{s}}}_{k+1}={\mathbf{\Theta }}^{(m)}{\tilde{\mathbf{u}}}_{k+1}\text{mbox}and{\tilde{\mathbf{p}}}_{k+1}={\mathbf{\Theta }}^{(n)}{\tilde{\mathbf{v}}}_{k+1},$$

$$\begin{array}{rl}(2K+1)C_{\text{mv}}& +(K+1)(m+n+C_{{\text{sk}}_n}+C_{{\text{sk}}_m}+3({\ell }_n+{\ell }_m)-2)\\ & +K(K+1)(m+n+2{\ell }_m+2{\ell }_n-1)\end{array}$$

$$(2K+1)C_{\text{mv}}+(K+1)(5(m+n)-2)-2(m+n)$$

Full text

{MSCcodes}

65F22 65F10 65K10 15A29

1 Introduction

Randomized approaches, in particular randomization for obtaining low-rank approximations or solving linear least-squares problems, have emerged as powerful tools in the numerical linear algebra community; see [25] and references therein. Randomized Krylov methods are iterative methods that incorporate randomization in various ways: their main benefit is the potential of dramatically reducing computational costs with respect to their deterministic counterparts, while maintaining strong theoretical justifications and properties. The broad goal of this work is to develop randomized Krylov methods for large-scale linear inverse problems, which may be used as purely iterative regularization methods or in a hybrid fashion, i.e., combining iterative and variational regularization.

Let us consider linear inverse problems of the form,

[TABLE]

where ${\bf b}\in\mathbb{R}^{m}$ collects observed measurements, ${\bf A}\in\mathbb{R}^{m\times n}$ represents the discretized forward model, ${\bf x}_{\text{\scriptsize true}}\in\mathbb{R}^{n}$ is an unknown quantity of interest, and ${\bf e}\in\mathbb{R}^{m}$ contains noise or errors in the data. Given ${\bf b}$ and ${\bf A},$ the goal of the inverse problem is to approximate ${\bf x}_{\text{\scriptsize true}}$. This is typically achieved by applying some kind of regularization to the original linear system (1). Since, in general, ${\bf A}$ is large and unstructured, a common approach to recover ${\bf x}_{\text{\scriptsize true}}$ is to apply an iterative (e.g., standard Krylov) solver to the least-squares problem

[TABLE]

where $\|\cdot\|$ denotes the vector 2-norm, sometimes also referred to as the $l_{2}$ norm. Early termination of many iterative solvers for (2) results in a regularized solution; see [20]. Another common approach to approximate ${\bf x}_{\text{\scriptsize true}}$ is to solve a Tikhonov-regularized problem, i.e., compute

[TABLE]

where $\lambda>0$ is the regularization parameter. Many iterative (e.g., standard Krylov) methods can be used to solve minimization problems Eq. 3 efficiently, but they can get computationally expensive for large-scale problems that require many iterations. Moreover, the regularization parameter $\lambda$ is typically not known a priori, and estimating it can be costly. By combining Tikhonov regularization (3) and iterative projection methods, hybrid projection methods provide a robust and computationally efficient approach for simultaneous solution approximation and regularization parameter selection; see [10] for a survey paper with more details.

Randomized methods have been considered for solving linear least-squares problems like Eqs. 2 and 3. In this setting there are various ways to incorporate randomization, the most common ones being the so-called “sketch-and-solve” and “sketch-and-precondition” paradigms. In the first instance and considering an overdetermined matrix ${\bf A}$, sketching is initially applied to reduce the dimensionality of the matrix ${\bf A}$ and the data vector ${\bf b}$; a standard algorithm is then used to solve the approximate problem. Such approaches can provide low-rank matrix approximations, but they are not suited for inverse problems with a slowly decaying spectrum. In the context of inverse problems, a randomized approach was proposed in [24], where the inner product free LSLU projection method was combined with “sketch-and-solve” on the projected problem. In the second instance, a sketching matrix is used to obtain a preconditioner for the considered iterative methods: for example, this underlies the LSNR [28] and the Blendenpik [2] algorithms. There is a vast body of literature (often within the statistical and data science remit), which explores randomization to solve regularized linear least-squares problems, formulated as damped least-squares problems; see, for instance, [21, 15, 8, 26]. Although these sketching-based strategies have gained widespread interest as tools for low-rank matrix approximation and for solving linear systems, they are numerically unstable for ill-conditioned, and even moderately ill-conditioned, problems; see [27]. Moreover, their utility for state-of-the-art inverse problems remains largely under-explored. Apart from the efficient solution of linear least-squares problems, other previous works on randomized Krylov methods focus on block Krylov methods for spectral computation and matrix approximation (e.g., for approximating extreme singular values and singular vectors); see, e.g., [29, 34]. More recently, randomized iterative solvers have been developed, where sketching techniques are used to reduce the computational cost of (re)orthogonalization [4], and have been used to speed up approximating matrix functions [13]. Specifically, a randomized Gram-Schmidt process was described in [4], which forms the basis for the randomized Arnoldi (rArnoldi) algorithm and the randomized GMRES (rGMRES) method. These approaches are restricted to square matrices.

In this work, we develop randomized Krylov methods for large-scale inverse problems where randomization is incorporated via sketched inner products. The main contributions of this work are the following.

• •

Extending the approach of [4], we develop a randomized Golub-Kahan (rGK) algorithm that can be applied to a rectangular matrix. This approach requires two different sketching matrices and generates an upper Hessenberg and an upper triangular matrix. We provide cost and storage comparisons to (standard) Golub-Kahan bidiagonalization.

• •

Then, building on the proposed rGK approach, we introduce various randomized iterative projection methods, where the main difference among the methods is the constraint imposed on the residuals or normal equations residuals. Specifically, we describe a randomized LSQR (rLSQR) method, a randomized CGLS (rCGLS) method, and a randomized LSMR (rLSMR) method.

• •

For all of the considered randomized Krylov methods (rGMRES for square systems, and rLSQR, rCGLS, and rLSMR for both square and rectangular systems), we describe their hybrid versions, where we estimate the regularization parameter during the iterative process by applying well-established regularization parameter selection methods.

An outline for this paper is as follows. In Section 2, we review the rArnoldi and rGMRES algorithms for square matrices, and we propose the rGK algorithm for rectangular matrices. We describe various randomized iterative methods for the linear least-squares problem (2) based on the rGK factorization, and draw connections to their standard counterparts. Hybrid projection methods based on these randomized projection methods and Tikhonov regularization are described in Section 3. Numerical results are provided in Section 4, and conclusions and future work are described in Section 5.

2 Randomized Krylov methods

Krylov subspace methods have been widely studied and used for solving large-scale inverse problems, where each iteration requires a matrix-vector product with ${\bf A}$ and, possibly, ${\bf A}^{\top}$. Various works have reduced the cost of matrix-vector products, e.g., by exploiting structure or sparsity in ${\bf A}$ or utilizing highly parallel distributed architectures; see [30, 35] just to cite a few. However, a remaining challenge is the computational cost of orthonormalization. In particular, computing inner products can be expensive if the vectors are very high-dimensional or the number of iterations is high. Moreover, inner products can present a communication bottleneck if vectors are distributed across multiple processors. Modifications and variants such as inner product free iterative methods [5, 6] have been considered, but often lead to quasi-minimum residual methods and can require expensive pivoting strategies.

Randomized approaches have been proposed within Krylov methods to reduce the cost of orthonormalization; see, e.g., [4, 31]. Most of these methods rely on the following random sketching technique to estimate inner products. Let $\mathcal{V}$ be a subspace of $\mathbb{R}^{n}$ of low dimension, and let ${\bm{\Theta}}^{(n)}\in\mathbb{R}^{\ell_{n}\times n}$, with $\ell_{n}\ll n$ be a sketching matrix. The $l_{2}$ inner product between vectors in subspaces of $\mathbb{R}^{n}$ is estimated by

[TABLE]

More precisely, we say that ${\bm{\Theta}}^{(n)}$ is an $\varepsilon$ embedding for the subspace $\mathcal{V}$ if

[TABLE]

The above relation implies that

[TABLE]

Adopting a randomized approach allows for savings in that inner products of high dimensional vectors are replaced by inner products of reduced vectors obtained by sketching. This is of course efficient as far as the random sketching matrix can be cheaply applied to a vector. For instance, ${\bm{\Theta}}^{(n)}$ can be chosen as a Gaussian matrix, although applying it to a vector costs $O(\ell_{n}n)$ flops; when ${\bm{\Theta}}^{(n)}$ is chosen as the more sophisticated Subsampled Randomized Hadamard Transform (SRHT), applying it to a vector costs as little as $2n\log_{2}(\ell_{n}+1)$ flops, depending on the implementation and the computational architecture; see [1]. We refer to [36] for a survey of other possible sketching matrices. Since, in the remaining part of the paper, the subspace $\mathcal{V}$ is a Krylov subspace generated as the iterations of the considered solver unfold, we need to consider sketching matrices that do not require a-priori knowledge of $\mathcal{V}$ except for its dimension. More precisely, we consider $(\varepsilon,\delta,K+1)$-oblivious embeddings such that (5) is satisfied with probability at least $1-\delta$ for any fixed $(K+1)$-dimensional subspace of $\mathbb{R}^{n}$. Both Gaussian and SRHT matrices have this property for big enough $\ell_{n}$.

To establish notation and as motivation, we begin in Section 2.1 with a review of the randomized Gram-Schmidt process that was described in [4], and the resulting rArnoldi and rGMRES algorithms. Then, in Section 2.2, we propose a new randomized Golub-Kahan (rGK) approach for rectangular matrices, where two different sketching matrices are used to efficiently project the problem using two different Krylov subspaces. Various randomized Krylov solvers based on the rGK approach are described in Section 2.3.

2.1 rArnoldi and rGMRES

Recall that standard GMRES and the underlying Arnoldi algorithm can be used to solve high-dimensional inverse problems of the form Eq. 1, where ${\bf A}\in\mathbb{R}^{n\times n}$ and ${\bf b}\in\mathbb{R}^{n}$; see, for instance, [22]. By taking the first basis vector $\bar{{\bf q}}_{1}$ as the normalized right-hand-side vector ${\bf b}$, $k\leq K$ iterations of the Arnoldi algorithm generate a matrix with orthonormal columns $\bar{\bf Q}_{k+1}=[\bar{\bf q}_{1},\dots,\bar{\bf q}_{k+1}]\in\mathbb{R}^{n\times(k+1)}$ and an upper Hessenberg matrix $\bar{\bf H}_{k}\in\mathbb{R}^{(k+1)\times k}$ that satisfy the Arnoldi identity

[TABLE]

The columns $\bar{{\bf q}}_{i}$, $i=1,\dots,k+1$, of $\bar{\bf Q}_{k+1}$ span the Krylov subspace $\mathcal{K}_{k+1}({\bf A},{\bf b})={\rm span}\{{\bf b},{\bf A}{\bf b},\ldots,{\bf A}^{k}{\bf b}\}$, i.e., $\mathcal{R}(\bar{\bf Q}_{k+1})=\mathcal{K}_{k+1}({\bf A},{\bf b})$. Here and in the following, $\mathcal{R}(\cdot)$ denotes the range or column space of a matrix. At the $k$th iteration of the GMRES method, the solution to Eq. 2 can be approximated by

[TABLE]

Here and in the following, ${\bf e}_{i}$ denotes the $i$th canonical basis vector, whose size should be clear from the context. The Arnoldi algorithm can be regarded as a specific instance of the Gram-Schmidt process, which is used to construct an orthonormal basis for $\mathcal{K}_{k}({\bf A},{\bf b})$ in an iterative process. Note that the orthonormalization process can be costly for very large dimensional vectors.

One remedy to mitigate the computational cost of the Gram-Schmidt process is to use the randomized Gram-Schmidt process that exploits sketching, as described in [4]. Given a matrix ${\bf W}\in\mathbb{R}^{n\times K}$ and an $\ell_{n}\times n$ sketching matrix ${\bm{\Theta}}^{(n)}$, with $\ell_{n}\ll n$, two matrices ${\bf Q}\in\mathbb{R}^{n\times K}$ and ${\bf R}\in\mathbb{R}^{K\times K}$ are generated, where the columns of ${\bf Q}$ are orthogonal with respect to the sketched inner product $\langle{\bm{\Theta}}^{(n)}\cdot,{\bm{\Theta}}^{(n)}\cdot\rangle$ and ${\bf R}$ is an upper triangular matrix. More specifically, at the $k$th step, $2\leq k\leq K$, we compute the strictly upper triangular entries ${\bf r}_{k-1}$ of ${\bf R}$ as

[TABLE]

where ${\bf w}_{k}$ is the $k$th column of ${\bf W}$ and ${\bf Q}_{k-1}\in\mathbb{R}^{n\times(k-1)}$ contains the first $(k-1)$ columns of ${\bf Q}$. This is followed by an update for ${\bf q}_{k}$ of the form

[TABLE]

In (9) we use the notation ${\bf S}_{k-1}:={\bm{\Theta}}^{(n)}{\bf Q}_{k-1}\in\mathbb{R}^{\ell_{n}\times{(k-1)}}$ to emphasize the reduction in the size of the vectors that are orthonormalized. Note that, in [4], (9) is replaced by the more general least-squares problem,

[TABLE]

Indeed, the vector ${\bf r}_{k-1}$ in (9) is a specific solution expression corresponding to the classical Gram-Schmidt implementation. Although more sophisticated ways of computing ${\bf r}_{k-1}$ (such as stable solvers based on Householder transformations or Richardson iterations for the associated normal equations) allow for greater stability of the randomized Gram-Schmidt process, we did not observe substantial differences when testing on a variety of inverse problems, so here we use the simpler version (9).

Now, the randomized Gram-Schmidt process can be used to compute an approximation of ${\bf x}_{\text{\scriptsize true}}$ in (1) in the subspaces $\mathcal{K}_{k}({\bf A},{\bf b})$, $k\leq K$. Let ${\bm{\Theta}}^{(n)}$ be an $(\varepsilon,\delta,K+1)$ oblivious subspace embedding (implying that it is an $\varepsilon$-embedding for any $\mathcal{K}_{k}({\bf A},{\bf b})$, $k\leq K+1$). The randomized Arnoldi (rArnoldi) approach is obtained by applying the randomized Gram-Schmidt algorithm to the Krylov subspace $\mathcal{K}_{k}({\bf A},{\bf b})$. More precisely, after $k$ iterations of rArnoldi, we have a basis $\{{\bf q}_{1},\dots,{\bf q}_{k+1}\}$ of $\mathcal{K}_{k+1}({\bf A},{\bf b})$ with columns ${\bf q}_{i}$, $i=1,\dots,k+1$, that are orthonormal with respect to the sketched inner product $\langle{\bm{\Theta}}^{(n)}\cdot,{\bm{\Theta}}^{(n)}\cdot\rangle$, and an upper Hessenberg matrix ${\bf H}_{k}\in\mathbb{R}^{(k+1)\times k}$. Let ${\bf Q}_{k}=[{\bf q}_{1},\dots,{\bf q}_{k}]\in\mathbb{R}^{n\times k}$ and ${\bf Q}_{k+1}=[{\bf Q}_{k},{\bf q}_{k+1}]\in\mathbb{R}^{n\times(k+1)}$, then we get an Arnoldi identity analogous to Eq. 7, i.e.,

[TABLE]

in infinite precision arithmetic. Building on the rArnoldi algorithm, the rGMRES solution at the $k$th iteration is given by

[TABLE]

In infinite precision, the $k$th rGMRES approximation ${\bf x}_{k}$ Eq. 11 minimizes the sketched (semi)norm of the residual over $\mathcal{K}_{k}({\bf A},{\bf b})$, since

[TABLE]

This optimality property of rGMRES can be also derived from more general facts about projection methods (as described, for instance, in [32, Chapter 5]): indeed, the $k$th rGMRES approximation is such that

[TABLE]

where, by $\perp_{{\bm{\Theta}}^{(n)}}$, we mean orthogonal with respect to the sketched inner product $\langle{\bm{\Theta}}^{(n)}\cdot,{\bm{\Theta}}^{(n)}\cdot\rangle$. Note that the above properties generalize the ones used by standard GMRES, which would hold with ${\bm{\Theta}}^{(n)}={\bf I}_{n}$. Using the fact that ${\bm{\Theta}}^{(n)}$ is an $(\varepsilon,\delta,K+1)$-oblivious embedding, we can prove the following bound for the rGMRES residual ${\bf b}-{\bf A}{\bf x}_{k}$ norm with respect to the GMRES residual ${\bf b}-{\bf A}\bar{\bf x}_{k}$ norm at the $k$th iteration

[TABLE]

which comes directly from the inequalities

[TABLE]

where we have used (6) for the first and third inequality, and the rGMRES optimality property (12) for the second inequality.

The randomized GMRES algorithm is provided in Algorithm 1. Note that, with respect to the version proposed in [4] and in addition to the simplified orthogonalization procedure in step 9, step 11 replaces the computation of $\widetilde{{\bf s}}_{k+1}$ with $\widetilde{{\bf s}}_{k+1}={\bm{\Theta}}^{(n)}\widetilde{{\bf q}}_{k+1}$. It can be easily shown by strong induction that the two expressions are equivalent. If applying sketching to an $n$-dimensional vector costs more than $\ell_{n}(k+1)$ flops, then our implementation is more efficient than the one in [4].

2.2 Randomized Golub-Kahan Approach

The aim of this section is to introduce a randomized Golub-Kahan (rGK) algorithm that works for general rectangular matrices ${\bf A}$. First recall that, given ${\bf A}\in\mathbb{R}^{m\times n}$ and ${\bf b}\in\mathbb{R}^{m}$, $k$ iterations of the standard Golub-Kahan bidiagonalization (GKB) algorithm generate two matrices $\bar{\bf U}_{k+1}\in\mathbb{R}^{m\times(k+1)}$ and $\bar{\bf V}_{k}\in\mathbb{R}^{n\times k}$ with numerically orthonormal columns such that

[TABLE]

i.e., the columns of $\bar{\bf V}_{k}$ and $\bar{\bf U}_{k+1}$ span the Krylov subspaces detailed above. The matrices $\bar{\bf V}_{k}$ and $\bar{\bf U}_{k+1}$ are linked by the relations

[TABLE]

where $\bar{\bf B}_{k}\in\mathbb{R}^{(k+1)\times k}$ is lower bidiagonal and $\bar{\bf B}_{k+1,k}$ is obtained by removing the last column of $\bar{\bf B}_{k+1}$.

Similar to the randomized approaches described in Section 2.1, the rGK method is obtained by replacing standard inner products appearing in GKB with sketched inner products. Since, in general, $m\neq n$, we need to consider two different random sketchings for estimating the inner products in the two Krylov subspaces in (13). Assume that we have identified a maximum number of iterations $K$ such that $k\leq K$. Let ${\bm{\Theta}}^{(n)}\in\mathbb{R}^{\ell_{n}\times n}$, with $\ell_{n}\ll n$, be an $(\varepsilon,\delta,K+1)$-oblivious sketching matrix for vectors in $\mathbb{R}^{n}$, and let ${\bm{\Theta}}^{(m)}\in\mathbb{R}^{\ell_{m}\times m}$, with $\ell_{m}\ll m$, be an $(\varepsilon,\delta,K+1)$-oblivious sketching matrix for vectors in $\mathbb{R}^{m}$. That is, the $l_{2}$ inner product between vectors in subspaces of dimensions up to $K$ in $\mathbb{R}^{n}$ and $\mathbb{R}^{m}$ are estimated by

[TABLE]

respectively. If $m=n$, one can simply take ${\bm{\Theta}}^{(n)}={\bm{\Theta}}^{(m)}$.

A first observation is that, while the traditional GKB algorithm enjoys short recurrences, its randomized version will require long (i.e., full Gram-Schmidt-like) recurrences. To understand why this is the case, in a simpler setting, one can consider generating an orthonormal basis for the subspace $\mathcal{K}_{k}(\widehat{{\bf A}},{\bf b})$, where $\widehat{{\bf A}}\in\mathbb{R}^{n\times n}$ is symmetric and ${\bf b}\in\mathbb{R}^{n}$. It is well known that, in this setting, the Arnoldi algorithm reduces to the symmetric Lanczos algorithm, whereby the leading $k\times k$ minor of the upper Hessenberg matrix $\bar{{\bf H}}_{k}$ is a symmetric tridiagonal matrix. One way of motivating why this happens is that, directly from (7) (with ${\bf A}=\widehat{{\bf A}}$),

[TABLE]

Since $\widehat{{\bf A}}$ is symmetric, $\bar{\bf q}_{j}^{\top}\widehat{{\bf A}}\bar{\bf q}_{i}=\bar{\bf q}_{i}^{\top}\widehat{{\bf A}}\bar{\bf q}_{j}$ and, therefore, $h_{i,j}=h_{j,i}=0$ unless $i=j-1,j,j+1$. In a randomized setting, when generating a ${\bm{\Theta}}^{(n)}$-orthonormal basis for $\mathcal{K}_{k}(\widehat{{\bf A}},{\bf b})$ we still get a relation of the kind (7), but relations of the kind (16) do not hold anymore because, due to the lack of orthogonality of the columns of ${\bf Q}_{k+1}=[{\bf q}_{1},\dots,{\bf q}_{k+1}]\in\mathbb{R}^{n\times(k+1)}$,

[TABLE]

in general. Indeed, in this setting, one has a relation of the kind

[TABLE]

so that the leading $k\times k$ minor of ${{\bf H}}_{k}$ is not symmetric unless the matrices $({\bm{\Theta}}^{(n)})^{\top}{\bm{\Theta}}^{(n)}$ and $\widehat{{\bf A}}$ commute (which is very unlikely).

After $k\leq K$ iterations of the rGK algorithm, we have the following two partial matrix factorizations

[TABLE]

where

• •

${\bf V}_{k}=[{\bf v}_{1},\dots,{\bf v}_{k}]\in\mathbb{R}^{n\times k}$ and ${\bf V}_{k+1}=[{\bf V}_{k},{\bf v}_{k+1}]\in\mathbb{R}^{n\times(k+1)}$ have ${\bm{\Theta}}^{(n)}$-orthogonal columns

• •

${\bf U}_{k+1}=[{\bf u}_{1},\dots,{\bf u}_{k+1}]\in\mathbb{R}^{m\times(k+1)}$ has ${\bm{\Theta}}^{(m)}$-orthogonal columns

• •

${\bf M}_{k}\in\mathbb{R}^{(k+1)\times k}$ is an upper Hessenberg matrix

• •

${\bf T}_{k+1}\in\mathbb{R}^{(k+1)\times(k+1)}$ is an upper triangular matrix

and the sketched vectors are the (orthonormal) columns of the matrices

[TABLE]

Algorithm 2 summarizes the main steps involved in the new rGK algorithm.

Note that, similarly to what was done for rArnoldi in Algorithm 1, steps 13 and 21 of Algorithm 2 can be replaced by the more general tasks of approximating the solution to the least-squares problems,

[TABLE]

respectively. Also, depending on the cost of performing sketching, one may replace steps 15 and 23 of Algorithm 2 by the sketched vectors

[TABLE]

respectively.

We now present an estimate of the cost, in terms of flops, of computing the rGK factorization (17), and compare it to the classical GKB. Let us denote by $C_{\text{\small{mv}}}$ the cost of matrix-vector products with ${\bf A}$ and let us assume this equals the cost of matrix-vector products with ${\bf A}^{\top}$; let us denote by $C_{{\text{\small{sk}}}_{n}}$ and $C_{{\text{\small{sk}}}_{m}}$ the cost of performing matrix-vector products with the sketching matrices ${\bm{\Theta}}^{(n)}$ and ${\bm{\Theta}}^{(m)}$, respectively. Then the cost of ($K$ iterations of) Algorithm 2 is estimated to be

[TABLE]

flops. In contrast, the cost of ($K$ iterations of) the (standard) GKB algorithm is

[TABLE]

flops. If GKB is implemented with reorthogonalization (ro-GKB), then its computational cost is

[TABLE]

flops. A few remarks are in order. All the above estimates have the term $(2K+1)C_{\text{\small{mv}}}$ in common, accounting for the matrix-vector multiplications with ${\bf A}$ and ${\bf A}^{\top}$ performed to expand the underlying Krylov subspaces (up to dimension $(K+1)$). The terms of the kind $O(K(m+n))$ account for rescaling and sums of (full-dimensional) vectors of length $m$ and $n$. In rGK and ro-GKB, the terms of the kind $O(K^{2}(m+n))$ account for the computations of the bases of the underlying Krylov subspaces (up to dimension $(K+1)$), and they appear with different scaling factors; in particular, the count for ro-GKB also includes scalar products computations. In rGK, the terms of the kind $O(K(\ell_{n}+\ell_{m}))$, and $O(K^{2}(\ell_{n}+\ell_{m}))$, account for rescaling and sums, and inner products, of sketched vectors, respectively. More precisely, when it comes to the cost of performing inner products, rGK requires

[TABLE]

ro-GKB requires

[TABLE]

standard GKB requires the computations of two norms at each iteration, totaling

[TABLE]

Therefore, in terms of inner product computations and recalling that $\ell_{m}\ll m$ and $\ell_{n}\ll n$, rGK is more convenient than ro-GKB when $C_{{\text{\small{sk}}}_{m}}$ and $C_{{\text{\small{sk}}}_{n}}$ are negligible.

To support the above reasoning about computational costs, we consider $m=n=10^{4}$ and we take ${\bm{\Theta}}^{(m)}={\bm{\Theta}}^{(n)}$ as a SRHT (also mentioned at the beginning of Section 2). It is well-known (see, e.g., [4]) that SRHT is an $(\varepsilon,\delta,K+1)$ oblivious $l_{2}$-subspace embedding if

[TABLE]

In Figure 1 we compare the estimated flops required to run rGK and ro-GKB when $\ell_{n}$ is selected according to (20), with $\varepsilon=\delta=0.5$, as well as with the more heuristic choice

[TABLE]

and when 20 approximately logarithmically equispaced values of $K$ between 10 and 5000 are considered. We also run comparisons when $K$ is fixed and 20 approximately logarithmically equispaced values of $\ell_{n}$ between $K$ and 5000 are selected. We can see that, when $\ell_{n}$ is kept relatively low, there is a clear computational gain in performing rGK rather than ro-GKB.

We conclude this subsection by mentioning that, when regularizing inverse problems, both in a variational and Bayesian setting, it is known that ro-GKB may perform better than GKB, as the former may deliver more accurate estimations of the regularized solution itself or of the so-called hyperparameters defining it (like the Tikhonov regularization parameter $\lambda$ in (3)), often requiring less iterations than standard GKB; see [3, 7, 33] for more details. Therefore, in this setting, one should adopt ro-GKB, and the new rGK may provide clear computational advantages over it.

2.3 Randomized rGK-based solvers for Eq. 2

Building upon the rGK algorithm, a number of solvers can be introduced to approximate a solution to problem Eq. 2. The solvers vary based on the conditions imposed on the residual.

We start by deriving a new randomized LSQR (rLSQR) method that, at the $k$th iteration, defines the $k$th approximate solution as

[TABLE]

and where ${\bf V}_{k}$ and ${\bf M}_{k}$ are the outputs of $k$ iterations of Algorithm 2, where also $\beta(=\|{\bm{\Theta}}^{(m)}{\bf b}\|)$ is defined. In exact arithmetic, it can be easily seen that ${\bf x}_{k}$ minimizes the norm of the sketched residual over the subspace $\mathcal{K}_{k}({\bf A}^{\top}{\bf A},{\bf A}^{\top}{\bf b})$. Indeed,

[TABLE]

where, in the second equality, we have exploited the orthonormality of the columns of ${\bm{\Theta}}^{(m)}{\bf U}_{k+1}$, and, in the third equality, we have used the first of the factorizations in (17). These derivations are analogous to the ones performed to define rGMRES starting from the rArnoldi process (see (12)); the optimality property for these two methods is analogous, too. In the general framework for iterative projection methods for (2) (see again [32, Chapter 5]), rLSQR can be regarded as a projection method that, at the $k$th iteration, computes

[TABLE]

where $\perp_{{\bm{\Theta}}^{(m)}}$ denotes orthogonality with respect to the second sketched inner product in (15). This generalizes LSQR, for which the constraints on the residual are given by ${\bf r}_{k}\perp{\bf A}\mathcal{K}_{k}({\bf A}^{\top}{\bf A},{\bf A}^{\top}{\bf b})$. According to the general theory of projection methods, (24) implies that the ${\bm{\Theta}}^{(m)}$-seminorm of ${\bf r}_{k}$ is minimized at each iteration, which agrees with (23). Eventually, one could also derive rLSQR starting from the normal equations associated to (23), i.e., by equaling the gradient of the objective function in (23) to zero, i.e.,

[TABLE]

The $k$th rLSQR approximate solution is such that

[TABLE]

Note that the standard $l_{2}$ inner product is used to impose the orthogonality condition above. Exploiting the rGK factorizations (17) we recover the projected problem

[TABLE]

which is equivalent to (22). For ${\bm{\Theta}}^{(m)}={\bf I}_{m}$, i.e., in a deterministic setting, conditions (24) and (25) coincide with the ones imposed in a deterministic setting to derive the mathematically equivalent LSQR and CGLS methods, respectively.

We next derive a rGK-based solver starting from the conditions

[TABLE]

where $\perp_{{\bm{\Theta}}^{(n)}}$ denotes orthogonality with respect to the first sketched inner product in (15). Exploiting the rGK factorizations (17), imposing the above conditions amounts to computing

[TABLE]

and where ${\bf V}_{k}$ and ${\bf M}_{k}$ are obtained running $k$ iterations of Algorithm 2, $\widetilde{{\bf T}}_{k+1}=[{\bf I}_{k},{\bf 0}]{\bf T}_{k+1}\in\mathbb{R}^{k\times(k+1)}$ is ${\bf T}_{k+1}$ without its last row, $t_{1,1}$ is the (1,1) entry of ${\bf T}_{k+1}$, and $\beta=\|{\bm{\Theta}}^{(m)}{\bf b}\|$. For ${\bm{\Theta}}^{(n)}={\bf I}_{n}$ and ${\bm{\Theta}}^{(m)}={\bf I}_{m}$, i.e., in a deterministic setting, conditions (26) would be equivalent to (25), which would underlie the derivation of CGLS (mathematically equivalent to LSQR). However, because of randomness, (26) leads to the solver (27) that is not mathematically equivalent to rLSQR. Because of the analogy with the deterministic setting, we regard ${\bf x}_{k}$ in (27) as the $k$th approximate solution computed by the new randomized CGLS (rCGLS) method. Note, however, that rCGLS does not bear other resemblances with standard CGLS, including the lack of short-recurrences for the solution updates. Equivalently, one could derive rCGLS starting from the equations

[TABLE]

and imposing similar conditions to the ones in (26), but with orthogonality of the residual in the standard $l_{2}$ inner product. Straightforward derivations involving the rGK partial factorizations (17) show that one obtains the same projected problem as in (27). The lack of optimality properties of this solver is evident (also) from the equation above, as the coefficient matrix $({\bm{\Theta}}^{(n)})^{\top}{\bm{\Theta}}^{(n)}{\bf A}^{\top}{\bf A}$ is not even symmetric (unless $({\bm{\Theta}}^{(n)})^{\top}({\bm{\Theta}}^{(n)})$ and ${\bf A}^{\top}{\bf A}$ commute, which is very unlikely). Indeed, in this setting, none of the well-known optimality results for general iterative projection methods apply; see [32, Chapter 5].

We finally derive a new randomized LSMR (rLSMR) method that, at the $k$th iteration, defines the $k$th approximate solution of (2) by imposing the conditions

[TABLE]

Running $k\leq K$ iterations of Algorithm 2 and taking ${\bf x}={\bf V}_{k}{\bf z}$, with ${\bf z}\in\mathbb{R}^{k}$, the second condition above and the rGK factorizations (17) lead to the linear system of equations

[TABLE]

which can be further simplified to compute the $k$th rLSMR solution as

[TABLE]

or, equivalently,

[TABLE]

also here, $\beta=\|{\bm{\Theta}}^{(m)}{\bf b}\|$ and $t_{1,1}$ is the (1,1) entry of ${\bf T}_{k+1}$. In exact arithmetic, it can be easily seen that rLSMR minimizes the norm of the sketched normal equations residual over the subspace $\mathcal{K}_{k}({\bf A}^{\top}{\bf A},{\bf A}^{\top}{\bf b})$. Indeed,

[TABLE]

We note that both conditions (28) and the optimality property (30) generalize to the randomized setting the standard LSMR ones, which would hold with ${\bm{\Theta}}^{(n)}={\bf I}_{n}$.

Although the new rGK-based solvers introduced in this section are natural generalizations of their classical counterparts, their implementations may significantly differ. Most notably, as a consequence of the lack of short recurrences for the updates of the basis vectors ${\bf v}_{i}$, ${\bf u}_{i}$, $i=1,\dots,k+1$ (and lack of bidiagonal structure in the matrices ${\bf M}_{k}$ and ${\bf T}_{k+1}$), the entire basis should be stored as the iterations proceed, so that the storage cost increases linearly with the number of iterations. Moreover, the sketched basis vectors collected in the matrices ${\bf P}_{k+1}\in\mathbb{R}^{\ell_{n}\times(k+1)}$ and ${\bf S}_{k+1}\in\mathbb{R}^{\ell_{m}\times(k+1)}$ should be stored (although, if $\ell_{n}\ll n$ and $\ell_{m}\ll m$, their overhead is minimal). These are the trade-offs of having an $l_{2}$ inner product free and randomized method, but in the inverse problems setting, they are mitigated by the following facts. First of all, as mentioned at the end of Section 2.2, in many situations it may be more appropriate to run ro-GKB rather than standard GKB, demanding the storage of all ${\bf v}_{i}$ and ${\bf u}_{i}$ vectors anyway. Second, when performing (standard) GKB-based hybrid projection methods (described for the randomized case in the next section), in order to allow adaptive regularization parameter choice strategies at each iteration, one has to store the solution basis vectors ${\bf v}_{i}$. We also mention that similar storage requirements occur in flexible Golub-Kahan methods to handle Tikhonov-regularized problems more general than (3); see [9]. Such memory demands may be mitigated, also in the randomized setting, by appropriately using some recycling strategies; see also [23].

SVD approximation

A common way to assess the regularization capabilities of iterative projection methods is to study the quality of the approximation of the dominating singular values of ${\bf A}$ obtained by taking the singular values of the corresponding projected matrices; see, for instance, [18, 19, 20]. Although a theoretical investigation of such properties for the new rGK-based solvers is outside the scope of this work, we provide an illustration to empirically show the different performance of such solvers.

We generate a sparse matrix ${\bf A}$ of size $8192\times 4096$ modelling a seismic tomography problem (see Section 4, Example 2, for a description of a larger instance of this test problem). In all the frames in Fig. 2 we plot the dominating singular values of ${\bf A}$ with solid vertical lines. Then, for $k=2,4,6,8,10$, $k$ different markers highlight the singular values of the $k$th projected matrices associated to different solvers (we show only the ones falling within the dominating singular values of ${\bf A}$). As $k$ increases, we observe that the rCGLS and rLSMR projected matrices provide better approximations than rLSQR to the larger singular values of ${\bf A}$. This may affect the behavior of these solvers when employed as iterative regularization methods, and when employed in a hybrid fashion (see the next section), especially when it comes to choosing the regularization parameter for the latter.

3 Hybrid randomized Krylov methods

In this section, we describe how the solvers in Section 2.1 and Section 2.3 can be applied to approximate the solution of the standard Tikhonov regularized problem Eq. 3, with possible adaptive choice of the regularization parameter $\lambda\geq 0$, in the framework of hybrid projection methods, as outlined in [10]. For most of this section, we will naturally consider the first formulation in (3).

For standard Krylov subspace methods, one starting point to introduce hybrid projection methods, which also justifies them from a theoretical point of view, is the fact that first projecting the original problem and then regularizing the projected problem is equivalent to first regularizing the original problem and then projecting the regularized problem; see [20]. In the following we show that, when considering the randomized Krylov subspace methods introduced in the previous subsections, this is only true when considering specific instances of using the sketched seminorms $\|{\bm{\Theta}}^{(m)}\cdot\|$ and $\|{\bm{\Theta}}^{(n)}\cdot\|$. We assume, for now, that $\lambda$ is fixed.

We begin by considering projection methods based on the rArnoldi factorization (10), and we show that regularizing the projected problem is not equivalent to projecting the regularized problem (3). That is, consider projecting the LS problem,

[TABLE]

onto the Krylov subspace $\mathcal{K}_{k}({\bf A},{\bf b})$, $k\leq K$, which amounts to solving

[TABLE]

We denote the project-then-regularize solution ${\bf x}_{k,\lambda}:={\bf V}_{k}{\bf z}_{k,\lambda}\in\mathcal{K}_{k}({\bf A},{\bf b})$, where

[TABLE]

On the other hand, if we start with the Tikhonov problem (3) and project it using rArnoldi (i.e., taking ${\bf x}\in\mathcal{R}({\bf V}_{k})$ in (3)), we get the regularize-then-project solution ${\bf x}_{\lambda,k}:={\bf V}_{k}{\bf z}_{\lambda,k}\in\mathcal{K}_{k}({\bf A},{\bf b})$, where

[TABLE]

Although both provide approximations to the regularized solution ${\bf x}_{\text{\scriptsize reg}}$ in (3), it is clear that ${\bf x}_{\lambda,k}\neq{\bf x}_{k,\lambda}$, as problems (31) and (32) are not equivalent and cannot be further simplified (because ${\bf V}_{k+1}$ does not have orthonormal columns). However, since ${\bf V}_{k}$ has ${\bm{\Theta}}^{(n)}$-orthonormal columns, replacing the regularization term in (3) and (32) with the ${\bm{\Theta}}^{(n)}$-seminorm results in the same projected problem,

[TABLE]

A further simplified projected problem can be obtained by replacing also the fit-to-data norm in (3), (31), and (32) with the ${\bm{\Theta}}^{(n)}$- seminorm. That is, the $k$th hybrid rGMRES iterate is given by ${\bf x}_{k}={\bf V}_{k}{\bf z}_{k}$, where

[TABLE]

and it holds that ${\bf x}_{k}$ is a minimizer of

[TABLE]

which is equivalent to damped LS problem,

[TABLE]

Moving to the rGK-based solvers, we can consider hybrid variants, where all the above rArnoldi derivations (including the lack of equivalence between the project-then-regularize and regularize-then-project approaches) generalize to this case. After running $k\leq K$ rGK iterations, we define the $k$th hybrid rLSQR iterate by ${\bf x}_{k}={\bf V}_{k}{\bf z}_{k}$, where

[TABLE]

and we have that

[TABLE]

or, equivalently,

[TABLE]

Alternatively, in the rGK setting, one could consider applying rLSQR directly to the second Tikhonov problem formulation in (3), i.e., Tikhonov formulated as a damped least-squares problem. At the $k$th rLSQR iteration, we get that

[TABLE]

is an approximation to ${\bf x}_{\text{\scriptsize reg}}$, where ${\bm{\Theta}}^{(m+n)}\in\mathbb{R}^{\ell_{m+n}\times(m+n)}$ is a $(\varepsilon,\delta,K+1)$-oblivious sketching matrix for vectors in $\mathbb{R}^{m+n}$. Note that problem (38) can be regarded as a special case of (39), as $\mathcal{K}_{k}(({\bf A}^{\top}{\bf A}+\lambda^{2}{\bf I}_{n}),{\bf A}^{\top}{\bf b})=\mathcal{K}_{k}({\bf A}^{\top}{\bf A},{\bf A}^{\top}{\bf b})$. However the advantage of the hybrid method defined in (38) over the sketched Tikhonov formulation (39) is that the former naturally allows to set the regularization parameter $\lambda$ adaptively; this is not so straightforward when considering the latter, as the matrices ${\bf M}_{k}$ and ${\bf T}_{k+1}$ appearing in the rGK decomposition (17) depend on $\lambda$; see also [14]. To the best of our knowledge, most of the other randomized methods for Tikhonov regularization consider a sketch-and-solve approach and are naturally formulated as (39) without the constraint on the approximate solution (i.e., taking ${\bf x}\in\mathbb{R}^{n}$), which is not necessarily computed through an iterative method. None of these solvers allow for easy tuning of the regularization parameter; see, e.g., [15, 21].

One could also consider a hybrid version of rCGLS and rLSMR that incorporates Tikhonov regularization for approximating ${\bf x}_{\text{\scriptsize reg}}$. Starting with rCGLS, at the $k$th iteration, $k\leq K$, condition (26) can be applied where we replace ${\bf A}$ with the matrix $[{\bf A}^{\top},\lambda{\bf I}_{n}]^{\top}$ and ${\bf r}_{k}$ with the vector $[({\bf b}-{\bf A}{\bf x}_{k})^{\top},\lambda{\bf x}_{k}^{\top}]^{\top}$, both naturally appearing in the second formulation of the Tikhonov problem (3). Then, using the properties of the matrices in the rGK factorization (17), we get the $k$th rCGLS approximation applied to the leftmost problem in (3) as

[TABLE]

where ${\bf z}_{k}$ is such that

[TABLE]

The above equation can be regarded as a regularization of the projected system appearing in (27), and we refer to the solver so obtained as hybrid rCGLS. Finally, one can regularize the rLSMR projected problem (29) by taking ${\bf x}_{k}={\bf V}_{k}{\bf z}_{k}\in\mathcal{K}_{k}({\bf A}^{\top}{\bf A},{\bf A}^{\top}{\bf b})$, where

[TABLE]

For the hybrid rLSMR method so defined it holds that

[TABLE]

which can be regarded as a Tikhonov-regularized version of problem (30). Contrary to hybrid rLSQR and hybrid rCGLS, hybrid rLSMR is not linked to the original Tikhonov formulation because the fit-to-data term is defined with respect to the normal equations residual. Indeed, even for (standard) LSMR, it is well-known that project-then-regularize is not equivalent to regularize-then-project; see [12] for more details. When rLSMR is concerned, we just briefly mention that one may obtain a different regularize-then-project rLSMR approximation of ${\bf x}_{\text{\scriptsize reg}}$ in (3) by applying rLSMR to the matrices appearing in the leftmost problem in (3) (i.e., following an approach similar to the one described for hybrid rCGLS).

All of the hybrid solvers developed in this section for the Tikhonov-regularized problem compute solutions that belong to approximation subspaces that are independent of the Tikhonov regularization parameter $\lambda$. This property allows us to cheaply approximate different solutions for different values of $\lambda$ and/or to adaptively set $\lambda$ at every iteration by solving the associated projected regularized problems for different values of $\lambda$ without having to recompute the corresponding approximation subspace; as far as the number of iterations $k\leq K\ll\min(m,n)$, the computational cost of doing this is negligible. Many parameter choice rules commonly used to set $\lambda$ for the full-dimensional Tikhonov problem have been successfully applied (with appropriate modifications) to adaptively set $\lambda$ at each iteration of standard (deterministic) hybrid methods; see [10, 17, 19] and the references therein. Here we focus on the discrepancy principle and (weighted) GCV only, and briefly describe how they should be modified for use together with the new hybrid randomized solvers. In order to highlight the dependence of the solution ${\bf x}_{k}$ on $\lambda$ and to establish a common framework for all the considered hybrid randomized formulations, we use the notation,

[TABLE]

where a summary of the new hybrid randomized methods and the corresponding definitions of ${\bf V}_{k}$, ${\bf P}_{k}$ and ${\bf c}_{k}$ are provided in Table 1.

The discrepancy principle can be applied if a good estimate of $\|{\bf e}\|$ (the norm of the noise affecting the data in (1)) is available. For hybrid randomized methods, at the $k$th iteration, $k\leq K$, we find the $\lambda$ that satisfies the nonlinear equation,

[TABLE]

by applying a zero finder. The constant $\tau>1$ appearing in the above equation is a safety factor that accounts for possible inaccuracies in the estimation of $\|{\bf e}\|$ and the fact that the data fit is evaluated in the sketched seminorm.

The weighted generalized cross-validation (WGCV) method can be applied when the noise level is unknown. At the $k$th iteration, $k\leq K$, of the considered hybrid randomized methods, WGCV prescribes to minimize, with respect to $\lambda$, the functional

[TABLE]

The numerator of $G_{w}(\lambda)$ can be expressed as in (45); the matrix ${\bf A}_{\text{\scriptsize reg}}^{\dagger}(\lambda)$ appearing in the denominator of $G_{w}(\lambda)$ is such that ${\bf x}_{k}(\lambda)={\bf A}_{\text{\scriptsize reg}}^{\dagger}(\lambda){\bf b}$ and is specific of each hybrid method. Namely,

[TABLE]

for hybrid rGMRES, rLSQR, rCGLS and rLSMR, respectively. For the denominator of Eq. 46, the trace of these matrices can be computed efficiently using the cyclic property. The weight $w$ in (46) can be adaptively set as described in [11]; choosing $w=1$ in (46) results in the (standard) GCV.

We conclude this section by emphasizing that the functionals in (45) and (46) can be efficiently evaluated for different values of $\lambda$, with negligible computational overhead.

4 Numerical Results

In this section, we illustrate the effectiveness of the randomized Krylov methods and their hybrid variants for solving inverse problems, compared to their non-randomized counterparts. We judge the reconstruction quality by the $l_{2}$-norm relative reconstruction error, i.e., by computing $\|{\bf x}_{\text{\scriptsize rec}}-{\bf x}_{\text{\scriptsize true}}\|/\|{\bf x}_{\text{\scriptsize true}}\|$, where ${\bf x}_{\text{\scriptsize rec}}$ is the reconstruction computed by each solver. We consider two test problems, generated via IR Tools [16]: namely, an image deblurring problem and a seismic tomography problem. Moreover, for Tikhonov regularization, we consider randomized methods with fixed $\lambda$ as well as hybrid randomized methods with automatic selection of the regularization parameter. When naming the methods, all randomized approaches have an ‘r’ in front. All hybrid methods are denoted with a dash and appended parameter or parameter choice method. For example, rGMRES-$\lambda$ is rGMRES with a fixed regularization parameter at each iteration, and rGMRES-opt is rGMRES where the optimal regularization parameter is selected at each iteration (i.e., the one minimizing the reconstruction error at each iteration). Other combinations include rGMRES-gcv and rGMRES-dp, where the GCV and the discrepancy principle are applied to select the regularization parameter at each iteration, respectively. LSQR and rLSQR can replace GMRES and rGMRES in the acronyms listed so far. Our implementations of the methods described in this paper, together with a selection of the tests illustrated in this section, will be available at https://github.com/juliannechung/randKrylov.

4.1 Experiment 1: rGMRES and its hybrid version

Here we consider an example from image deblurring and investigate the performance of rGMRES and hybrid rGMRES. The deblurring example is given in Fig. 3, where the goal is to reconstruct the true image of size 512 $\times$ 512 pixels (displayed in the leftmost frame), which was corrupted by mild blur defined by the point spread function (PSF, displayed in the middle frame), and additive Gaussian noise, so that the observation (displayed in the rightmost frame) contains $1\%$ noise level.

Since the forward matrix ${\bf A}$ is square of order $262144$, we compare the rGMRES and hybrid rGMRES methods to the GMRES and hybrid GMRES methods for this problem. We perform at most $K=50$ iterations of these solvers. The relative reconstruction error norms per iteration are displayed in Figure 4. We observe that both GMRES and rGMRES exhibit the so-called semiconvergence behavior, where the relative error norm decreases in early iterations but then increases rapidly. rGMRES-opt and rGMRES-gcv follow GMRES-opt and GMRES-gcv respectively. Although, in later iterations, the randomized approaches begin to deviate from the nonrandomized counterparts, the error is small. This difference is not noticeable from the approximate relative residual norms (not shown here). For all of these experiments, we use the SHRT with $\ell_{n}=13106.$ This is 5% of the vector dimension $n$.

Next, we investigate the impact of $\ell_{n}$ on the rGMRES results. Specifically, in Fig. 5, we consider different values for $\ell_{n}$ and provide relative reconstruction error norms per iteration for the case where the optimal regularization parameter (left plot) and the GCV computed parameter (right plot) are selected at each iteration. For the former we observe that, if $\ell_{n}$ is too small, since we can expect poor inner product approximations during the iterative method, the error behavior is unstable. For large values of $\ell_{n}$, the algorithm performs similarly to GMRES-opt (blue dashed line in Fig. 4). Nevertheless, there are a range of values for $\ell_{n}$ that give reasonable results. We note that $\ell_{n}=319$ corresponds to the default value (21), also used in [4]. The other values for $\ell_{n}$ are selected to be $0.01\%,0.5\%,1\%,5\%$ and $10\%$ of the dimension $n$. Similar results are observed for rGMRES-GCV (shown in the right plot of Fig. 4) and rGMRES without Tikhonov regularization (not shown due to early semiconvergence behavior).

Although the GMRES and the randomized GMRES approaches only work for square problems, they can be useful for inverse problems where ${\bf A}^{\top}$ is not accessible or where ${\bf A}$ is normal and close to being symmetric; see [22]. In the randomized context, we have shown that rGMRES and hybrid variants can be used to solve large-scale inverse problems, where the main benefit compared to standard methods is that sketched inner products are replaced by approximate inner products. We observe that for very small $\ell_{n}$ (relative to the size of the problem), reconstruction errors in the early iterations of rGMRES follow the error curve for GMRES, and that for hybrid rGMRES approaches, the regularization parameter can be selected using a standard method (e.g., GCV). Further numerical investigations for rGK based solvers are provided in the next section.

4.2 Experiment 2: rGK-based solvers and their hybrid versions

The second experiment consists of a seismic travel-time tomography example, where the aim is to reconstruct an image of size 256 $\times$ 256 pixels representing slowness (the reciprocal of the sound speed), from observed data collected along waves modeled as straight lines. Noise was added to the data, in such a way that the observation contained a $4\%$ noise level. The true image and observations are provided in Fig. 6. The corresponding matrix ${\bf A}$ is of size $131,072\times 65,536$.

First we investigate the performance of rLSQR with different sample sizes. The maximum number of iterations is set to $100$, so the sampling parameters automatically selected according to (21) are $\ell_{m}=512$ and $\ell_{n}=482$. We consider sample sizes of $\ell_{m}/\ell_{n}=40/20$, $92/46$, and $5243/2622$, corresponding to $0.03\%$, $0.07\%$ and $4\%$ of $m/n$, respectively. The relative reconstruction error norms per iteration are provided in the left frame of Figure 7. Notice the expected semiconvergence behavior in all plots, where early iterations perform similarly, even for very small sample sizes. However, for very small sample sizes, errors increase quickly and exhibit erratic behavior. Somewhat similarly, for the hybrid rLSQR method with the optimal regularization parameter selected at each iteration, in the right frame of Figure 7 we see that early reconstructions are good even for small sample sizes, but errors start to increase quickly after, and eventually, minimal improvements are made. On the other hand, for a sufficiently large sample size, reconstruction errors remain small for all iterations: this behavior mimics the hybrid LSQR method with optimal regularization parameter (LSQR-opt), which is provided for comparison. For a visual comparison, in Figure 8 we provide the reconstruction from rLSQR-opt at 100 iterations for $\ell_{n}=482$ and $\ell_{n}=2622$, along with the LSQR-opt reconstruction. We observe that, as anticipated, some details are lost in the rLSQR-opt reconstructions, especially for small sample parameters, but both reconstructions are still reasonable.

Next, we investigate the effect of the choice of the regularization parameter $\lambda$ on the hybrid rGK-based solvers. For these experiments, we allow the automatically selected sample size $\ell_{m}/\ell_{n}=512/482$.

We first consider the case where a fixed regularization parameter is used at each iteration. In the left frame of Figure 9 we provide relative reconstruction errors for rLSQR-$\lambda$ for different values of $\lambda$. We note that $22.3951$ is the optimal regularization parameter for rLSQR-opt at iteration $100$, and $30.6996$ is the optimal regularization parameter for LSQR-opt at iteration $100$. For regularization parameters that are too small, we observe semiconvergence, and for regularization parameters that are too large we obtain overly regularized solutions (again corresponding to large errors). Still in Figure 9, in the right frame, for comparison, we investigate how rLSQR performs when applied to the Tikhonov problem, for a fixed regularization parameter and as described in Eq. 39. We use the same set of regularization parameters as in the left plot of Figure 9. Note that, since the augmented matrix in Eq. 39 has $m+n$ rows, the sample sizes selected according to (21) are $\ell_{m+n}/\ell_{n}=530/482$. We observe that applying rLSQR to the Tikhonov problem results in more irregular behavior and overall larger errors than rLSQR-$\lambda$.

To back up the results displayed in the left frame of Fig. 9, since there is randomness involved in the sketching process, we run the same solver $24$ times for two fixed values of the regularization parameter (among the ones considered in Figure 9) and, in Figure 10, we provide box plots for the relative reconstruction errors at each iteration, highlighting median and the $25$th and $75$th percentiles. Outliers are plotted separately. We see that randomness in the sketching can affect the results, but the general trends remain the same.

Next, we consider different adaptive regularization parameter selection methods. In the left frame of Fig. 11 we show the relative reconstruction errors per iteration of the hybrid rLSQR method. Namely, we provide results for GCV, DP, WGCV, and when the optimal regularization parameter is selected at each iteration. The corresponding regularization parameters computed at each iteration are provided in the right frame of Fig. 11. Note that the optimal regularization parameter is zero in early iterations, since the problem is still well-conditioned and no regularization is needed. The discrepancy principle and the WGCV method provide good approximations of the regularization parameter at later iterations (i.e., close to its optimal values), while the GCV method overestimates the regularization parameter.

Finally, we investigate the performance of the different solvers based on the rGK algorithm. In the left plot of Fig. 12 we provide relative reconstruction errors per iteration for rLSQR, rCGLS, rLSMR and their hybrid versions with the optimal regularization parameter selected at each iteration. In the right frame of Fig. 12, we use the regularization parameter selected by the discrepancy principle at each iteration. We observe that all of the considered methods exhibit similar behavior, with the rLSMR methods having slightly smaller reconstruction errors.

With the seismic tomography example, we have shown that rGK-based solvers can be used to solve inverse problems with rectangular discretized forward operator ${\bf A}$, where regularized solutions can be obtained with early termination (e.g., before the errors increase) or via a hybrid approach, where various parameter selection methods can be included. Compared to standard Krylov methods, all of the randomized Krylov approaches (including the ones tested for the image deblurring example) require the selection of a sampling parameter. For problems where only a few iterations are needed to obtain a reasonable solution (e.g., problems with very high noise levels or severe ill-conditioning), smaller values of the sampling parameter may be sufficient. Choosing a good regularization parameter is also important, although it may not be able to compensate for errors introduced in the basis generation.

5 Conclusions

In this work, we develop randomized Krylov methods for large-scale inverse problems, where sketched inner products are used within a randomized Gram-Schmidt process to generate suitable projected problems. For square problems, we consider randomized GMRES and introduce its hybrid variant for solving Tikhonov-regularized problems. For rectangular matrices, we develop a randomized Golub-Kahan approach, where each iteration requires a matrix-vector multiplication with ${\bf A}$ and ${\bf A}^{\top}$, and all the inner products are approximated with sketched inner products. We describe new iterative solvers based on the the randomized Golub-Kahan method, and describe hybrid variants of all approaches, with automatic Tikhonov regularization parameter selection at each iteration. Numerical results show that randomized Krylov methods can be used to solve inverse problems, taking into account that semiconvergence appears when using randomized GMRES and LSQR, and when the sample size for the randomized inner products approximation is too small. In the latter case, even a good regularization parameter may not be able to provide a good reconstruction if not terminated early.

There are various promising directions for future work. For example, there are many choices for selecting the sketching matrix, and exploiting any matrix structure may improve the results and significantly improve run time. We are currently working on a highly efficient computational implementation and an analysis of the computational complexity of the algorithm. Moreover, there are potential extensions to settings such as Tikhonov regularization in general form and flexible Krylov methods to handle $p$-norm regularization, $0<p\leq 1$. Finally, further research is needed to understand the theoretical properties of the low-rank matrix approximations associated with the new rGK algorithm, which may eventually be useful for estimating matrix functions or for uncertainty quantification.

6 Acknowledgments

We would like to thank Leonardo Robol for many helpful and stimulating discussions about randomized numerical linear algebra.