# Two New Efficient Iterative Regularization Methods for Image Restoration Problems

Authors: Chao Zhao; Ting-Zhu Huang; Xi-Le Zhao; Liang-Jian Deng
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

Iterative regularization methods are efficient regularization tools for image restoration problems. The IDR(s) 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(s) and LSMR methods for image restoration problems. Comparative numerical experiments show that IDR(s) 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 L-curve method to choose the regularization parameter.


---

## Body

## 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:(1)∫ℝ2‍K(s-s′,t-t′)f(s′,t′)ds′dt′=g(s,t),K(s-s′,t-t′)f(s′,t′)ds′ddt′=(s,t)∈ℝ2,
where f and g are the original image and the blurred image, respectively. K 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
(2)Axtrue=btrue,A∈ℝm2×m2,x,b∈ℝm2,
where the matrix A is ill-conditioned since it has many singular values close to zero [1]. For simplicity, we assume that A in this paper is nonsingular.However, the right-hand sideb 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
(3)Ax=b,b=btrue+e,e∈ℝm2,
where e represents the noise.Our goal is to obtain a good approximation of the original imagextrue by solving the system (3) instead of the system (2) since btrue is not known. However, the solution of (3) is not a good approximation of the solution of (2) because of ill-conditioned A. In fact, there is a quite remarkable disparity among the corresponding solutions of (3) and (2) even if the norm of e is small.The singular value decomposition (SVD) is a powerful tool for the analysis of discrete ill-posed problems. The SVD takes the form(4)A=UΣVT=∑i=1m2‍uiσiviT,
where Σ is a diagonal matrix with the singular values σ1≥σ2≥⋯≥σm2≥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 |uiTb| decay faster than the σi.Regularization methods are used to produce better solution. One of the most famous regularization methods is truncated SVD (TSVD) [2], which yields solution xk=∑i=1k‍σi-1(uiTb)vi. Another regularization method is Tikhonov regularization which consists of solving the minimization problem
(5)minx∈Rm2{∥Ax-b∥22+λ∥Lx∥22},
where L 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 {x1,x2,…,xk,…} 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 k can be considered as a regularization parameter. We also introduce different methods for choosing an effective k 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 A is not SPD, the method can be applied to the normal equations ATAx=ATb. 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(s) [10] and LSMR [11] iterative methods have been proposed for solving large linear system. IDR(s) is a new family of iterative methods which is fast and requires low memory. IDR(s) 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 ∥Ark∥ and ∥rk∥ decrease monotonically for LSMR, while only ∥rk∥ 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 discreteL-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 L-curve criterion is used to choose the regularization parameter by using the curve of (log∥rk∥,log∥xk∥). For more information about methods for choosing the regularization parameter, we suggest [6, Chapter  5] and references therein.Inspired by the advantages of IDR(s) and LSMR, we apply them to solve the image restoration problems. We compare the performance of the IDR(s) 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 L-curve criterion is used to choose the regularization parameter.The paper is organized as follows. In Section2, we review the IDR(s) and LSMR algorithms briefly and present analysis of the regularizing properties of the IDR(s) and LSMR algorithms. In Section 3, methods for choosing the regularization parameter especially the DP method and the discrete L-curve criterion are introduced. IDR(s) and LSMR methods are applied to solve image restoration problems and compared with with LSQR in Section 4. Section 5 summarizes conclusions.


## 2. IDR(s) and LSMR

### 2.1. The IDR(s) Method

