Two New Efficient Iterative Regularization Methods for Image Restoration Problems

Abstract

Iterative regularization methods are efficient regularization tools for image restoration problems. The IDR(𝑠) and LSMR methods are state-of-the-arts iterative methods for solving large linear systems. Recently, they have attracted considerable attention. Little is known of them as iterative regularization methods for image restoration. In this paper, we study the regularization properties of the IDR(𝑠) and LSMR methods for image restoration problems. Comparative numerical experiments show that IDR(𝑠) can give a satisfactory solution with much less computational cost in some situations than the classic method LSQR when the discrepancy principle is used as a stopping criterion. Compared to LSQR, LSMR usually produces a more accurate solution by using the 𝐿-curve method to choose the regularization parameter.

1. Introduction

Image deblurring is one of the most classic linear inverse problems. Blurring in images can arise from many sources, such as limitations of the optical system, camera and object motion, astigmatism, and environmental effects. The blurring process of an image can be formulated as a Fredholm integral equation of the first kind which has the following classic form: _ℝ2𝐾𝑠−𝑠^′,𝑡−𝑡^′𝑓𝑠^′,𝑡^′𝑑𝑠^′𝑑𝑡^′(),()=𝑔𝑠,𝑡𝑠,𝑡∈ℝ², (1) where 𝑓 and 𝑔 are the original image and the blurred image, respectively. 𝐾 is a given point spread function (PSF). PSF describes the blurring and the resulting image of the point source. The degree of blurring of the point object is a measure for the quality of an imaging system. More information about different PSFs can be found in [1].

Equation (1) can be discretized to form a linear system 𝐴𝑥_true=𝑏_true,𝐴∈ℝ^𝑚2×𝑚2,𝑥,𝑏∈ℝ^𝑚2, (2) where the matrix 𝐴 is ill-conditioned since it has many singular values close to zero [1]. For simplicity, we assume that 𝐴 in this paper is nonsingular.

However, the right-hand side 𝑏 is not available in many practical applications of image restoration because of the contamination of noise, so the linear system (2) can be reformulated as 𝐴𝑥=𝑏,𝑏=𝑏_true+𝑒,𝑒∈ℝ^𝑚2, (3) where 𝑒 represents the noise.

Our goal is to obtain a good approximation of the original image 𝑥_true by solving the system (3) instead of the system (2) since 𝑏_true is not known. However, the solution of (3) is not a good approximation of the solution of (2) because of ill-conditioned 𝐴. In fact, there is a quite remarkable disparity among the corresponding solutions of (3) and (2) even if the norm of 𝑒 is small.

The singular value decomposition (SVD) is a powerful tool for the analysis of discrete ill-posed problems. The SVD takes the form 𝐴=𝑈Σ𝑉^𝑇=^𝑚2_𝑖=1𝑢_𝑖𝜎_𝑖𝑣^𝑇_𝑖, (4) where Σ is a diagonal matrix with the singular values 𝜎₁≥𝜎₂≥⋯≥𝜎_𝑚2≥0. For most image deblurring problems, the discrete Picard condition is satisfied: for all singular values larger than the norm of noise, the corresponding coefficients |𝑢^𝑇_𝑖𝑏| decay faster than the 𝜎_𝑖.

Regularization methods are used to produce better solution. One of the most famous regularization methods is truncated SVD (TSVD) [2], which yields solution 𝑥_𝑘=∑^𝑘_𝑖=1𝜎_𝑖⁻¹(𝑢^𝑇_𝑖𝑏)𝑣_𝑖. Another regularization method is Tikhonov regularization which consists of solving the minimization problem min_{𝑥∈𝑅𝑚2}‖𝐴𝑥−𝑏‖²₂+𝜆‖𝐿𝑥‖²₂, (5) where 𝐿 is a regularization operator, which is often chosen to be the identity matrix ([3, Chapter  5]) and 𝜆 is the regularization parameter that controls the weighting between the two terms. However, the previously mentioned techniques are usually computationally expensive for large-scale problems like image deblurring. Iterative methods, especially Krylov subspace iterative methods, are used to solve these problems due to their inexpensive computational cost and easy implementation (see [4] for more details about Krylov subspace methods).

