A novel modified TRSVD method for large-scale linear discrete ill-posed problems

Abstract

The truncated singular value decomposition (TSVD) is a popular method for solving linear discrete ill-posed problems with a small to moderately sized matrix A. This method replaces the matrix A by the closest matrix $A_k$ of low rank k, and then computes the minimal norm solution of the linear system of equations with a rank-deficient matrix so obtained. The modified TSVD (MTSVD) method improves the TSVD method, by replacing A by a matrix that is closer to A than $A_k$ in a unitarily invariant matrix norm and has the same spectral condition number as $A_k$. Approximations of the SVD of a large matrix A can be computed quite efficiently by using a randomized SVD (RSVD) method. This paper presents a novel modified truncated randomized singular value decomposition (MTRSVD) method for computing approximate solutions to large-scale linear discrete ill-posed problems. The rank, k, is determined with the aid of the discrepancy principle, but other techniques for selecting a suitable rank also can be used. Numerical examples illustrate the effectiveness of the proposed method and compare it to the truncated RSVD method.

Introduction

This paper is concerned with the computation of approximate solutions of large minimization problems of the form

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

where $A\in {\mathbb{R}}^{m\times n}$ is a large matrix whose singular values gradually decay to zero without a significant gap. In particular, A is severely ill-conditioned and may be rank-deficient. Minimization problems (1.1) with a matrix of this kind often are referred to as discrete ill-posed problems. They arise, for example, from the discretization of linear ill-posed problems, such as Fredholm integral equations of the first kind with a smooth kernel. Throughout this paper ${\Vert \cdot \Vert }_2$ denotes the Euclidean vector norm or the spectral matrix norm. We allow $m\ge n$ as well as $m<n$.

The vector $b\in {\mathbb{R}}^m$ in (1.1) represents measured data that are contaminated by an error $e\in {\mathbb{R}}^m$, which may stem from measurement or discretization errors. Let $\hat{b}\in {\mathbb{R}}^m$ denote the unknown error-free vector associated with b, i.e.,

$$b=\hat{b}+e.$$

We will assume $\hat{b}$ to be in the range of A, and that a fairly accurate bound for the relative error

$$\frac{{\Vert e\Vert }_2}{{\Vert \hat{b}\Vert }_2}\le ϵ,$$

in b is known. Then a regularized solution of (1.1) can be determined with the aid of the discrepancy principle; see below.

We are interested in computing an approximation of the solution $\hat{x}$ of minimal Euclidean norm of the unknown error-free least-squares problem

$$\underset{x\in {\mathbb{R}}^n}{\min }{\Vert Ax-\hat{b}\Vert }_2.$$

Let $A^{†}\in {\mathbb{R}}^{n\times m}$ denote the Moore-Penrose pseudoinverse of A. Then

$$\hat{x}=A^{†}\hat{b}.$$

Because A has many positive singular values close to zero, the matrix $A^{†}$ is of very large norm, and the solution of the available least-squares problem (1.1), given by

$$\check{x}=A^{†}b=A^{†}(\hat{b}+e)=\hat{x}+A^{†}e,$$

typically is dominated by the propagated error $A^{†}e$, and then is meaningless. This difficulty can be mitigated by replacing the matrix A by a nearby matrix that does not have tiny positive singular values. This replacement commonly is referred to as regularization. One of the most popular regularization methods for discrete ill-posed problem (1.1) of small to moderate size is the truncated singular value decomposition (TSVD); see, e.g., [2], [3], [7]. Introduce the singular value decomposition

$$A=U\mathrm{\Sigma }{V}^{\top },$$

where $U=[u_1,u_2,\dots ,u_m]\in {\mathbb{R}}^{m\times m}$ and $V=[v_1,v_2,\dots ,v_n]\in {\mathbb{R}}^{n\times n}$ are orthogonal matrices, the superscript ^⊤ denotes transposition, and the nontrivial entries (known as the singular values) of the (possibly rectangular) diagonal matrix