The IDR(s) methods are a new family of iterative methods based on the induced dimension reduction (IDR) method proposed by Sonneveld in 1980. IDR(s) 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(6)𝒢j=(I-ωjA)(𝒢j-1∩𝒮),
where 𝒮 is a fixed proper subspace of ℝN and ωj are nonzero scalars. The space 𝒮 is usually assumed as the left null space of P∈ℝN×s:
(7)P=(p1p2…ps),𝒮=𝒩(PH).The columns ofP are chosen as orthogonalization of a set of random vectors.Suppose that aftern iterations we have s vectors gi∈𝒢j, i=n-1,…,n-s and s corresponding vectors ui with gi=Aui. Define the matrices
(8)Gn=(gn-sgn-s+1⋯gn-1),Un=(un-sun-s+1⋯un-1).The algorithm producess+1 residuals in 𝒢j; then the next residual can be found in 𝒢j+1. Two steps are needed to obtain the s+1 residuals.Step  1
Compute the first residual in𝒢j+1.
Calculatec from (PHGn)c=PHrn
v
n
=
r
n
-
G
n
c
r
n
+
1
=
r
n
-
ω
j
+
1
v
n
-
G
n
c
x
n
+
1
=
x
n
+
ω
j
+
1
A
v
n
+
U
n
c.Step  2
Compute the remaining residuals in𝒢j+1, for k=1,…s.
Calculatec from (PHGn+k-1)c=PHrn+k
v
n
+
k
=
r
n
+
k
-
G
n
+
k
-
1
c
u
n
+
k
=
ω
j
+
1
v
n
+
k
+
U
n
+
k
-
1
c
g
n
+
k
=
A
u
n
+
k
selectαi and βi
g
n
+
k
=
g
n
+
k
-
∑
i
=
1
k
-
1
‍
α
i
g
n
+
i; un+k=un+k-∑i=1k-1‍αiun+i;
rn+k+1=rn+k-∑i=1k‍βign+i; xn+k+1=xn+k+∑i=1k‍βiun+i;
updateG:,k=gn+k to get Gn+k+1; U:,k=un+k to get Un+k+1.Different choices ofαi and βi would result in different variants of IDR(s). Take the same notations as in [12]: IDR(s)-proto represents the algorithm in [10] and IDR(s)-biortho represents the variant of the IDR(s) in [12]. IDR(s)-proto just computes
(9)rn+k+1=rn+k-gn+k,xn+k+1=xn+k+uu+k,
that is, αi=0(i=1,…,k-1), βi=0(i=1,…,k-1), and βi=1(i=k).IDR(s)-biortho selects αi and βi to construct vectors that satisfy the biorthogonality conditions:
(10)gn+k⊥pi,i=1,…,k-1,k=2,…,s,rn+k+1⊥pi,i=1,…,k,k=1,…,s.
The IDR(s)-biortho algorithm is more stable and more accurate than the IDR(s)-proto algorithm, so we apply the IDR(s)-biortho algorithm to solve the problems in this paper.As described in [10, 16], we have the following expression after n iterations:
(11)rn=Φn(A)r0,Φn(A) is a polynomial of degree n. Φn(A) can be explicitly written as a product of two polynomials
(12)Φn(A)=Ωi(A)Ψn-i(A),
where Ωi(t)=∏k=1i‍(1-ωkt), Ψn-i(t)=∏j=0n-i‍cjtj. The factor Ωi(A) is called damping factor and Ψn-i(A) 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 P, whereas the damping factor depends on the choice of the parameters ωi. Actually, in our experiments, different choices of parameter s (which controls the dimension of P) 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:(13)β1u1=b,α1v1=ATu1fork=1,2,…,βk+1uk+1=Avk-αkuk,αk+1vk+1=ATuk-βk+1vk.
The scalars αk≥0 and βk≥0 are chosen so that the norms of uk and vk are both equal to one.We define thatVk=(v1v2⋯vk), Uk=(u1u2⋯uk), and
(14)Bk=(α1β2α2⋱⋱βkαkβk+1).
Then, the following relations hold:
(15)Uk+1(β1e1)=b,AVk=Uk+1Bk,ATUk+1=Vk+1(BkTαk+1ek+1T),ATAVk=ATUk+1Bk=Vk+1(BkTαk+1ek+1T)Bk=Vk+1(BkTBkαk+1βk+1ek+1T).While LSQR computes at each iteration the vectoryk which minimizes ∥rk∥, LSMR chooses yk differently, trying to minimize ∥ATrk∥. Since
(16)rk=b-AVkyk=β1u1-Uk+1Bkyk=Uk+1(β1e1-Bkyk),ATrk=ATb-ATAxk=β1α1v1-ATAVkyk=β1α1v1-Vk+1(BkTBkαk+1βk+1ek+1T)yk=Vk+1(β1α1e1-(BkTBkαk+1βk+1ek+1T)yk).
So, we have the following two least squares problems:
(17)minyk∥β1e1-Bkyk∥,forLSQRminyk∥(β1α1e1-(BkTBkαk+1βk+1ek+1T)yk)∥,forLSMR.
Both ∥rk∥ and ∥ATrk∥ decrease monotonically for LSMR algorithm, while only ∥rk∥ 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:
(18)𝒦=span{ATb,(ATA)ATb,(ATA)2ATb,…,(ATA)k-1ATb}.
We insert the SVD of A into the expansion
(19)𝒦=span{VΣUTb,VΣ3UTb,…,VΣ2k-1UTb}.
The Krylov vectors for LSMR are VΣ2i-1UTb, i=1,2,…,k. The diagonal elements of Σ and the coefficients in UTb 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, xk 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.