When we apply a Krylov subspace method to solve (3), the semiconvergence property is observed. Iterative methods can produce a sequence of iterates {𝑥₁,𝑥₂,…,𝑥_𝑘,…} that initially tends to get closer to the exact solution but diverges again from the exact solution in later stages because the influence of the noise starts to dominate the solution. In this situation, the iteration number 𝑘 can be considered as a regularization parameter. We also introduce different methods for choosing an effective 𝑘 in this paper.

The conjugate gradient (CG) method is a classical Krylov subspace iterative method for solving the linear systems with symmetric positive definite (SPD) matrix, and its regularizing effects are well known (see [5] and [3, Chapter  6] for more details). When 𝐴 is not SPD, the method can be applied to the normal equations 𝐴^𝑇𝐴𝑥=𝐴^𝑇𝑏. CGLS is the most stable way to implementation of CG algorithm for the normal equations ([3, Chapter  6], [6, Chapter  6]). Other Krylov subspace methods also have been used to solve image restoration problems. The regularizing properties of GMRES have been studied in [7]. The regularizing properties of BiCG and QMR have been tested in [8], which shows that BiCG and QMR are much faster than CGLS. The comparison of the regularizing properties of Krylov subspace methods CGLS, GMRES, CGS, BiCG, QMR, and BiCGSTAB has been investigated in [9]: CGS, BiCG, and BiCGSTAB are not sufficiently consistent with the discrepancy principle (we say that a method is consistent with the discrepancy principle if the residual norm with discrepancy principle as the stopping rule is linear dependence to the noise norm [9]. See Section 3.1 for details of the discrepancy principle); CGLS method is both efficient and consistent with the discrepancy principle, but it needs a large number of iterations; GMRES and QMR show intermediate features.

Recently, IDR(𝑠) [10] and LSMR [11] iterative methods have been proposed for solving large linear system. IDR(𝑠) is a new family of iterative methods which is fast and requires low memory. IDR(𝑠) has attracted considerable attention, and different variants have been proposed in [12, 13]. LSMR is similar to LSQR by applying the MINERS method to the normal equation. The regularizing effects of MINERS have been analysed in [14, 15]. It has been proved that both ‖𝐴𝑟_𝑘‖ and ‖𝑟_𝑘‖ decrease monotonically for LSMR, while only ‖𝑟_𝑘‖ is monotonic for LSQR. So, it is safer to terminate LSMR early.

Many methods have been proposed for choosing the regularization parameter. But a robust method appropriate for different situations is to be found. Two methods are used in this paper: the discrepancy principle (DP) and the discrete 𝐿-curve criterion. DP is the most famous and frequently used method. It has distinct advantages and disadvantages: it is simple and easy to be utilized, but it needs to know the norm of the noise in advance. The discrete 𝐿-curve criterion is used to choose the regularization parameter by using the curve of (log‖𝑟_𝑘‖,log‖𝑥_𝑘‖). For more information about methods for choosing the regularization parameter, we suggest [6, Chapter  5] and references therein.

Inspired by the advantages of IDR(𝑠) and LSMR, we apply them to solve the image restoration problems. We compare the performance of the IDR(𝑠) method with that of the LSQR method for solving (3) with the iteration terminated by the DP method. We also make the comparison between the LSMR method and the LSQR method when the discrete 𝐿-curve criterion is used to choose the regularization parameter.

The paper is organized as follows. In Section 2, we review the IDR(𝑠) and LSMR algorithms briefly and present analysis of the regularizing properties of the IDR(𝑠) and LSMR algorithms. In Section 3, methods for choosing the regularization parameter especially the DP method and the discrete 𝐿-curve criterion are introduced. IDR(𝑠) and LSMR methods are applied to solve image restoration problems and compared with with LSQR in Section 4. Section 5 summarizes conclusions.

2. IDR(𝑠) and LSMR

2.1. The IDR(𝑠) Method

The IDR(𝑠) methods are a new family of iterative methods based on the induced dimension reduction (IDR) method proposed by Sonneveld in 1980. IDR(𝑠) algorithm can generate a sequence of nested subspaces that corresponding residuals are forced to be in. The dimension of these nested subspaces is monotonically decreasing, and so in exact arithmetic, the solution can be obtained within a not very large number of iterations.