$$\mathrm{\Sigma }=\mathrm{diag}[{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_{\min \{m,n\}}]\in {\mathbb{R}}^{m\times n}$$

are ordered according to

$${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_r>{\sigma }_{r+1}=\dots ={\sigma }_{\min \{m,n\}}=0.$$

Here r is the rank of A. Define the matrix

$${\mathrm{\Sigma }}_k=\mathrm{diag}[{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_k,0,\dots ,0]\in {\mathbb{R}}^{m\times n}$$

for $k\le r$ by setting the singular values ${\sigma }_{k+1},{\sigma }_{k+2},\dots ,{\sigma }_{\min \{m,n\}}$ in (1.6) to zero. Then the matrix

$$A_k=U{\mathrm{\Sigma }}_k{V}^{\top }$$

is a best rank-k approximation of A in any unitarily invariant matrix norm, such as the spectral and Frobenius norms. We have

$${\Vert A_k-A\Vert }_2={\sigma }_{k+1},{\Vert A_k-A\Vert }_F=\sqrt{\sum _{i=k+1}^r{{\sigma }_i}^2},0\le k\le r,$$

where ${\Vert \cdot \Vert }_F$ denotes the Frobenius norm, and we define $A_0=0$ and ${\sigma }_{n+1}=0$.

The TSVD method replaces the matrix A in (1.1) by $A_k$. Let $x_k\in {\mathbb{R}}^n$ denote the solution of minimal Euclidean norm of

$$\underset{x\in {\mathbb{R}}^n}{\min }{\Vert A_{k}x-b\Vert }_2.$$

It is given by $x_k=A_k^{†}b$. The discrepancy principle prescribes that the truncation index $k\ge 0$ in (1.7) be chosen as the smallest integer such that

$$\frac{{\Vert A_k{x}_k-b\Vert }_2}{{\Vert \hat{b}\Vert }_2}\le \tau ϵ,$$

where $\tau >1$ is a user-chosen constant that is independent of the bound ϵ in (1.3); see [3], [7]. We will use the discrepancy principle in the computed examples reported in Section 4. However, other methods for determining k also can be used when no accurate estimate of $\Vert e\Vert$ is available, such as the L-curve criterion and generalized cross validation; see, e.g., [2], [7], [11], [12], [14].

The modified truncated singular value decomposition (MTSVD) method proposed in [13] replaces the matrix $A_k$ in (1.8) by a matrix that (generally) is closer to A in a unitarily invariant matrix norm with the same condition number as $A_k$. This replacement often results in more accurate approximations of $\hat{x}$ when the discrepancy principle is used to determine the truncation index k; see [13] for computed examples.

The TSVD and MTSVD methods are not suitable for application to the solution of large-scale discrete ill-posed problems (1.1) due to the high cost of evaluating the SVD of a large matrix A; see, e.g., [5, p. 493] for counts of arithmetic floating point operations. However, the computational effort can be reduced by computing an approximation of the singular value decomposition by a randomized method and applying the modified truncated singular value decomposition to the computed approximate SVD. This yields the modified truncated randomized singular value decomposition (MTRSVD).

In recent years several randomized algorithms have been proposed for computing approximate factorizations of a large matrix, such as an approximate SVD; see, e.g., Halko et al. [6]. These factorizations have been used to compute approximate solutions to ill-posed problems (1.1) by Tikhonov regularization; see, e.g., Jia and Yang [10] and Xiang and Zou [15], [16].

It is the purpose of the present paper to discuss the use of a randomized singular value decomposition (RSVD) in conjunction with the MTSVD described in [13]. We refer to this scheme as the MTRSVD method. Our reason for regularizing by a truncated singular value decomposition instead of using Tikhonov regularization is that the former method is easier to implement. Moreover, we base our method on the MTSVD instead of on the standard TSVD, because the former typically yields approximations of the desired solution (1.4) of higher quality when the required singular vectors have been computed with high enough accuracy.

The remainder of this paper is organized as follows. Section 2 provides a brief review of the MTSVD, RSVD, and TRSVD methods. The proposed MTRSVD method is described in Section 3, where also its computational complexity and error bounds are presented. Section 4 shows several numerical examples that illustrate the efficiency of the MTRSVD method. Section 5 contains concluding remarks.