## 2.1. The IDR(s) Method

The IDR(s) methods are a new family of iterative methods based on the induced dimension reduction (IDR) method proposed by Sonneveld in 1980. IDR(s) 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(6)𝒢j=(I-ωjA)(𝒢j-1∩𝒮),
where 𝒮 is a fixed proper subspace of ℝN and ωj are nonzero scalars. The space 𝒮 is usually assumed as the left null space of P∈ℝN×s:
(7)P=(p1p2…ps),𝒮=𝒩(PH).The columns ofP are chosen as orthogonalization of a set of random vectors.Suppose that aftern iterations we have s vectors gi∈𝒢j, i=n-1,…,n-s and s corresponding vectors ui with gi=Aui. Define the matrices
(8)Gn=(gn-sgn-s+1⋯gn-1),Un=(un-sun-s+1⋯un-1).The algorithm producess+1 residuals in 𝒢j; then the next residual can be found in 𝒢j+1. Two steps are needed to obtain the s+1 residuals.Step  1
Compute the first residual in𝒢j+1.
Calculatec from (PHGn)c=PHrn
v
n
=
r
n
-
G
n
c
r
n
+
1
=
r
n
-
ω
j
+
1
v
n
-
G
n
c
x
n
+
1
=
x
n
+
ω
j
+
1
A
v
n
+
U
n
c.Step  2
Compute the remaining residuals in𝒢j+1, for k=1,…s.
Calculatec from (PHGn+k-1)c=PHrn+k
v
n
+
k
=
r
n
+
k
-
G
n
+
k
-
1
c
u
n
+
k
=
ω
j
+
1
v
n
+
k
+
U
n
+
k
-
1
c
g
n
+
k
=
A
u
n
+
k
selectαi and βi
g
n
+
k
=
g
n
+
k
-
∑
i
=
1
k
-
1
‍
α
i
g
n
+
i; un+k=un+k-∑i=1k-1‍αiun+i;
rn+k+1=rn+k-∑i=1k‍βign+i; xn+k+1=xn+k+∑i=1k‍βiun+i;
updateG:,k=gn+k to get Gn+k+1; U:,k=un+k to get Un+k+1.Different choices ofαi and βi would result in different variants of IDR(s). Take the same notations as in [12]: IDR(s)-proto represents the algorithm in [10] and IDR(s)-biortho represents the variant of the IDR(s) in [12]. IDR(s)-proto just computes
(9)rn+k+1=rn+k-gn+k,xn+k+1=xn+k+uu+k,
that is, αi=0(i=1,…,k-1), βi=0(i=1,…,k-1), and βi=1(i=k).IDR(s)-biortho selects αi and βi to construct vectors that satisfy the biorthogonality conditions:
(10)gn+k⊥pi,i=1,…,k-1,k=2,…,s,rn+k+1⊥pi,i=1,…,k,k=1,…,s.
The IDR(s)-biortho algorithm is more stable and more accurate than the IDR(s)-proto algorithm, so we apply the IDR(s)-biortho algorithm to solve the problems in this paper.As described in [10, 16], we have the following expression after n iterations:
(11)rn=Φn(A)r0,Φn(A) is a polynomial of degree n. Φn(A) can be explicitly written as a product of two polynomials
(12)Φn(A)=Ωi(A)Ψn-i(A),
where Ωi(t)=∏k=1i‍(1-ωkt), Ψn-i(t)=∏j=0n-i‍cjtj. The factor Ωi(A) is called damping factor and Ψn-i(A) 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 P, whereas the damping factor depends on the choice of the parameters ωi. Actually, in our experiments, different choices of parameter s (which controls the dimension of P) 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:(13)β1u1=b,α1v1=ATu1fork=1,2,…,βk+1uk+1=Avk-αkuk,αk+1vk+1=ATuk-βk+1vk.
The scalars αk≥0 and βk≥0 are chosen so that the norms of uk and vk are both equal to one.We define thatVk=(v1v2⋯vk), Uk=(u1u2⋯uk), and
(14)Bk=(α1β2α2⋱⋱βkαkβk+1).
Then, the following relations hold:
(15)Uk+1(β1e1)=b,AVk=Uk+1Bk,ATUk+1=Vk+1(BkTαk+1ek+1T),ATAVk=ATUk+1Bk=Vk+1(BkTαk+1ek+1T)Bk=Vk+1(BkTBkαk+1βk+1ek+1T).While LSQR computes at each iteration the vectoryk which minimizes ∥rk∥, LSMR chooses yk differently, trying to minimize ∥ATrk∥. Since
(16)rk=b-AVkyk=β1u1-Uk+1Bkyk=Uk+1(β1e1-Bkyk),ATrk=ATb-ATAxk=β1α1v1-ATAVkyk=β1α1v1-Vk+1(BkTBkαk+1βk+1ek+1T)yk=Vk+1(β1α1e1-(BkTBkαk+1βk+1ek+1T)yk).
So, we have the following two least squares problems:
(17)minyk∥β1e1-Bkyk∥,forLSQRminyk∥(β1α1e1-(BkTBkαk+1βk+1ek+1T)yk)∥,forLSMR.
Both ∥rk∥ and ∥ATrk∥ decrease monotonically for LSMR algorithm, while only ∥rk∥ 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:
(18)𝒦=span{ATb,(ATA)ATb,(ATA)2ATb,…,(ATA)k-1ATb}.
We insert the SVD of A into the expansion
(19)𝒦=span{VΣUTb,VΣ3UTb,…,VΣ2k-1UTb}.
The Krylov vectors for LSMR are VΣ2i-1UTb, i=1,2,…,k. The diagonal elements of Σ and the coefficients in UTb 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, xk 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 L-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 ∥xk∥ and ∥rk∥ show monotonic behavior, the discrete L-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 discreteL-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 k such that the residual norm equals the discrepancy, which is the norm of the noise e. Generally, we terminate the iteration as soon as we have a solution xk, such that
(20)∥b-Axk∥≤η,∥b-Axk-1∥>η,
where η is the norm of e.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 k if we use an estimate of the η.