The nested subspaces are defined by 𝒢_𝑗=𝐼−𝜔_𝑗𝐴𝒢_𝑗−1∩𝒮, (6) where 𝒮 is a fixed proper subspace of ℝ^𝑁 and 𝜔_𝑗 are nonzero scalars. The space 𝒮 is usually assumed as the left null space of 𝑃∈ℝ^𝑁×𝑠: 𝑃=𝑝₁𝑝₂…𝑝_𝑠,𝒮=𝒩𝑃^𝐻. (7)

The columns of 𝑃 are chosen as orthogonalization of a set of random vectors.

Suppose that after 𝑛 iterations we have 𝑠 vectors 𝑔_𝑖∈𝒢_𝑗, 𝑖=𝑛−1,…,𝑛−𝑠 and 𝑠 corresponding vectors 𝑢_𝑖 with 𝑔_𝑖=𝐴𝑢_𝑖. Define the matrices 𝐺_𝑛=𝑔_𝑛−𝑠𝑔_{𝑛−𝑠+1}⋯𝑔_𝑛−1,𝑈_𝑛=𝑢_𝑛−𝑠𝑢_{𝑛−𝑠+1}⋯𝑢_𝑛−1. (8)

The algorithm produces 𝑠+1 residuals in 𝒢_𝑗; then the next residual can be found in 𝒢_𝑗+1. Two steps are needed to obtain the 𝑠+1 residuals.

Step  1 Compute the first residual in 𝒢_𝑗+1. Calculate 𝑐 from (𝑃^𝐻𝐺_𝑛)𝑐=𝑃^𝐻𝑟_𝑛 𝑣_𝑛=𝑟_𝑛−𝐺_𝑛𝑐 𝑟_𝑛+1=𝑟_𝑛−𝜔_𝑗+1𝑣_𝑛−𝐺_𝑛𝑐 𝑥_𝑛+1=𝑥_𝑛+𝜔_𝑗+1𝐴𝑣_𝑛+𝑈_𝑛𝑐.

Step  2 Compute the remaining residuals in 𝒢_𝑗+1, for 𝑘=1,…𝑠. Calculate 𝑐 from (𝑃^𝐻𝐺_{𝑛+𝑘−1})𝑐=𝑃^𝐻𝑟_𝑛+𝑘 𝑣_𝑛+𝑘=𝑟_𝑛+𝑘−𝐺_{𝑛+𝑘−1}𝑐 𝑢_𝑛+𝑘=𝜔_𝑗+1𝑣_𝑛+𝑘+𝑈_{𝑛+𝑘−1}𝑐 𝑔_𝑛+𝑘=𝐴𝑢_𝑛+𝑘 select 𝛼_𝑖 and 𝛽_𝑖 𝑔_𝑛+𝑘=𝑔_𝑛+𝑘−∑^𝑘−1_𝑖=1𝛼_𝑖𝑔_𝑛+𝑖; 𝑢_𝑛+𝑘=𝑢_𝑛+𝑘−∑^𝑘−1_𝑖=1𝛼_𝑖𝑢_𝑛+𝑖;  𝑟_𝑛+𝑘+1=𝑟_𝑛+𝑘−∑^𝑘_𝑖=1𝛽_𝑖𝑔_𝑛+𝑖; 𝑥_𝑛+𝑘+1=𝑥_𝑛+𝑘+∑^𝑘_𝑖=1𝛽_𝑖𝑢_𝑛+𝑖;  update 𝐺_∶,𝑘=𝑔_𝑛+𝑘 to get 𝐺_𝑛+𝑘+1; 𝑈_∶,𝑘=𝑢_𝑛+𝑘 to get 𝑈_𝑛+𝑘+1.

Different choices of 𝛼_𝑖 and 𝛽_𝑖 would result in different variants of IDR(𝑠). Take the same notations as in [12]: IDR(𝑠)-proto represents the algorithm in [10] and IDR(𝑠)-biortho represents the variant of the IDR(𝑠) in [12]. IDR(𝑠)-proto just computes 𝑟_𝑛+𝑘+1=𝑟_𝑛+𝑘−𝑔_𝑛+𝑘,𝑥_𝑛+𝑘+1=𝑥_𝑛+𝑘+𝑢_𝑢+𝑘, (9) that is, 𝛼_𝑖=0(𝑖=1,…,𝑘−1), 𝛽_𝑖=0(𝑖=1,…,𝑘−1), and 𝛽_𝑖=1(𝑖=𝑘).

