Yin Wu, Lei Zhao^*
School of Mathematical Sciences, Gansu Minzu Normal University, Hezuo, China.
DOI: 10.4236/ojapps.2026.162034   PDF    HTML   XML   94 Downloads   276 Views   Citations

Abstract

Truncated singular value decomposition (TSVD) and Golub-Kahan diagonalization are two elementary techniques for solving a least squares problem from a linear discrete ill-posed problems. For small to medium sized ill-posed problems, we propose a novel Modified Truncated Singular Value Decomposition (NMTSVD) based on the MTSVD method. This method presents three approaches to select truncation indices, leading to three approximate matrices of the coefficient matrix $A$ . And the relationships of these approximate matrices are established in terms of the spectral condition number. For large-scale unstable discrete ill-posed problems, we propose a truncated Reorthogonalized Golub-Kahan Bidiagonalization (RGKB) method by integrating the Golub-Kahan bidiagonalization process, truncation index, and the Gram-Schmidt reorthogonalization procedure. Numerical experiments are presented to illustrate the effectiveness of the methods NMTSVD and RGKB.

Keywords

Discrete Ill-Posed Problems, Truncated Singular Value Decomposition, Reorthogonalized Golub-Kahan Bidiagonalization, Gram-Schmidt

Share and Cite:

1. Introduction and Main Results

Consider the computation of an approximate solution of the minimization problem:

$\underset{x\in {ℝ}^n}{\text{min}}\Vert Ax-b\Vert ,A\in {ℝ}^{m\times n},b\in {ℝ}^m.$ (1)