### 3.2. The DiscreteL-Curve Criterion

The discreteL-curve criterion is more complicated than DP. The L-curve tries to get the regularization parameter by the analysis of the norm of regularized solution ∥xk∥ and corresponding residual norm ∥b-Axk∥. The L-curve is a log-log plot of ∥xk∥ versus ∥b-Axk∥, and it has two different parts, an approximately horizontal part and an approximately vertical part.(1)
Whenk is small,
(21)∥xk∥≈∥xtrue∥approximates a constant,∥b-Axk∥decreases askincreases.(2)
Whenk is large enough, the solution is dominated by the error caused by noise,
(22)∥xk∥increases askincreases,∥b-Axk∥≈∥e∥approximates a constant.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 continuousL-curves. It is similar to determine the corner for discrete L-curves. However, it is not straightforward to find the corner, because the discrete L-curve may have many local corners. In this paper, we adopt the adaptive pruning algorithm in [24] for the discrete L-curve criterion, which uses a sequence of pruned L-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 discreteL-curve criterion does not need information of the norm of noise, which is an advantage in comparison with the DP. The L-curve also has its limitations [25–27]: it will fail when the solution is very smooth; the regularization parameter computed by the L-curve may not behave consistently with the problem size increasing.The discreteL-curve criterion has a requirement that the norm of xk increases monotonically as k increases and the norm of rk decreases monotonically as k increases. CGLS and LSQR are mathematically equivalent, and both of them can use the L-curve criterion as a stopping rule. The monotonicities of ∥rk∥ in CGLS and LSQR are apparent. Starting with the zero vector, the monotonicities of ∥xk∥ in CGLS and LSQR have been proved in [3, Theorem 6.3.1] and [28, Theorem 3.3.1]. The monotonicities of ∥rk∥ and ∥xk∥ in LSMR have been proved in [28, Theorem 3.3.11 and Theorem 3.3.6] accordingly. So, L-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∥rk∥<tol (0<tol<1) is not a good one. The iteration may be terminated before the “corner” of L-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 A. We choose
(23)∥xk∥>2∥b∥,∥xk-1∥≤2∥b∥
as a stopping criterion thanks to the monotonicity of ∥xk∥. This stopping criterion is used with other stopping criterions together:
(24)∥rk∥<tol(0<tol<1),for a small tol,
and k≥Maxitn, where Maxitn is an explicit limit on the iteration number.