IDR(𝑠)-biortho selects 𝛼_𝑖 and 𝛽_𝑖 to construct vectors that satisfy the biorthogonality conditions: 𝑔_𝑛+𝑘⟂𝑝_𝑖𝑟,𝑖=1,…,𝑘−1,𝑘=2,…,𝑠,_𝑛+𝑘+1⟂𝑝_𝑖,𝑖=1,…,𝑘,𝑘=1,…,𝑠. (10) The IDR(𝑠)-biortho algorithm is more stable and more accurate than the IDR(𝑠)-proto algorithm, so we apply the IDR(𝑠)-biortho algorithm to solve the problems in this paper.

As described in [10, 16], we have the following expression after 𝑛 iterations: 𝑟_𝑛=Φ_𝑛(𝐴)𝑟₀, (11)Φ_𝑛(𝐴) is a polynomial of degree 𝑛. Φ_𝑛(𝐴) can be explicitly written as a product of two polynomials Φ_𝑛(𝐴)=Ω_𝑖(𝐴)Ψ_𝑛−𝑖(𝐴), (12) where Ω_𝑖(𝑡)=∏^𝑖_𝑘=1(1−𝜔_𝑘𝑡), Ψ_𝑛−𝑖(𝑡)=∏^𝑛−𝑖_𝑗=0𝑐_𝑗𝑡^𝑗. The factor Ω_𝑖(𝐴) is called damping factor and Ψ_𝑛−𝑖(𝐴) is called Lanczos factor. The damping factor is the residual polynomial of damped Richardson iterative method. Damped Richardson iterative method is used as an iterative regularizing method in [17] and sometimes referred to as Van Citter method in the literature; compare [18]. According to the conclusion of [16], the Lanczos factor is determined completely by the choice of 𝑃, whereas the damping factor depends on the choice of the parameters 𝜔_𝑖. Actually, in our experiments, different choices of parameter 𝑠 (which controls the dimension of 𝑃) did not show notable influence to the performance in solving image restoration problems. So, the damping factor plays a role of iterative regularizing method.

2.2. The LSMR Method

The LSMR algorithm is equivalent to the MINRES applied to the normal equation. The algorithm is implemented by the Lanczos bidiagonalization algorithm. One form of the bidiagonalization procedure is the Golub-Kahan process: 𝛽₁𝑢₁=𝑏,𝛼₁𝑣₁=𝐴^𝑇𝑢₁𝛽for𝑘=1,2,…,_𝑘+1𝑢_𝑘+1=𝐴𝑣_𝑘−𝛼_𝑘𝑢_𝑘,𝛼_𝑘+1𝑣_𝑘+1=𝐴^𝑇𝑢_𝑘−𝛽_𝑘+1𝑣_𝑘. (13) The scalars 𝛼_𝑘≥0 and 𝛽_𝑘≥0 are chosen so that the norms of 𝑢_𝑘 and 𝑣_𝑘 are both equal to one.

We define that 𝑉_𝑘=(𝑣₁𝑣₂⋯𝑣_𝑘), 𝑈_𝑘=(𝑢₁𝑢₂⋯𝑢_𝑘), and 𝐵_𝑘=𝛼₁𝛽₂𝛼₂𝛽⋱⋱_𝑘𝛼_𝑘𝛽_𝑘+1. (14) Then, the following relations hold: 𝑈_𝑘+1𝛽₁𝑒₁=𝑏,𝐴𝑉_𝑘=𝑈_𝑘+1𝐵_𝑘,𝐴^𝑇𝑈_𝑘+1=𝑉_𝑘+1𝐵^𝑇_𝑘𝛼_𝑘+1𝑒^𝑇_𝑘+1,𝐴^𝑇𝐴𝑉_𝑘=𝐴^𝑇𝑈_𝑘+1𝐵_𝑘=𝑉_𝑘+1𝐵^𝑇_𝑘𝛼_𝑘+1𝑒^𝑇_𝑘+1𝐵_𝑘=𝑉_𝑘+1𝐵^𝑇_𝑘𝐵_𝑘𝛼_𝑘+1𝛽_𝑘+1𝑒^𝑇_𝑘+1. (15)