where and throughout this paper, $\Vert \cdot \Vert$ denotes the Euclidean vector norm or the associated induced matrix norm. The singular values of the matrix $A$ are assumed of different size closed to the origin. It follows that $A$ is severely ill-conditioned and may be singular. Minimization problems (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 problem, such as Fredholm integral equations of the first kind with a smooth kernel [1]. The application background of discrete ill-posed problems is very extensive, such as signal processing, image denoising and deblurring [2]. We will for notational simplicity assume that $m\ge n$ ; however, the methods discussed also can be applied when $m<n$ .

There are measurement or discretization errors $e$ in $b$ , here $e\in {ℝ}^m$ is called noise, and it is generally Gaussian white noise. Let $b_{true}\in {ℝ}^m$ denotes the (unknown) noise-free vector associated with $b$ , i.e., $b=b_{true}+e$ . We are interested in computing an approximation of the solution $\widehat{x}$ of the minimal Euclidean norm of the error-free least-squares problem:

$\underset{x\in {ℝ}^n}{\text{min}}\Vert Ax-b_{true}\Vert .$ (2)

Let $A^{†}$ denote the Moore-Penrose pseudoinverse of $A$ . Due to the ill-conditioning of $A$ , the solution of (1)

$\widehat{x}=A^{†}b=A^{†}\left(b_{true}+e\right)=x_{true}+A^{†}e,$

often does not yield a meaningful approximation of $x_{true}$ , because $A^{†}e$ causes a relatively large error. We would like to determine an approximate of $\widehat{x}$ by computing a suitable approximate solution of (1) [3].

For small to medium sized discrete ill-posed problems, the Tikhonov regularization method is most commonly used for numerical solution [4]. This method replaces (1) by a penalized least-squares problem

$\underset{x\in {ℝ}^n}{\text{min}}\left\{{\Vert Ax-b\Vert }^2+{\mu }^{-2}{\Vert Lx\Vert }^2\right\},$ (3)

where $L\in {ℝ}^{p\times n}$ ($p\le n$ ) is the regularization matrix, and the scalar $\mu >0$ is the regularization parameter. The normal equations associated with (3) are given by $\left(A^{\text{T}}A+{\mu }^{-2}{L}^{\text{T}}L\right)x=A^{\text{T}}b$ . The matrix $L$ is assumed to satisfy

$\mathcal{N}\left(A\right)\cap \mathcal{N}\left(L\right)=\left\{0\right\},$ (4)

where $\mathcal{N}\left(M\right)$ denotes the null space of matrix $M$ , then (3) has a unique solution

$x_{\mu }={\left(A^{\text{T}}A+{\mu }^{-2}{L}^{\text{T}}L\right)}^{-1}{A}^{\text{T}}b.$

Throughout this paper $A^{\text{T}}$ denotes the transpose of matrix $A$ . The value of $\mu$ derermines how sensitive $x_{\mu }$ is to the error $e$ and how close $x_{\mu }$ is to $\widehat{x}$ , and Hansen for discussions on Tiknonov regularization [5]. Both the choice of the regularization parameter $\mu$ and the regularization matrix $L$ are crucial, because they affect the approximation degree of the solution to (3) with the true solution $x_{true}$ and the sensitivity of the solution of (3) to noise. $L$ is generally selected as an orthogonal matrix or a banded matrix with known null space, such as the identity matrix, first-order or second-order derivative operator matrices. When $L$ is the identity matrix, the Tikhonov minimization problem (3) is said to be in standard form, otherwise it is said to be in general form. When the magnitude of the noise is known, the discrepancy principle is used to determine $\mu$ ; when the magnitude of the noise is unknown, the L-curve criterion and Generalized Cross Validation (GCV) are used to determine $\mu$ , see, e.g., [6]-[9]. There are also other methods, such as the vector extrapolation method and the three-point interpolation zero-finding method [10]-[12].

Regularization methods and Singular Value Decomposition methods are often used for numerical solution of the system (1). Hansen et al. proposed the Truncated Singular Value Decomposition (TSVD) in [5] [13]. This method based on the singular value decomposition of matrix $A$ , find an optimal rank-$k$ matrix $A_k$ to approximate $A$ , so that the approximate problem $A_{k}x=b$ approximates (1). The truncated index $k$ is determined by the discrete Picard condition [14]. Morigi et al. proposed the Truncated Projection SVD method (TPSVD) in [15]. This method applied orthogonal projection to matrix $A$ and the right-hand side vector $b$ to obtain a projection problem, and then apply TSVD to solve it. Reichel proposed the Modified Truncated Singular Value Decomposition (MTSVD) in [16]. This method modify the truncation index and some singular values, so that the coefficient matrices determined by TSVD and MTSVD have the same spectral condition number, but MTSVD can obtain a better approximate solution. Zhongxiao Jia [17] proposed the Modified Truncated Randomized Singular Value Decomposition (MTRSVD) Algorithms for large scale discrete ill-posed problems with general-form regularization.

In this paper, for small-to-medium sized discrete ill-posed problems, we propose a new Modified Truncated Singular Value Decomposition (NMTSVD) base on the MTSVD method. This method presents three approaches to select truncation index and three new partial singular value correction methods, leading to three approximate matrices of the coefficient matrix $A$ . And the relationships of these approximate matrices are established in terms of the spectral condition number. Compared with TSVD and MTSVD, NMTSVD can obtain better approximate solutions.

For the large-scale nostable iterative Tikhonov regularization problem [18].

$\underset{x\in {ℝ}^n}{\text{min}}\left\{{\Vert Ax-b\Vert }^2+{\mu }_k^{-2}{\Vert L\left(x-x_{k-1}\right)\Vert }^2\right\},k=1,2,\cdot \cdot \cdot .$ (5)

where $x_0\in {ℝ}^n$ is the initial approximate solution, generally taken as $x_0=0$ . Only when $A$ and $L$ satisfy (4) can be ensured the system (5) to have the unique solution. The regularization parameter ${\mu }^{-2}$ has advantages over $\mu$ . The selection of the regularization matrix and parameter is crucial, leading to affect the approximation accuracy between the approximate solution of (5) and the exact solution. The unstable iterative Tikhonov regularization method yields solutions with higher approximation accuracy [7] [19]-[21]. The normal equations associated with (5) are given by

$\left(A^{\text{T}}A+{\mu }_k^{-2}{L}^{\text{T}}L\right)x_k={\mu }_k^{-2}{L}^{\text{T}}L{x}_{k-1}+A^{\text{T}}b,k=1,2,\cdot \cdot \cdot .$ (6)

The computed approximation solution $x_{\mu }^k$ , lives in the series of low-rank $k$ -dimensional Krylov subspace

${\mathcal{K}}_k\left(A^{\text{T}}A,A^{\text{T}}b\right)=span\left\{A^{\text{T}}b,\left(A^{\text{T}}A\right)A^{\text{T}}b,\cdots ,{\left(A^{\text{T}}A\right)}^{k-1}{A}^{\text{T}}b\right\}.$ (7)

The Golub-Kahan bidiagonalization (GKB) is popularly applied to reduce large-scale unstable discrete ill-posed problems to small dimension minimization problems. However, during the solution process of the iterative method, a semi-convergence phenomenon may occur in the approximate solution, i.e., the convergence effect of the approximate solution initially improves, but when the number of iterations exceeds a certain value, the approximate solution tends to diverge [22] [23].

The columns of orthogonal matrices $U_k,V_k$ of the $k$ steps of Golub-Kahan bidiagonalization to $A$ tend to lose orthogonality, therefore necessary to apply the reorthogonalization strategy of the columns of $V_k$ to preserve the orthogonality of the orthogonal matrices and convergence of the singular value. Many researchers concern about the reorthogonalization strategy [24] [25] for deeply studying the Golub-Kahan bidiagonal matrix factorization of $A$ . Paige [26] pointed out that the lose of orthogonality in lanczos reductions is structured in the sense that it is coincident with the convergence of approximate eigenvalues and eigenvectors (calls Ritz and vectors). Parlett and Scott [27] used this observation to develop partial reorthogonalization procedures, and gave a good summary of the surrounding issues. Barlow supported three reorthogonalization strategies: complete reorthogonalization, selective reorthogonalization and parameterrized reorthogonalization, see [28] for details. Ramlau [29] proposed the error estimates for Golub-Kahan-digiagonalization for linear ill-posed problems. Reichel [30] proposed the iterated Golub-Kahan-Tikhonov method for large discrete linear ill-posed problems.

This paper combines the GKB method with the reorthogonalization method and truncated the dimension of Krylov subspace, referred to as the Truncated Reorthogonalized Golub-Kahan Bidiagonalization method (RGKB). Compared with GKB, RGKB yields a better approximate solution.

This paper is organized as follows. In Section 2, we propose a novel Modified Truncated Singular Value Decomposition (NMTSVD) methods based on MTSVD for small-to-medium sized discrete ill-posed problems, presenting the filter factor representations for their approximate solutions. In Section 3, introduce the Truncated Reorthogonalized Golub-Kahan Bidiagonalization (RGKB) method for large-scale discrete ill-posed problems, and the RGKB algorithm is presented. In Section 4, present numerical experiments demonstrating the feasibility and effectiveness of NMTSVD and RGKB.

2. Novel Modified Truncated Singular Value Decomposition (NMTSVD)

2.1. NMTSVD Method

First, let’s review the truncated singular value decomposition (TSVD) method and Modified TSVD (MTSVD) method.

TSVD method: Given the singular value decomposition (SVD) of matrix $A\in {ℝ}^{m\times n}$ .

$A=U\left(\begin{array}{c}\Sigma \\ 0\end{array}\right)V^{\text{T}},$ (8)

where $U=\left[u_1,u_2,\cdot \cdot \cdot ,u_m\right]\in {ℝ}^{m\times m}$ and $V=\left[v_1,v_2,\cdot \cdot \cdot ,v_n\right]\in {ℝ}^{n\times n}$ are both orthogonal matrices, i.e., $U^{\text{T}}U=I$ and $V^{\text{T}}V=I$ , where $I$ is the identity matrix. The diagonal elements ${\sigma }_i$ of the diagonal matrix $\Sigma =diag\left[{\sigma }_1,{\sigma }_2,\cdot \cdot \cdot ,{\sigma }_n\right]\in {ℝ}^{n\times n}$ are all the singular values of the coefficient matrix $A$ . The rank of matrix $A$ is $\ell$ , and the singular values of $A$ are ordered as follows:

${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_{\ell }>{\sigma }_{\ell +1}=\cdots ={\sigma }_n=0.$

Define the matrix,

${\Sigma }_k=diag\left[{\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_k,0,\cdots ,0\right]\in {ℝ}^{n\times n}$

by setting the singular values ${\sigma }_{k+1},{\sigma }_{k+2},\cdots ,{\sigma }_n$ to zero. The matrix $A=U\left(\begin{array}{c}{\Sigma }_k\\ 0\end{array}\right)V^{\text{T}}$ is the best rank-$k$ approximation of $A$ in any unitarily invariant matrix norm ,the spectral and Frobenius norm, i.e.,

${\Vert A-A_k\Vert }_2={\sigma }_{k+1},{\Vert A-A_k\Vert }_F=\sqrt{{\displaystyle \sum }_{i=k+1}^n{\left({\sigma }_i\right)}^2}.$

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

The TSVD method replaces the matrix $A$ in (1) by $A_k$ and determines the least-square solution $x_{TSVD}$ of minimal Euclidean norm.

$x_{TSVD}={\displaystyle \sum }_{i=1}^{\ell }{\varphi }_T\frac{u_i^{\text{T}}b}{{\sigma }_i}{v}_i,$

where $u_i$ and $v_i$ are columns of the matrices $U$ and $V$ in (8), respectively, and ${\varphi }_T$ are filter factors defined by

${\varphi }_T=\{\begin{array}{ll}1,\hfill & 1\le j\le k,\hfill \\ 0,\hfill & k<j\le \ell .\hfill \end{array}$

MTSVD method: A closest matrix to $A$ in the spectral or Frobenius norms with smallest singular value ${\sigma }_k$ is given by

$A_{\overline{k}}=U{\overline{\Sigma }}_{\overline{k}}V,$

where $U$ and $V$ are the orthogonal matrix in the SVD (8) of $A$ and ${\overline{\Sigma }}_{\overline{k}}$ has the entries

${\overline{\sigma }}_j={\sigma }_j,1\le j\le k,$

${\overline{\sigma }}_j={\sigma }_k,k\le j\le \overline{k},$

${\overline{\sigma }}_j={\sigma }_j,\overline{k}\le j\le n,$

where $\overline{k}$ is determined by the inequalities ${\sigma }_{\overline{k}}\ge \frac{{\sigma }_k}{2}$ and ${\sigma }_{\overline{k}+1}\le \frac{{\sigma }_k}{2}$ .

We have in the spectral norm and in the Frobenius norm

${\Vert A-A_{\overline{k}}\Vert }_2\le \frac{{\sigma }_k}{2},{\Vert A-A_k\Vert }_2={\sigma }_{k+1},$

${\Vert A-A_{\overline{k}}\Vert }_F\le \frac{{\sigma }_k}{2}\sqrt{n-k},{\Vert A-A_k\Vert }_F=\sqrt{{\displaystyle \sum }_{i=k+1}^n{\left({\sigma }_i\right)}^2}.$

Moreover,

${\Vert A-A_{\overline{k}}\Vert }_2\le {\Vert A-A_k\Vert }_2,{\Vert A-A_{\overline{k}}\Vert }_F\le {\Vert A-A_k\Vert }_F,$

${\kappa }_2\left(A_{\overline{k}}\right)={\kappa }_2\left(A_k\right).$

The approximate solution can be expressed as,

$x_{MTSVD}={\displaystyle \sum }_{i=1}^{\ell }{\varphi }_M\frac{u_i^{\text{T}}b}{{\sigma }_i}{v}_i,$

with the filter factors

${\varphi }_M=\{\begin{array}{ll}1,\hfill & 1\le j\le k,\hfill \\ \frac{{\sigma }_j}{{\sigma }_k},\hfill & k\le j\le \overline{k},\hfill \\ 0,\hfill & \overline{k}<j\le \ell .\hfill \end{array}$

Then,

${\varphi }_T\le {\varphi }_M,1\le j\le \ell ,$

and,

$\frac{1}{2}\le {\varphi }_M\le 1,k<j\le \overline{k}.$

Finally, this paper, we propose a novel modified truncated singular value decomposition (NMTSVD). First, a new selection method for the truncation index is given. The truncation index $k_1$ determined from MTSVD method:

$\begin{array}{l}{k}_1:{\sigma }_{k_1}\ge \frac{{\sigma }_k}{2}\text{and}{\sigma }_{k_1+1}\le \frac{{\sigma }_k}{2}.\\ k_2:{\sigma }_{k_2}\ge \frac{\sqrt{\left({\sigma }_k^2+{\mu }^2\right)/2}}{2}\text{and}{\sigma }_{k_2+1}\le \frac{\sqrt{\left({\sigma }_k^2+{\mu }^2\right)/2}}{2}.\\ k_3:{\sigma }_{k_3}\ge \frac{{\sigma }_k+\mu }{4}\text{and}{\sigma }_{k_3+1}\le \frac{{\sigma }_k+\mu }{4}.\\ k_4:{\sigma }_{k_4}\ge \frac{\mu }{2}\text{and}{\sigma }_{k_4+1}\le \frac{\mu }{2}.\end{array}$ (9)

where $\mu$ is the regularization parameter determined via the discrepancy principle, and the truncation index $k$ is determined by the following MATLAB command,

$k=\max \left(find\left(s>\mu \right)\right),$

where $s=\left[{\sigma }_1,{\sigma }_2,\cdot \cdot \cdot ,{\sigma }_n\right]$ . That is, the truncation index $k$ satisfies:

${\sigma }_k\ge \mu \ge {\sigma }_{k+1}.$ (10)

Next, the modification methods for partial singular values are presented. Four approaches to selecting partial singular values are outlined below. The first is the diagonal matrix ${\Sigma }_1$ from MTSVD after modifying partial singular values. The diagonal matrices from the NMTSVD methods after modifying partial singular values are denoted as ${\Sigma }_2,{\Sigma }_3,{\Sigma }_4$ . The specific modification methods are as follows:

${\Sigma }_1=\{\begin{array}{ll}{\sigma }_j,\hfill & 1\le j\le k,\hfill \\ {\sigma }_k,\hfill & k<j\le k_1,\hfill \\ 0,\hfill & k_1<j\le n.\hfill \end{array}$

${\Sigma }_2=\{\begin{array}{ll}{\sigma }_j,\hfill & 1\le j\le k,\hfill \\ \sqrt{\left({\sigma }_k^2+{\mu }^2\right)/2},\hfill & k<j\le k_2,\hfill \\ 0,\hfill & k_2<j\le n.\hfill \end{array}$

${\Sigma }_3=\{\begin{array}{ll}{\sigma }_j,\hfill & 1\le j\le k,\hfill \\ \left({\sigma }_k+\mu \right)/2,\hfill & k<j\le k_3,\hfill \\ 0,\hfill & k_3<j\le n.\hfill \end{array}$

${\Sigma }_4=\{\begin{array}{ll}{\sigma }_j,\hfill & 1\le j\le k,\hfill \\ \mu ,\hfill & k<j\le k_4,\hfill \\ 0,\hfill & k_4<j\le n.\hfill \end{array}$

From (10), the modified partial singular values ${\sigma }_k,\sqrt{\frac{{\sigma }_k^2+{\mu }^2}{2}},\frac{{\sigma }_k+\mu }{2},\mu$ satisfy the relationship,

${\sigma }_k\ge \sqrt{\frac{{\sigma }_k^2+{\mu }^2}{2}}\ge \frac{{\sigma }_k+\mu }{2}\ge \mu .$ (11)

The corresponding relationships among the four truncation indices are red

$k\le k_1\le k_2\le k_3\le k_4.$ (12)

The four approximate matrices of $A$ are,

$A_1=U\left(\begin{array}{c}{\Sigma }_1\\ 0\end{array}\right)V^{\text{T}},A_2=U\left(\begin{array}{c}{\Sigma }_2\\ 0\end{array}\right)V^{\text{T}},$

$A_3=U\left(\begin{array}{c}{\Sigma }_3\\ 0\end{array}\right)V^{\text{T}},A_4=U\left(\begin{array}{c}{\Sigma }_4\\ 0\end{array}\right)V^{\text{T}}.$

where $U,$ $V$ are the orthogonal matrices in the singular value decomposition of $A$ .

Under the spectral norm, the relationships between $A$ and the approximate matrices $A_i,i=1,2,3,4$ are respectively:

${\Vert A-A_1\Vert }_2={\Vert \Sigma -{\Sigma }_1\Vert }_2={\sigma }_{k_1+1}\le \frac{{\sigma }_k}{2},$

${\Vert A-A_2\Vert }_2={\Vert \Sigma -{\Sigma }_2\Vert }_2=\max \left\{\sqrt{\frac{{\sigma }_k^2+{\mu }^2}{2}}-{\sigma }_{k_2},{\sigma }_{k_2+1}\right\}\le \frac{\sqrt{\left({\sigma }_k^2+{\mu }^2\right)/2}}{2},$

${\Vert A-A_3\Vert }_2={\Vert \Sigma -{\Sigma }_3\Vert }_2=\max \left\{\frac{{\sigma }_k+\mu }{2}-{\sigma }_{k_3},{\sigma }_{k_3+1}\right\}\le \frac{{\sigma }_k+\mu }{4},$

${\Vert A-A_4\Vert }_2={\Vert \Sigma -{\Sigma }_4\Vert }_2=\max \left\{\mu -{\sigma }_{k_4},{\sigma }_{k_4+1}\right\}\le \frac{\mu }{2}.$

Then, under the Frobenius norm, the relationships between $A$ and the approximate matrices $A_i,i=1,2,3,4$ are respectively:

$\begin{array}{c}{\Vert A-A_1\Vert }_F={\Vert \Sigma -{\Sigma }_1\Vert }_F\\ =\sqrt{{\displaystyle \sum }_{i=k+1}^{k_1}{\left({\sigma }_i-{\sigma }_k\right)}^2+{\displaystyle \sum }_{i=k_1+1}^{\ell }{\sigma }_i^2}\\ \le \sqrt{\ell -k_1}\frac{{\sigma }_k}{2},\end{array}$

$\begin{array}{c}{\Vert A-A_2\Vert }_F={\Vert \Sigma -{\Sigma }_2\Vert }_F\\ =\sqrt{{\displaystyle \sum }_{i=k+1}^{k_2}{\left({\sigma }_i-\sqrt{\frac{{\sigma }_k^2+{\mu }^2}{2}}\right)}^2+{\displaystyle \sum }_{i=k_2+1}^{\ell }{\sigma }_i^2}\\ \le \sqrt{\ell -k_2}\frac{\sqrt{\left({\sigma }_k^2+{\mu }^2\right)/2}}{2},\end{array}$

$\begin{array}{c}{\Vert A-A_3\Vert }_F={\Vert \Sigma -{\Sigma }_3\Vert }_F\\ =\sqrt{{\displaystyle \sum }_{i=k+1}^{k_3}{\left({\sigma }_i-\frac{{\sigma }_k+\mu }{2}\right)}^2+{\displaystyle \sum }_{i=k_3+1}^{\ell }{\sigma }_i^2}\\ \le \sqrt{\ell -k_3}\frac{{\sigma }_k+\mu }{4},\end{array}$

$\begin{array}{c}{\Vert A-A_4\Vert }_F={\Vert \Sigma -{\Sigma }_4\Vert }_F\\ =\sqrt{{\displaystyle \sum }_{i=k+1}^{k_4}{\left({\sigma }_i-\mu \right)}^2+{\displaystyle \sum }_{i=k_4+1}^{\ell }{\sigma }_i^2}\\ \le \sqrt{\ell -k_4}\frac{\mu }{2}.\end{array}$

From (11), the relationship between the Frobenius norms of matrix $A$ and the approximate matrices is:

${\Vert A-A_1\Vert }_F\ge {\Vert A-A_2\Vert }_F\ge {\Vert A-A_3\Vert }_F\ge {\Vert A-A_4\Vert }_F.$

The spectral condition numbers of the approximate matrices $A_i,i=1,2,3,4$ are:

${\kappa }_2\left(A_1\right)=\frac{{\sigma }_1}{{\sigma }_k},{\kappa }_2\left(A_2\right)=\frac{\sqrt{2}{\sigma }_1}{\sqrt{{\sigma }_k^2+{\mu }^2}},$

${\kappa }_2\left(A_3\right)=\frac{2{\sigma }_1}{{\sigma }_k+\mu },{\kappa }_2\left(A_4\right)=\frac{{\sigma }_1}{\mu }.$

The relationship between the spectral condition numbers of the four approximate matrices is:

${\kappa }_2\left(A_1\right)\le {\kappa }_2\left(A_2\right)\le {\kappa }_2\left(A_3\right)\le {\kappa }_2\left(A_4\right).$

2.2. Filter Representation of the Solution

The filtered solution of NMTSVD is:

$x_{filt}^i={\displaystyle \sum }_{j=1}^{\ell }{\varphi }_i\frac{u_j^{\text{T}}b}{{\sigma }_j}{v}_j,i=1,2,3,4.$

The corresponding four filter factors ${\varphi }_i,i=1,2,3,4$ are:

${\varphi }_1=\{\begin{array}{ll}1,\hfill & 1\le j\le k,\hfill \\ \frac{{\sigma }_j}{{\sigma }_k},\hfill & k<j\le k_1,\hfill \\ 0,\hfill & k_1<j\le \ell .\hfill \end{array}$

${\varphi }_2=\{\begin{array}{ll}1,\hfill & 1\le j\le k,\hfill \\ \frac{{\sigma }_j}{\sqrt{\left({\sigma }_k^2+{\mu }^2\right)/2}},\hfill & k<j\le k_2,\hfill \\ 0,\hfill & k_2<j\le \ell .\hfill \end{array}$

${\varphi }_3=\{\begin{array}{ll}1,\hfill & 1\le j\le k,\hfill \\ \frac{{\sigma }_j}{\left({\sigma }_k+\mu \right)/2},\hfill & k<j\le k_3,\hfill \\ 0,\hfill & k_3<j\le \ell .\hfill \end{array}$

${\varphi }_4=\{\begin{array}{ll}1,\hfill & 1\le j\le k,\hfill \\ \frac{{\sigma }_j}{\mu },\hfill & k<j\le k_4,\hfill \\ 0,\hfill & k_4<j\le \ell .\hfill \end{array}$ (13)

From the relationship between the four filter factors and the four modified partial singular values (11), the following Proposition 2.1 is obtained.

Proposition 1 The four filter factors ${\varphi }_i,i=1,2,3,4$ are defined by (13), respectively, for $1\le j\le \ell$ . Then,

${\varphi }_1\le {\varphi }_2\le {\varphi }_3\le {\varphi }_4,1\le j\le \ell ,$

And,

$\frac{1}{2}\le {\varphi }_j^i\le 1,k<j\le k_i,i=1,2,3,4.$

Proof. The inequalities follow from (9), (11) and (13).

3. Truncated Reorthogonalized Golub-Kahan Bidiagonalization (RGKB)

In this paper, we first propose the Gram-Schmidt reorthogonalization process to reorthogonalize $V_k$ , which is summarized as function 1; $U_k$ is still calculated using the GKB recurrence formula (14), and the orthogonality of $U_k$ will not be lost, which is summarized as function 2. Second, in order not to increase too much computational complexity, the iteration truncate the dimension $\overline{k}$ of the Krylov subspace, which is summarized as Algorithm 1. Finally, Algorithm 2 for the truncated reorthogonalized Golub-Kahan bidiagonalization method (RGKB) is presented, and the regularization parameter $\mu$ is obtained by minimizing the GCV function.

The recursion formulas for the GKB decomposition of $A$ with the initial vector $b$ are given by,

${\beta }_1{u}_1=b,{\alpha }_1{v}_1=A^{\text{T}}{u}_1;$

${\beta }_{k+1}{u}_{k+1}=A{v}_k-{\alpha }_k{u}_k,{\alpha }_{k+1}{v}_{k+1}=A^{\text{T}}{u}_{k+1}-{\beta }_{k+1}{v}_k.$

The $k$ steps of Golub-Kahan bidiagonalization to $A$ with the initial vector $b$ givens the decompositions,

$A{V}_k=U_{k+1}{\tilde{C}}_k,A^{\text{T}}{U}_k=V_k{C}_k^{\text{T}},U_{k+1}{e}_1=b/\Vert b\Vert ,$ (14)

where the matrices $U_k=\left[u_1,u_2,\cdot \cdot \cdot ,u_k\right]\in {ℝ}^{m\times k}$ , $U_{k+1}=\left[U_k,u_{k+1}\right]\in {ℝ}^{m\times \left(k+1\right)}$ , and $V_k=\left[v_1,v_2,\cdot \cdot \cdot ,v_k\right]\in {ℝ}^{n\times k}$ have orthogonal columns.

$C_k=\left[\begin{array}{lllll}{\alpha }_1\hfill & \hfill & \hfill & \hfill & 0\hfill \\ {\beta }_2\hfill & {\alpha }_2\hfill & \hfill & \hfill & \hfill \\ \hfill & \ddots \hfill & \ddots \hfill & \hfill & \hfill \\ \hfill & \hfill & {\beta }_{k-1}\hfill & {\alpha }_{k-1}\hfill & \hfill \\ 0\hfill & \hfill & \hfill & {\beta }_k\hfill & {\alpha }_k\hfill \end{array}\right]\in {ℝ}^{k\times k}$

is lower bidiagonal matrix, and ${\tilde{C}}_k={\left[C_k^{\text{T}},{\beta }_{k+1}{e}_k\right]}^{\text{T}}$. Here,

$e_j={\left[0,\cdot \cdot \cdot ,0,1,0,\cdot \cdot \cdot ,0\right]}^{\text{T}}$ is the $j$ -th axis vector. The columns of the matrix $V_k$ span the Krylov subspace (7).

Introduce the QR factorization of $L{V}_k$ ,

$L{V}_k=Q_k{R}_k,$ (15)

where $Q_k\in {ℝ}^{p\times k}$ has orthogonal columns and $R_k\in {ℝ}^{k\times k}$ is upper triangular.

From the GKB decomposition (14) of $A$ and the QR decomposition (15) of $L{V}_k$ , we have

$A^{\text{T}}A=V_k{\tilde{C}}_k^{\text{T}}{U}_{k+1}^{\text{T}}{U}_{k+1}{\tilde{C}}_k{V}_k^{\text{T}}=V_k{\tilde{C}}_k^{\text{T}}{\tilde{C}}_k{V}_k^{\text{T}}.$

${\left(L{V}_k\right)}^{\text{T}}\left(L{V}_k\right)={\left(Q_k{R}_k\right)}^{\text{T}}\left(Q_k{R}_k\right)=R_k^{\text{T}}{R}_k.$

The influence matrix is given by $A^{♯}={\left(A^{\text{T}}A+{\mu }_k^{-2}{L}^{\text{T}}L\right)}^{-1}{A}^{\text{T}}$ . Let $x_{\mu }=V_k{y}_k$ , where $y_k\in {ℝ}^k$ , then the Generalized Cross-Validation (GCV) function is defined as follows:

$\text{GCV}\left(\mu \right):=\frac{{\Vert b-A{x}_{\mu }\Vert }^2}{{\left[trace\left(I-A{A}^{♯}\right)\right]}^2},$

where the numerator part is

$\begin{array}{c}{\Vert b-A{x}_{\mu }\Vert }^2={\Vert b-A{V}_k{y}_k\Vert }^2\\ ={\Vert b-U_{k+1}{\tilde{C}}_k{y}_k\Vert }^2\\ ={\Vert U_{k+1}^{\text{T}}b-{\tilde{C}}_k{y}_k\Vert }^2\\ ={\Vert \Vert b\Vert e_1-{\tilde{C}}_k{y}_k\Vert }^2.\end{array}$

The matrix part in the denominator is:

$\begin{array}{c}I-A{A}^{♯}=I-A{\left(A^{\text{T}}A+{\mu }_k^{-2}{L}^{\text{T}}L\right)}^{-1}{A}^{\text{T}}\\ =I-A{V}_k{\left({\left(A{V}_k\right)}^{\text{T}}A{V}_k+{\mu }_k^{-2}{\left(L{V}_k\right)}^{\text{T}}\left(L{V}_k\right)\right)}^{-1}{V}_k^{\text{T}}{A}^{\text{T}}\\ =I-U_{k+1}{\tilde{C}}_k{\left(V_k{\tilde{C}}_k^{\text{T}}{\tilde{C}}_k{V}_k^{\text{T}}+{\mu }_k^{-2}{R}_k^{\text{T}}{R}_k\right)}^{-1}{\tilde{C}}_k^{\text{T}}{U}_{k+1}^{\text{T}},\end{array}$

where $r_k=b-A{x}_k$ , $r_{k-1}=b-A{x}_{k-1}$ , and $tol$ is a user-defined very small error margin.

Left-multiplying both sides of equation (6) by the matrix $V_k^{\text{T}}$ , we have

$\left({\left(A{V}_k\right)}^{\text{T}}\left(A{V}_k\right)+{\mu }_k^{-2}{\left(L{V}_k\right)}^{\text{T}}\left(L{V}_k\right)\right)y_k={\mu }_k^{-2}{\left(L{V}_k\right)}^{\text{T}}L{V}_{k-1}{y}_{k-1}+{\left(A{V}_k\right)}^{\text{T}}b.$ (16)

Since $A$ and $L$ satisfy (4), the matrix $\left[A{V}_k,L{V}_k\right]$ is of full rank. For any $0<{\mu }_k<\infty$ , the matrix on the left-hand side of (16) is nonsingular. Therefore,

$\left({\tilde{C}}_k^{\text{T}}{\tilde{C}}_k+{\mu }_k^{-2}{R}_k^{\text{T}}{R}_k\right)y_k={\mu }_k^{-2}{R}_k^{\text{T}}{R}_k{y}_{k-1}+{\tilde{C}}_k^{\text{T}}\Vert b\Vert e_1$ (17)

has a unique solution $y_k$ .

The termination condition for the iteration is:

$\frac{\Vert r_k-r_{k-1}\Vert }{\Vert b\Vert }\le tol.$ (18)

As the number of iteration steps $k$ increases, the orthogonal matrices $U_k$ and $V_k$ in (14) tend to lose their orthogonality. To ensure the orthogonality of $U_k$ and $V_k$ , the Gram-Schmidt reorthogonalization process is used to reorthogonalize $V_k$ . $U_k$ is still calculated using the recursive formula (14), and its orthogonality will not be lost. For the detailed proof of the reorthogonalization process, refer to Barlow’s paper [28].

Let $r_k=A^{\text{T}}{u}_{k+1}-{\beta }_{k+1}{v}_k$ , and the Gram-Schmidt reorthogonalization process is summarized into the following subfunction GS_reorth.

function 1: $\left[v_{k+1},h,{\alpha }_{k+1}\right]=\text{GS\_reorth}\left(V_k,r_k,setzero\right)$

$h^1=V_k^{\text{T}}{r}_k;r_k^1=r_k-V_k{h}^1$ ;

if ${\Vert r_k^1\Vert }_2\ge \sqrt{0.8}{\Vert r_k\Vert }_2$

${\alpha }_{k+1}={\Vert r_k^1\Vert }_2;v_{k+1}=r_k^1/{\alpha }_{k+1};h_k=h^1$ ;

else

$h^2=V_k^T{r}_k^1;r_k^2=r_k^1-V_k{h}^2$ ;

$h=h^1+h^2$ ;

if ${\Vert r_k^2\Vert }_2\ge \sqrt{0.8}{\Vert r_k^1\Vert }_2$

${\alpha }_{k+1}={\Vert r_k^2\Vert }_2;v_{k+1}=r_k^2/{\alpha }_{k+1}$ ;

else

Find $e_J$ such that

${\Vert V_k^{\text{T}}{e}_J\Vert }_2=\underset{1\le j\le n}{\text{min}}{\Vert V_k^{\text{T}}{e}_j\Vert }_2$ .

$c^1=V_k^{\text{T}}{e}_J;t^1=e_J-V_k{c}^1$ ;

$c^2=V_k^{\text{T}}{t}^1;t^2=t^1-V_k{c}^2$ ;

$V_{k+1}=t^2/{\Vert t^2\Vert }_2$ ;

if setzero

${\alpha }_{k+1}=0$

else

${\alpha }_{k+1}=v_{k+1}^{\text{T}}{r}_k^2$

end

end

end

end GS_reorth

First, the GKB method with the GS process is presented and summarized as the subfunction GKB_reorth.

function 2 $\left[U,C,V\right]=\text{GKB\_reorth}\left(A,b,k,setzero\right)$

$U_1=b/\Vert b\Vert$ ;

$C_{1,1}={\Vert A^{\text{T}}{U}_1\Vert }_2;V_1=A^{\text{T}}{U}_1/C_{1,1}$ ;

for $j=1,\cdots ,k$

$\beta =A{V}_j-C_{j,j}{U}_j$ ;

$C_{j+1,j}=\Vert \beta \Vert$ ;

$U_{j+1}=\beta /C_{j+1,j}$ ;

$\alpha =A^{\text{T}}{U}_{j+1}-C_{j+1,j}{V}_j$ ;

$\left[V_{j+1},h,C_{j+1,j+1}\right]=\text{GS\_reorth}\left(V_j,\alpha ,setzero\right)$ ;

end

end GKB_reorth.

Next, to ensure the orthogonality of the orthogonal matrices $U_k$ and $V_k$ without adding excessive computational complexity, this paper truncates the iteration process at subscript $\overline{k}$ , meaning the dimension of the solution space is $\overline{k}$ . This is summarized as Algorithm 0. Finally, the process of the RGKB method proposed in this paper is summarized as Algorithm 0.

RGKB method proposed in this paper is summarized as Algorithm 2.

Ultimately, only the $\overline{k}$ -dimensional solution space $V_{\overline{k}}$ is considered. Due to the semi-convergence property of iterative methods for large-scale discrete ill-posed problems, the convergence within the $\overline{k}$ -dimensional solution space is excellent, while the convergence beyond this space deteriorates. This phenomenon will be demonstrated in the experimental section of Section 4.2.

4. Numerical Experiments

This chapter presents two experiments to verify the effectiveness of the proposed NMTSVD methods. The experiments are conducted using MATLAB R2024b on an INTEL(R) Core(TM) i3-2310M CPU (2.1GHz) with 2 GB RAM.

We use examples from the MATLAB package Regularization Tools [31]. The error $e$ in the vector $b$ is Gaussian white noise, with the noise level defined as:

$\theta :=\frac{\Vert e\Vert }{\Vert b_{true}\Vert }.$

The regularization matrix is,

$L_1=\left[\begin{array}{llll}1\hfill & -1\hfill & \hfill & \hfill \\ \hfill & \ddots \hfill & \ddots \hfill & \hfill \\ \hfill & \hfill & 1\hfill & -1\hfill \end{array}\right]\in {ℝ}^{\left(n-1\right)\times n}.$

The measure of approximation between the approximate solution $x$ and the true solution $x_{true}$ is:

$err:=\frac{\Vert x-x_{true}\Vert }{\Vert x_{true}\Vert }.$

The free parameter $tol$ in the termination condition (18) is set to 1 × 10⁻⁶.

4.1. The Experimental Results of NMTSVD

Example 0.1: The test problems use the code shaw, gravity, foxgood, and deriv2 from [31], with the problem dimension set as $n=200$ , $m=n$ . The regularization parameter is solved by the discrepancy principle with $\eta =1.01$ . Under different noise levels $\theta =1\times {10}^{-2},1\times {10}^{-3},1\times {10}^{-4}$ , the truncation indices $k$ of TSVD [5], $k_1$ of MTSVD [16], and $k_2,k_3,k_4$ of NMTSVD methods determined by the four problems are given below.

First, Table 1 presents the values of the truncation indices $k,k_1,k_2,k_3,k_4$ and their relationships for several noise levels. In most cases, there is $k\le k_1=k_2=k_3<k_4$ , while in individual cases, there is $k=k_1=k_2=k_3=k_4$ .

Second, Figures 1-4 display the singular value distributions of the four problems within a specific index.

Table 1. Truncation indices and relationship 0.1.

Problems	$\theta$	$k$	$k_1$	$k_2$	$k_3$	$k_4$	Relation
3*shaw	1⋅10−2	5	6	6	6	7	$k<k_1=k_2=k_3<k_4$
	1⋅10−3	7	7	7	7	7	$k=k_1=k_2=k_3=k_4$
	1⋅10−4	8	8	8	8	9	$k=k_1=k_2=k_3<k_4$
3*gravity	1⋅10−2	7	8	8	8	8	$k<k_1=k_2=k_3=k_4$
	1⋅10−3	8	9	9	9	10	$k<k_1=k_2=k_3<k_4$
	1⋅10−4	10	11	11	11	12	$k<k_1=k_2=k_3<k_4$
3*foxgood	1⋅10−2	2	2	2	2	2	$k=k_1=k_2=k_3=k_4$
	1⋅10−3	2	3	3	3	3	$k<k_1=k_2=k_3=k_4$
	1⋅10−4	3	3	3	3	4	$k=k_1=k_2=k_3<k_4$
3*deriv2	1⋅10−2	8	11	11	11	12	$k<k_1=k_2=k_3<k_4$
	1⋅10−3	18	25	25	25	26	$k<k_1=k_2=k_3<k_4$
	1⋅10−4	39	54	54	54	55	$k<k_1=k_2=k_3<k_4$

Figure 1. Shaw, the singular value of No.4-12.

Figure 2. Gravity, the singular value of No.6-13.

Figure 3. Foxgood, the singular value of No.1-6.

Figure 4. Deriv2, the singular value of No.8-55.

Finally, Table 2 presents the comparative experiment results of TSVD Method, MTSVD Method and NMTSVD Method, including the relative errors of the approximate solutions obtained by each method for the four problems under different noise levels $\theta =1\times {10}^{-2},1\times {10}^{-3},1\times {10}^{-4}$ , as well as the running time under specific noise levels $\theta =1\times {10}^{-4}$ . For the four problems under the same noise level, in most cases, the NMTSVD Method yield better approximate solutions than TSVD and MTSVD. The data for foxgood and deriv2 in Table 2 indicate that the NMTSVD methods produce close-to-optimal approximations, with the third NMTSVD method performing particularly better. This is because the singular values of foxgood and deriv2 decay gently (as shown in Figure 3 and Figure 4), and the differences between modified partial singular values are minimal, leading to similar approximation accuracies. The data for shaw and gravity show that the three NMTSVD methods can either obtain better approximations or perform slightly worse than MTSVD but still outperform TSVD. This is attributed to the steep decay of singular values in shaw and gravity (Figure 1 and Figure 2), where smaller partial singular values with $\mu$ may not necessarily guarantee optimal solutions.

Table 2. Errors and Running time under the TSVD, MTSVD, and NMTSVD methods 0.1.

Problems	$\theta$	TSVD	MTSVD	NMTSVD_1	NMTSVD_2	NMTSVD_3
6*shaw	1⋅10−2	0.14619	0.12142	0.11906	0.11622	0.09992
	1⋅10−3	0.04835	0.04835	0.04835	0.04835	0.04910
	1⋅10−4	0.04708	0.04708	0.04708	0.04708	0.03908
	time (10−2)	0.01479	0.01474	0.01474	0.01474	0.01474
	time (10−3)	0.00059	0.00055	0.00055	0.00055	0.00055
	time (10−4)	0.00066	0.00062	0.00062	0.00062	0.00062
6*gravity	1⋅10−2	0.03775	0.03475	0.03486	0.03488	0.03575
	1⋅10−3	0.01762	0.01506	0.01492	0.01490	0.01492
	1⋅10−4	0.00851	0.00700	0.00688	0.00686	0.00671
	time (10−2)	0.00978	0.00974	0.00974	0.00974	0.00974
	time (10−3)	0.00065	0.00061	0.00061	0.00061	0.00061
	time (10−4)	0.00057	0.00053	0.00053	0.00053	0.00053
6*foxgood	1⋅10−2	0.03148	0.03148	0.03148	0.03148	0.03142
	1⋅10−3	0.01072	0.01072	0.01072	0.01072	0.01038
	1⋅10−4	0.00627	0.00627	0.00627	0.00627	0.00457
	time (10−2)	0.01182	0.01178	0.01177	0.01177	0.01177
	time (10−3)	0.00066	0.00062	0.00062	0.00062	0.00062
	time (10−4)	0.00057	0.00053	0.00053	0.00053	0.00053
6*deriv2	1⋅10−2	0.27058	0.24832	0.24741	0.24741	0.24636
	1⋅10−3	0.18437	0.16796	0.16754	0.16753	0.16736
	1⋅10−4	0.12440	0.11242	0.11240	0.11240	0.11235
	time (10−2)	0.01133	0.01129	0.01129	0.01129	0.01129
	time (10−3)	0.00062	0.00059	0.00059	0.00059	0.00059
	time (10−4)	0.00064	0.00061	0.00061	0.00061	0.00061

Collectively, Examples 0.1 demonstrate that for solving small-to-medium sized discrete ill-posed problems, compared with TSVD and MTSVD, we proposed NMTSVD methods achieve superior approximate solutions with fewer iterations, especially when singular values decay gently. The third modified truncated singular value decomposition method is particularly effective. Additionally, the smaller the noise level, the better the approximation accuracy of NMTSVD.

4.2. The Experimental Results of RGKB

Example 0.2: The test problems uses the code shaw, phillips, foxgood, and gravity from [31] with dimensions $n=1000$ , $m=n$ . The noise level is set to $\theta =1\times {10}^{-3}$ , and the parameter $setzero=false$ . For these four problems, the GKB [28] and RGKB algorithms presented in Section 3 are applied to compare the dimension of the final solution space (dimension of the Krylov subspace), the approximation accuracy of the approximate solutions, and the approximation accuracy of the singular value matrices. For each problem under the same noise level, 1000 trials are conducted to compute the average values.

First, Table 3 presents the dimension ($k_0$ and $k$ ) of the Krylov subspace, the relative errors (GKB-S and RGKB-S) of the solution, the running time (GKB-T and RGKB-T), and the approximation accuracy (GKB-M and RGKB-M) of the coefficient matrix for the GKB and RGKB methods. We can observe that RGKB outperforms GKB in terms of efficiency. When the dimension $\ell$ of the initial solution space is known, the final solution space dimension $k$ obtained by RGKB is slightly larger than $k_0$ obtained by GKB, i.e., $k_0<k$ . The approximation accuracy of the solution (RGKB-S) is higher than that of GKB-S, the running time of the solution (RGKB-S) is shorter than that of GKB-S, and the approximation accuracy of the coefficient matrix (RGKB-M) is higher than that of GKB-M. Even with occasional loss of orthogonality, RGKB achieves better approximations for both solutions and matrices with minimal reorthogonalization efforts.

Second, Table 4 verifies the effectiveness of the dimension of the Krylov subspace Algorithm 1 in the RGKB method. The relative errors of the approximate solutions for each problem under different dimensions $k,k+1,k+2,k+3$ are presented. It is found that, for each problem, the RGKB method yields the best approximation of the solution with the dimension $k$ of the Krylov subspace, which also demonstrates the effectiveness of Algorithm 1.

Finally, Figures 5-8 show a comparison of the approximation performance of the approximate solutions for each problem between the GKB and RGKB methods. All results indicate that the RGKB method yields better approximate solutions.

Table 3. The dimension of the Krylov subspace, the relative errors of the solution, the running time , and the approximation accuracy of the coefficient matrix for the GKB and RGKB methods 0.2.

Problems	$\ell$	$k_0$	$k$	GKB-E	RGKB-E	GKB-T	RGKB-T	GKB-M	RGKB-M
shaw	6	6	8	0.0593	0.0477	0.1485	0.0786	1.4350	0.0013
phillips	6	9	10	0.0082	0.0069	0.1464	0.0787	0.3920	0.1305

Continued

foxgood	2	3	4	0.0114	0.0087	0.1470	0.0787	0.0841	0.0002
gravity	6	7	11	0.0319	0.0189	0.1376	0.0748	0.1317	0.0123

Table 4. The effectiveness of dimension of the Krylov subspace Algorithm 1 in the RGKB method.

Problems	$\theta$	$k_0$	$k$	$k+1$	$k+2$	$k+3$
shaw	err	0.0606	0.0476	0.0477	0.0477	0.0477
phillips	err	0.0086	0.0068	0.0114	0.0135	0.0223
foxgood	err	0.0102	0.0085	0.0089	0.0090	0.0089
gravity	err	0.0302	0.0192	0.0265	0.0325	0.0365

Figure 5. Shaw.

Figure 6. Phillips.

Figure 7. Foxgood.

Figure 8. Gravity.

Experimental Examples 0.2 show that: for solving large-scale unstable discrete ill-posed problems, the RGKB proposed in this paper can obtain better approximate solutions compared with GKB. The dimension of the solution space with the best convergence is the dimension $k$ of the solution space determined by RGKB. In solution spaces with dimensions larger or smaller than $k$ , better approximate solutions cannot be obtained.

Example 0.3: We employ the 256 × 256 Cameraman benchmark image as our experimental test subject. First, we corrupt the image with additive Gaussian white noise characterized by a noise level of $\sigma =0.05$ to emulate realistic image degradation in practical scenarios. We then configure the algorithmic parameter $setzero=false$ , and proceed to apply both the GKB and RGKB methods to perform deblurring on the noisy, blurred image. Finally, we conduct an objective quantitative evaluation of the two methods performance using the Peak Signal-to-Noise Ratio (PSNR) as the primary quality assessment metric.

Figure 9. Cameraman.

Table 5. PSNR of deblurred images by GKB and RGKB.

Method	Noisy Image	GKB Reconstructed Image	RGKB Reconstructed Image
PSNR(dB)	26.04	19.74	26.43

Figure 9 presents the visual deblurring results of the GKB and RGKB methods on the 256 × 256 Cameraman benchmark image. The original image (top-left) shows clear details of the cameraman and background architecture. The noisy image (top-right) exhibits visible Gaussian noise, which obscures fine textures. Quantitative results in Table 5 demonstrate that: The GKB-reconstructed image (bottom-left) suffers from severe over-smoothing, with a significant loss of edge details, leading to a PSNR drop to 19.74 dB. The RGKB-reconstructed image (bottom-right) effectively suppresses noise while preserving critical structural details, achieving a PSNR of 26.43 dB, which is 0.41 dB higher than that of the noisy image.

These results demonstrate that the RGKB method outperforms GKB by a substantial margin, with a PSNR improvement of 6.69 dB. This confirms that RGKB achieves a better balance between noise suppression and detail preservation, whereas GKB introduces excessive smoothing and information loss.

Acknowledgements

Sincere thanks to the members of JAMP for their professional performance, and special thanks to managing editor Hellen XU for a rare attitude of high quality.

Funding

This work was supported by the Gansu Province Innovation Fund for College Teachers (Grant No. 2024A-173).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Björck, Å. (1988) A Bidiagonalization Algorithm for Solving Large and Sparse Ill-Posed Systems of Linear Equations. BIT, 28, 659-670.[CrossRef]
[2]	Lawson, C.L. and Hanson, R.J. (1974) Solving Least Squares Problems, Classics in Applied Mathematics, Series Number 15. Prentice-Hall Inc.
[3]	Noschese, S. and Reichel, L. (2016) Some Matrix Nearness Problems Suggested by Tikhonov Regularization. Linear Algebra and its Applications, 502, 366-386.[CrossRef]
[4]	Golub, G.H., Hansen, P.C. and O'Leary, D.P. (1999) Tikhonov Regularization and Total Least Squares. SIAM Journal on Matrix Analysis and Applications, 21, 185-194.[CrossRef]
[5]	Hansen, P.C. (1990) Truncated Singular Value Decomposition Solutions to Discrete Ill-Posed Problems with Ill-Determined Numerical Rank. SIAM Journal on Scientific and Statistical Computing, 11, 503-518.[CrossRef]
[6]	Bauer, F. and Lukas, M.A. (2011) Comparingparameter Choice Methods for Regularization of Ill-Posed Problems. Mathematics and Computers in Simulation, 81, 1795-1841.[CrossRef]
[7]	Calvetti, D., Lewis, B. and Reichel, L. (2002) GMRES, L-Curves, and Discrete Ill-Posed Problems. BIT Numerical Mathematics, 42, 44-65.[CrossRef]
[8]	Golub, G.H., Heath, M. and Wahba, G. (1979) Generalized Cross-Validation as a Method for Choosing a Good Ridge Parameter. Technometrics, 21, 215-223.[CrossRef]
[9]	Kilmer, M.E. and O'Leary., D.P. (2001) Choosing Regularization Parameters in Iterative Methods for Ill-Posed Problems. SIAM Journal on Matrix Analysis and Applications, 22, 1204-1221.[CrossRef]
[10]	Brezinski, C., Redivo-Zaglia, M., Rodriguez, G. and Seatzu, S. (1998) Extrapolation Techniques for Ill-Conditioned Linear Systems. Numerische Mathematik, 81, 1-29.[CrossRef]
[11]	Reichel, L. and Shyshkov, A. (2008) A New Zero-Finder for Tikhonov Regularization. BIT Numerical Mathematics, 48, 627-643.[CrossRef]
[12]	Reichel, L. and Rodriguez, G. (2012) Old and New Parameter Choice Rules for Discrete Ill-Posed Problems. Numerical Algorithms, 63, 65-87.[CrossRef]
[13]	Li, Z., Huang, H. and Wei, Y. (2011) Ill-Conditioning of the Truncated Singular Value Decomposition, Tikhonov Regularization and Their Applications to Numerical Partial Differential Equations. Numerical Linear Algebra with Applications, 18, 205-221.[CrossRef]
[14]	Hansen, P.C. (1990) The Discrete Picard Condition for Discrete Ill-Posed Problems. BIT, 30, 658-672.[CrossRef]
[15]	Morigi, S., Reichel, L. and Sgallari, F. (2007) A Truncated Projected SVD Method for Linear Discrete Ill-Posed Problems. Numerical Algorithms, 43, 197-213.[CrossRef]
[16]	Noschese, S. and Reichel, L. (2014) A Modified Truncated Singular Value Decomposition Method for Discrete Ill-Posed Problems. Numerical Linear Algebra with Applications, 21, 813-822.[CrossRef]
[17]	Jia, Z. and Yang, Y. (2018) Modified Truncated Randomized Singular Value Decomposition (MTRSVD) Algorithms for Large Scale Discrete Ill-Posed Problems with General-Form Regularization. Inverse Problems, 34, Article ID: 055013.[CrossRef]
[18]	Hanke, M. and Groetsch, C.W. (1998) Nonstationary Iterated Tikhonov Regularization. Journal of Optimization Theory and Applications, 98, 37-53.[CrossRef]
[19]	Hansen, P.C. and O’Leary, D.P. (1993) The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems. SIAM Journal on Scientific Computing, 14, 1487-1503.[CrossRef]
[20]	Huang, G., Reichel, L. and Yin, F. (2015) On the Choice of Solution Subspace for Nonstationary Iterated Tikhonov Regularization. Numerical Algorithms, 72, 1043-1063.[CrossRef]
[21]	Novati, P. and Russo, M.R. (2013) A GCV Based Arnoldi-Tikhonov Regularization Method. BIT Numerical Mathematics, 54, 501-521.[CrossRef]
[22]	Gazzola, S., Novati, P. and Russo, M.R. (2015) On Krylov Projection Methods and Tikhonov Regularization. Electronic Transactions on Numerical Analysis, 44, 83-123.
[23]	Gazzola, S. (2014) Regularization Techniques Based on Krylov Methods for Ill-Posed Linear Systems. Ph.D. Thesis, University of Padua.
[24]	Larsen, R.M. (1998) Lanczos Bidiagonalization with Partial Reorthogonalization. DAIMI Report Series, 27, 1-101.[CrossRef]
[25]	Daniel, J.W., Gragg, W.B., Kaufman, L. and Stewart, G.W. (1976) Reorthogonalization and Stable Algorithms for Updating the Gram-Schmidt QR Factorization. Mathematics of Computation, 30, 772-795.[CrossRef]
[26]	Paige, C.C. (1971) The Computation of Eigenvalues and Eigenvectors of Very Large Sparse Matrices. Ph.D. Thesis, University of London.
[27]	Parlett, B.N. and Scott, D.S. (1979) The Lanczos Algorithm with Selective Orthogonalization. Mathematics of Computation, 33, 217-238.[CrossRef]
[28]	Barlow, J.L. (2013) Reorthogonalization for the Golub-Kahan-Lanczos Bidiagonal Reduction. Numerische Mathematik, 124, 237-278.[CrossRef]
[29]	Alqahtani, A., Ramlau, R. and Reichel, L. (2022) Error Estimates for Golub-Kahan Bidiagonalization with Tikhonov Regularization for Ill-Posed Operator Equations. Inverse Problems, 39, Article ID: 025002.[CrossRef]
[30]	Bianchi, D. and Donatelli, M. and Furchi, D. and Reichel, L. (2025) The Iterated Golub-Kahan-Tikhonov Method. arXiv: 2507.12307.
[31]	Hansen, P.C. (2007) Regularization Tools Version 4.0 for Matlab 7.3. Numerical Algorithms, 46, 189-194.[CrossRef]