## 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 k such that the residual norm equals the discrepancy, which is the norm of the noise e. Generally, we terminate the iteration as soon as we have a solution xk, such that
(20)∥b-Axk∥≤η,∥b-Axk-1∥>η,
where η is the norm of e.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 k if we use an estimate of the η.


## 3.2. The DiscreteL-Curve Criterion

The discreteL-curve criterion is more complicated than DP. The L-curve tries to get the regularization parameter by the analysis of the norm of regularized solution ∥xk∥ and corresponding residual norm ∥b-Axk∥. The L-curve is a log-log plot of ∥xk∥ versus ∥b-Axk∥, and it has two different parts, an approximately horizontal part and an approximately vertical part.(1)
Whenk is small,
(21)∥xk∥≈∥xtrue∥approximates a constant,∥b-Axk∥decreases askincreases.(2)
Whenk is large enough, the solution is dominated by the error caused by noise,
(22)∥xk∥increases askincreases,∥b-Axk∥≈∥e∥approximates a constant.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 continuousL-curves. It is similar to determine the corner for discrete L-curves. However, it is not straightforward to find the corner, because the discrete L-curve may have many local corners. In this paper, we adopt the adaptive pruning algorithm in [24] for the discrete L-curve criterion, which uses a sequence of pruned L-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 discreteL-curve criterion does not need information of the norm of noise, which is an advantage in comparison with the DP. The L-curve also has its limitations [25–27]: it will fail when the solution is very smooth; the regularization parameter computed by the L-curve may not behave consistently with the problem size increasing.The discreteL-curve criterion has a requirement that the norm of xk increases monotonically as k increases and the norm of rk decreases monotonically as k increases. CGLS and LSQR are mathematically equivalent, and both of them can use the L-curve criterion as a stopping rule. The monotonicities of ∥rk∥ in CGLS and LSQR are apparent. Starting with the zero vector, the monotonicities of ∥xk∥ in CGLS and LSQR have been proved in [3, Theorem 6.3.1] and [28, Theorem 3.3.1]. The monotonicities of ∥rk∥ and ∥xk∥ in LSMR have been proved in [28, Theorem 3.3.11 and Theorem 3.3.6] accordingly. So, L-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∥rk∥<tol (0<tol<1) is not a good one. The iteration may be terminated before the “corner” of L-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 A. We choose
(23)∥xk∥>2∥b∥,∥xk-1∥≤2∥b∥
as a stopping criterion thanks to the monotonicity of ∥xk∥. This stopping criterion is used with other stopping criterions together:
(24)∥rk∥<tol(0<tol<1),for a small tol,
and k≥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(s) and LSMR methods with the LSQR method. All examples are problems of the restoration of blurred and noisy images. The blurring operator A is constructed by blur.m 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(s) 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 L-curve criterion. The relative error norm (REN) is defined by
(25)REN=∥xtrue-xk∥2∥xtrue∥2.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(s) and LSQR