While LSQR computes at each iteration the vector 𝑦_𝑘 which minimizes ‖𝑟_𝑘‖, LSMR chooses 𝑦_𝑘 differently, trying to minimize ‖𝐴^𝑇𝑟_𝑘‖. Since 𝑟_𝑘=𝑏−𝐴𝑉_𝑘𝑦_𝑘=𝛽₁𝑢₁−𝑈_𝑘+1𝐵_𝑘𝑦_𝑘=𝑈_𝑘+1𝛽₁𝑒₁−𝐵_𝑘𝑦_𝑘,𝐴^𝑇𝑟_𝑘=𝐴^𝑇𝑏−𝐴^𝑇𝐴𝑥_𝑘=𝛽₁𝛼₁𝑣₁−𝐴^𝑇𝐴𝑉_𝑘𝑦_𝑘=𝛽₁𝛼₁𝑣₁−𝑉_𝑘+1𝐵^𝑇_𝑘𝐵_𝑘𝛼_𝑘+1𝛽_𝑘+1𝑒^𝑇_𝑘+1𝑦_𝑘=𝑉_𝑘+1𝛽₁𝛼₁𝑒₁−𝐵^𝑇_𝑘𝐵_𝑘𝛼_𝑘+1𝛽_𝑘+1𝑒^𝑇_𝑘+1𝑦_𝑘. (16) So, we have the following two least squares problems: min_𝑦𝑘𝛽₁𝑒₁−𝐵_𝑘𝑦_𝑘,forLSQRmin_𝑦𝑘𝛽₁𝛼₁𝑒₁−𝐵^𝑇_𝑘𝐵_𝑘𝛼_𝑘+1𝛽_𝑘+1𝑒^𝑇_𝑘+1𝑦_𝑘,forLSMR. (17) Both ‖𝑟_𝑘‖ and ‖𝐴^𝑇𝑟_𝑘‖ decrease monotonically for LSMR algorithm, while only ‖𝑟_𝑘‖ is monotonic for LSQR algorithm. So, compared to LSQR, it is safer to terminate LSMR early.

LSMR works on the normal equation, so a simple analysis for the regularizing properties of the LSMR algorithm can be done, which is similar to the analysis of the CGLS method in [15]. Indeed LSMR and LSQR share the same Krylov subspace: 𝒦=span𝐴^𝑇𝑏,𝐴^𝑇𝐴𝐴^𝑇𝑏,𝐴^𝑇𝐴²𝐴^𝑇𝑏,…,𝐴^𝑇𝐴^𝑘−1𝐴^𝑇𝑏. (18) We insert the SVD of 𝐴 into the expansion 𝒦=span𝑉Σ𝑈^𝑇𝑏,𝑉Σ³𝑈^𝑇𝑏,…,𝑉Σ^2𝑘−1𝑈^𝑇𝑏. (19) The Krylov vectors for LSMR are 𝑉Σ^2𝑖−1𝑈^𝑇𝑏, 𝑖=1,2,…,𝑘. The diagonal elements of Σ and the coefficients in 𝑈^𝑇𝑏 decay, due to the discrete Picard condition. So, the Krylov subspace 𝒦 can be considered as an approximation to the subspace spanned by the first right singular vectors. Hence, 𝑥_𝑘 can be considered as an approximation to the TSVD solution.

From another perspective, LSMR is equivalent to the MINRES method applied to the normal equation. The regularizing properties of MINRES have been analysed in [14, 15]. So, LSMR has the regularizing properties apparently.

3. Choosing the Regularization Parameter

A reliable, efficient, and robust method for choosing the regularization parameter has not yet been found. Many techniques have been proposed to choose a suitable regularization parameter that works well under certain conditions. New methods are being developed too.

Several well-known methods, including the DP [19], the 𝐿-curve [20], the generalized cross validation (GCV) [21], and the NCP method, are described in Hansen’s book [6], and a comparison of these methods is involved. DP is often used to choose the regularization parameter of Krylov subspace methods. When ‖𝑥_𝑘‖ and ‖𝑟_𝑘‖ show monotonic behavior, the discrete 𝐿-curve criterion can be utilized. Hnětynková et al. recently proposed to use the information from the Golub-Kahan iterative bidiagonalization to estimate the norm of the noise, which can be used to construct efficient stopping criteria [22]. For more information about choosing the regularization parameter, we suggest [23] and references therein.

Two frequently used methods are utilized to choose the regularization parameter in this paper: the DP and the discrete 𝐿-curve criterion.

3.1. The Discrepancy Principle

One simple and frequently used method for choosing the regularization parameter is the DP ([14, Chapter 3.3], [6, Chapter 5.2]). It suggests that we should choose parameter 𝑘 such that the residual norm equals the discrepancy, which is the norm of the noise 𝑒. Generally, we terminate the iteration as soon as we have a solution 𝑥_𝑘, such that 𝑏−𝐴𝑥_𝑘≤𝜂,𝑏−𝐴𝑥_𝑘−1>𝜂, (20) where 𝜂 is the norm of 𝑒.

The main disadvantage of the DP method is that 𝜂 is often not available in advance. Because the choice of regularization parameter is sensitive to the accuracy of the estimate for 𝜂, it is not easy to get a good parameter 𝑘 if we use an estimate of the 𝜂.

3.2. The Discrete 𝐿-Curve Criterion

The discrete 𝐿-curve criterion is more complicated than DP. The 𝐿-curve tries to get the regularization parameter by the analysis of the norm of regularized solution ‖𝑥_𝑘‖ and corresponding residual norm ‖𝑏−𝐴𝑥_𝑘‖. The 𝐿-curve is a log-log plot of ‖𝑥_𝑘‖ versus ‖𝑏−𝐴𝑥_𝑘‖, and it has two different parts, an approximately horizontal part and an approximately vertical part.(1)When 𝑘 is small, 𝑥_𝑘≈𝑥_trueapproximatesaconstant,𝑏−𝐴𝑥_𝑘decreasesas𝑘increases. (21)(2)When 𝑘 is large enough, the solution is dominated by the error caused by noise, 𝑥_𝑘increasesas𝑘increases,𝑏−𝐴𝑥_𝑘≈‖𝑒‖approximatesaconstant. (22)

This also gives an explanation why DP is effective to choose the regularization parameter.

The effective regularization parameter can be determined by the “corner" that separates the two parts. The corner is the point with maximum curvature for continuous 𝐿-curves. It is similar to determine the corner for discrete 𝐿-curves. However, it is not straightforward to find the corner, because the discrete 𝐿-curve may have many local corners. In this paper, we adopt the adaptive pruning algorithm in [24] for the discrete 𝐿-curve criterion, which uses a sequence of pruned 𝐿-curves at different scales to get the global corner. The MATLAB implementations of this algorithm are available at the website (http://www.mathworks.com/matlabcentral/fileexchange/52-regtools/content/regu/corner.m).

The discrete 𝐿-curve criterion does not need information of the norm of noise, which is an advantage in comparison with the DP. The 𝐿-curve also has its limitations [25–27]: it will fail when the solution is very smooth; the regularization parameter computed by the 𝐿-curve may not behave consistently with the problem size increasing.

The discrete 𝐿-curve criterion has a requirement that the norm of 𝑥_𝑘 increases monotonically as 𝑘 increases and the norm of 𝑟_𝑘 decreases monotonically as 𝑘 increases. CGLS and LSQR are mathematically equivalent, and both of them can use the 𝐿-curve criterion as a stopping rule. The monotonicities of ‖𝑟_𝑘‖ in CGLS and LSQR are apparent. Starting with the zero vector, the monotonicities of ‖𝑥_𝑘‖ in CGLS and LSQR have been proved in [3, Theorem 6.3.1] and [28, Theorem 3.3.1]. The monotonicities of ‖𝑟_𝑘‖ and ‖𝑥_𝑘‖ in LSMR have been proved in [28, Theorem 3.3.11 and Theorem 3.3.6] accordingly. So, 𝐿-curve is allowed to choose the regularization parameter when LSMR is used as an iterative regularization method.