In this section, the IDR(s)-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 Figure1(a). The image has 128×128 pixels. The parameters of blur.m 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 δ=η/∥b∥=10-2,10-3,10-4 are added to the noise-free image in our experiments.(a) True image, (b) blurred and noisy image withσ=1, (c) restored image with LSQR, and (d) restored image with IDR(1).
(a)
(b)
(c)
(d)The IDR(s) method is applied to solve this problem, while the popular iterative regularization method LSQR is used for comparison. The parameter s in IDR(s) is chosen as s=1, 2, 4, and 8. The DP method is used to terminate the iteration.The cost for each IDR(s) step is dominated by one matrix-vector product, while LSQR needs two matrix-vector products. Actually, IDS(s) need more memory and computation as s increases. More information about the cost of IDR(s) can be found in [12].The data in Table1 indicates that IDR(s) 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 s 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(s) (s>1) for saving storage and computational load.Table 1
Image deblurring using LSQR and IDR(s) for “Lena” (σ=1).
Method
LSQR
IDR(1)
IDR(2)
IDR(4)
IDR(8)
δ
=
1
0
-
2
Iter
8
4
6
5
9
MV
16
4
6
5
9
REN
0.0625
0.0705
0.1352
0.0799
0.1832
δ
=
1
0
-
3
Iter
32
10
10
9
13
MV
64
10
10
9
13
REN
0.0347
0.0323
0.0314
0.0313
0.0319
δ
=
1
0
-
4
Iter
145
32
33
30
35
MV
290
32
33
30
35
REN
0.0170
0.0135
0.0175
0.0140
0.0177However, IDR(s) does not give a desirable solution when the noise level is high. Observing the relative error norm when δ=10-2, IDR(s) 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(s) 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 A is about 109). Although IDR(s) can obtain the solution much faster than LSQR, we should use these methods carefully. In our experiments, IDR(s) can always give a satisfactory solution when A is not very ill-conditioned and the noise level is low.Table 2
Image deblurring using LSQR and IDR(s) for “Lena” (σ=1.5).
δ
=
1
0
-
2
δ
=
1
0
-
3
δ
=
1
0
-
4
Iter
MV
REN
Iter
MV
REN
Iter
MV
REN
LSQR
12
24
0.0919
58
116
0.0670
329
658
0.0486
IDR(s)
6
6
0.1344
36
36
0.7551
40
40
0.0488Figure1 shows the true image and blurred and noisy image with δ=10-3, σ=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 commandphantom(m) which can create a head phantom image with m×m pixels. In this experiment, we set the size m=64 and the bandwidth band=5, σ=2. L-curve is fit for choosing the regularization parameter in this situation. The Gaussian white noise is added to the true image with noise level δ=η/∥b∥=10-2. 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 discreteL-curve criterion is employed to choose the regularization parameter. We generate 100 instances of Gaussian white noise with noise level δ=10-2. For each instance, we use the discrete L-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) 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.
(a)
(b)
(c)
(d)Figure 4
Comparison of the REN produced by LSQR and LSMR.(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.
(a)
(b)
(c)
(d)TheL-curve works well for large noise level. Actually, we found that the regularization parameter obtained from the L-curve is almost optimal when the level of noise is 5×10-2. But the L-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].