In practical implementations, we need stopping criterions to terminate the iteration (after the corner appears). Then, we can get the regularization parameter using the information we have. The iteration for choosing the regularization parameter needs effective stopping criterion. The frequently-used stopping criterion ‖𝑟_𝑘‖<tol (0<tol<1) is not a good one. The iteration may be terminated before the “corner” of 𝐿-curve appears for a bigger tol, while the number of iterations may be very large for a smaller tol because of the ill-condition of 𝐴. We choose 𝑥_𝑘>2‖𝑏‖,𝑥_𝑘−1≤2‖𝑏‖ (23) as a stopping criterion thanks to the monotonicity of ‖𝑥_𝑘‖. This stopping criterion is used with other stopping criterions together: 𝑟_𝑘<tol(0<tol<1),forasmalltol, (24) and 𝑘≥Maxitn, where Maxitn is an explicit limit on the iteration number.

4. Numerical Experiments

In this section we present some numerical examples to compare the IDR(𝑠) and LSMR methods with the LSQR method. All examples are problems of the restoration of blurred and noisy images. The blurring operator 𝐴 is constructed by 𝑏𝑙𝑢𝑟.𝑚 in regularization tools package by Hansen [29], which is a block Toeplitz with Toeplitz block matrix that models the blurring of an image by a Gaussian PSF. This section includes two subsections: the first part is the numerical examples for comparison of the IDR(𝑠) method and LSQR method where the DP is used to terminate the iteration; the second part is the numerical examples for comparison of the LSMR method and the LSQR method where the regularization parameters are chosen by the discrete 𝐿-curve criterion. The relative error norm (REN) is defined by 𝑥REN=^true−𝑥_𝑘2𝑥_true2. (25)

All experiments presented in this section are implemented on a computer with an Intel(R) Pentium(R) D CPU (3.00 GHz) and 1 GB of RAM using MATLAB R2010a.

4.1. The Comparison of IDR(𝑠) and LSQR

In this section, the IDR(𝑠)-biortho method is used due to its stability and accuracy.

Example 1. The example is the restoration of a noisy and blurred image “Lena.” The true image is shown in Figure 1(a). The image has 128×128 pixels. The parameters of 𝑏𝑙𝑢𝑟.𝑚 are set as the bandwidth band=10, 𝜎=1. The parameter 𝜎 controls the width of the Gaussian PSF, and the larger the 𝜎 is, the more ill-conditioned the problem is. Different Gaussian white noises with noise levels 𝛿=𝜂/‖𝑏‖=10⁻²,10⁻³,10⁻⁴ are added to the noise-free image in our experiments.

(a)

(b)

(c)

(d)

Figure 1 

(a) True image, (b) blurred and noisy image with 𝜎=1, (c) restored image with LSQR, and (d) restored image with IDR(1).

The IDR(𝑠) method is applied to solve this problem, while the popular iterative regularization method LSQR is used for comparison. The parameter 𝑠 in IDR(𝑠) is chosen as 𝑠=1, 2, 4, and 8. The DP method is used to terminate the iteration.

The cost for each IDR(𝑠) step is dominated by one matrix-vector product, while LSQR needs two matrix-vector products. Actually, IDS(𝑠) need more memory and computation as 𝑠 increases. More information about the cost of IDR(𝑠) can be found in [12].

The data in Table 1 indicates that IDR(𝑠) can produce a satisfying solution with much less iterations than LSQR (“Iter” is the number of iterations and “MV” is the number of matrix-vector products). Moreover, the parameter 𝑠 has a small influence on the convergence and the number of iterations. This suggests that we should use IDR(1) to solve the problem rather than IDR(𝑠) (𝑠>1) for saving storage and computational load.

Table 1 

Image deblurring using LSQR and IDR(𝑠) for “Lena” (𝜎=1).

However, IDR(𝑠) does not give a desirable solution when the noise level is high. Observing the relative error norm when 𝛿=10⁻², IDR(𝑠) cannot give a good enough solution compared to LSQR. The performance of IDR(s) is worse if the noise level is higher. Another shortcoming of IDR(𝑠) is that it cannot obtain a desirable solution or even cannot converge for very ill-condition problems which can be obtained with a larger 𝜎. The data in Table 2 shows the performance of LSQR and IDR(1) for Example 1 with 𝜎=1.5 (the condition number of 𝐴 is about 10⁹). Although IDR(𝑠) can obtain the solution much faster than LSQR, we should use these methods carefully. In our experiments, IDR(𝑠) can always give a satisfactory solution when 𝐴 is not very ill-conditioned and the noise level is low.

Table 2 

Image deblurring using LSQR and IDR(𝑠) for “Lena” (𝜎=1.5).

Figure 1 shows the true image and blurred and noisy image with 𝛿=10⁻³, 𝜎=1 and corresponding images restored by LSQR and IDR(1).

4.2. The Comparison of LSMR and LSQR

Example 2. The true image is constructed by MATLAB command phantom(𝑚) which can create a head phantom image with 𝑚×𝑚 pixels. In this experiment, we set the size 𝑚=64 and the bandwidth band=5, 𝜎=2. 𝐿-curve is fit for choosing the regularization parameter in this situation. The Gaussian white noise is added to the true image with noise level 𝛿=𝜂/‖𝑏‖=10⁻². The true image is Figure 3(a), and the blurred and noisy image is Figure 3(b).

Example 3. The true image is obtained from regularization tools package by Hansen [29]. The remaining parameters are the same as the parameters of Example 2. The true image and the blurred and noisy image are shown in Figures 5(a) and 5(b).

The LSMR method is applied to solve this problem, while the LSQR method is used for comparison. The discrete 𝐿-curve criterion is employed to choose the regularization parameter. We generate 100 instances of Gaussian white noise with noise level 𝛿=10⁻². For each instance, we use the discrete 𝐿-curve criterion to get the regularization parameters of LSMR and LSQR. Then, we obtain the solutions by running corresponding LSMR algorithm and LSQR algorithm. The comparison of the relative error norm is shown in Figure 2 (Example 2) and Figure 4 (Example 3). It is evident that LSMR more likely produces a better solution than LSQR. Figure 3 shows one instance of Example 2 including the true image, blurred and noisy image, and images restored by LSQR and LSMR. Similarly, Figure 5 shows one instance of Example 3. Furthermore, if we compare the optimal solution produced by LSMR and LSQR, the corresponding REN of the solution produced by LSMR is usually the smaller one.

Figure 2 

Comparison of the REN produced by LSQR and LSMR.

(a)

(b)

(c)

(d)

Figure 3 

(a) True image, (b) blurred and noisy image, (c) restored image with LSQR, Iter = 104, and REN = 0.5004, and (d) restored image with LSMR, Iter = 128, and REN = 0.4919.

Figure 4 

Comparison of the REN produced by LSQR and LSMR.

(a)

(b)

(c)

(d)

Figure 5 

(a) True image, (b) blurred and noisy image, (c) restored image with LSQR, Iter = 95, and REN = 0.3519, and (d) restored image with LSMR, Iter = 129, and REN = 0.3417.

The 𝐿-curve works well for large noise level. Actually, we found that the regularization parameter obtained from the 𝐿-curve is almost optimal when the level of noise is 5×10⁻². But the 𝐿-curve may fail or need a large number of iterations to find the corner for small noise level.

Besides, we have done the comparison of LSMR and LSQR with the regularization parameter chosen by DP. LSMR method can get nearly the same REN as LSQR method, while LSQR method could be terminated somewhat sooner. The memory and computational costs of every step in LSMR are also slightly more expensive than that in LSQR [11].

5. Conclusion

In this paper, we briefly introduced the IDR(𝑠) method and the LSMR method, and then we applied the IDR(𝑠) and LSMR methods to solve image deblurring problems. These methods show some superiorities when compared to the classic iterative regularization method LSQR.

When DP is used to terminate the iteration, IDR(𝑠) can give a satisfactory solution with a low computational cost when 𝐴 is not very ill-condition and the noise level is low. The parameter 𝑠 has little influence on the performance of IDR(𝑠), so we can use IDR(1) for saving storage and computational cost.

LSMR performs as well as LSQR but needs more computational cost when DP is used to terminate the iteration. When the 𝐿-curve method is utilized to choose the regularization parameter, LSMR more likely produces a more accurate solution than LSQR.

Acknowledgments

This research is supported by 973 Program (2013CB329404), NSFC (61170311), Chinese Universities Specialized Research Fund for the Doctoral Program (20110185110020), and Sichuan Province Sci. & Tech. Research Project (2012GZX0080). The authors would like to express their thankfulness to the referees and the editor for their much helpful suggestions for revising this paper.