## 4.1. The Comparison of IDR(s) and LSQR

In this section, the IDR(s)-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 Figure1(a). The image has 128×128 pixels. The parameters of blur.m 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 δ=η/∥b∥=10-2,10-3,10-4 are added to the noise-free image in our experiments.(a) True image, (b) blurred and noisy image withσ=1, (c) restored image with LSQR, and (d) restored image with IDR(1).
(a)
(b)
(c)
(d)The IDR(s) method is applied to solve this problem, while the popular iterative regularization method LSQR is used for comparison. The parameter s in IDR(s) is chosen as s=1, 2, 4, and 8. The DP method is used to terminate the iteration.The cost for each IDR(s) step is dominated by one matrix-vector product, while LSQR needs two matrix-vector products. Actually, IDS(s) need more memory and computation as s increases. More information about the cost of IDR(s) can be found in [12].The data in Table1 indicates that IDR(s) 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 s 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(s) (s>1) for saving storage and computational load.Table 1
Image deblurring using LSQR and IDR(s) for “Lena” (σ=1).
Method
LSQR
IDR(1)
IDR(2)
IDR(4)
IDR(8)
δ
=
1
0
-
2
Iter
8
4
6
5
9
MV
16
4
6
5
9
REN
0.0625
0.0705
0.1352
0.0799
0.1832
δ
=
1
0
-
3
Iter
32
10
10
9
13
MV
64
10
10
9
13
REN
0.0347
0.0323
0.0314
0.0313
0.0319
δ
=
1
0
-
4
Iter
145
32
33
30
35
MV
290
32
33
30
35
REN
0.0170
0.0135
0.0175
0.0140
0.0177However, IDR(s) does not give a desirable solution when the noise level is high. Observing the relative error norm when δ=10-2, IDR(s) 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(s) 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 A is about 109). Although IDR(s) can obtain the solution much faster than LSQR, we should use these methods carefully. In our experiments, IDR(s) can always give a satisfactory solution when A is not very ill-conditioned and the noise level is low.Table 2
Image deblurring using LSQR and IDR(s) for “Lena” (σ=1.5).
δ
=
1
0
-
2
δ
=
1
0
-
3
δ
=
1
0
-
4
Iter
MV
REN
Iter
MV
REN
Iter
MV
REN
LSQR
12
24
0.0919
58
116
0.0670
329
658
0.0486
IDR(s)
6
6
0.1344
36
36
0.7551
40
40
0.0488Figure1 shows the true image and blurred and noisy image with δ=10-3, σ=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 commandphantom(m) which can create a head phantom image with m×m pixels. In this experiment, we set the size m=64 and the bandwidth band=5, σ=2. L-curve is fit for choosing the regularization parameter in this situation. The Gaussian white noise is added to the true image with noise level δ=η/∥b∥=10-2. 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 discreteL-curve criterion is employed to choose the regularization parameter. We generate 100 instances of Gaussian white noise with noise level δ=10-2. For each instance, we use the discrete L-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) 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.
(a)
(b)
(c)
(d)Figure 4
Comparison of the REN produced by LSQR and LSMR.(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.
(a)
(b)
(c)
(d)TheL-curve works well for large noise level. Actually, we found that the regularization parameter obtained from the L-curve is almost optimal when the level of noise is 5×10-2. But the L-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(s) method and the LSMR method, and then we applied the IDR(s) 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(s) can give a satisfactory solution with a low computational cost when A is not very ill-condition and the noise level is low. The parameter s has little influence on the performance of IDR(s), 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 theL-curve method is utilized to choose the regularization parameter, LSMR more likely produces a more accurate solution than LSQR.



---

*Source: 129652-2013-06-26.xml*