# Global Quasi-Minimal Residual Method for Image Restoration

Authors: Jun Liu; Ting-Zhu Huang; Xiao-Guang Lv; Hao Xu; Xi-Le Zhao
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The global quasi-minimal residual (QMR) method is a popular iterative method for the solution of
linear systems with multiple right-hand sides. In this paper, we consider the application of the global
QMR method to classical ill-posed problems arising from image restoration. Since the scale of the problem
is usually very large, the computations with the blurring matrix can be very expensive. In this regard,
we use a Kronecker product approximation of the blurring matrix to benefit the computation. In order to
reduce the disturbance of noise to the solution, the Tikhonov regularization technique is adopted to produce
better approximation of the desired solution. Numerical results show that the global QMR method
outperforms the classic CGLS method and the global GMRES method.


---

## Body

## 1. Introduction

In the area of remote sensing, materials science, medical and astronomical imaging, and so on, image restoration plays an important role in preprocessing and postprocessing the image [1]. Many image restoration tasks can be posed as problems of the form (1)∬Ωhx,y;s,tfs,tdsdt=gx,y,where the functions f, g represent the original and blurred images, respectively. The kernel h is apoint spread function (PSF) which is a function that specifies the degree of blurring. PSFs are often classified as either spatially variant or spatially invariant [2, 3]. For simplicity, we take into account spatially variant PSF in this paper. By means of discretization methods such as the Galerkin method or quadrature method [4], (1) can be discretized to the following linear equations: (2)Ax=b~,A∈Rn2×n2,x,b~∈Rn2,where x is a vector representing the true image and b~ is a vector representing the blurred image, which are the discretized versions of f and g in (1), respectively. The matrix A is the blurring matrix constructed from the discretized version of the PSF h. It should be noted that the PSF is assumed to be known here. In fact, if the PSF is unknown, there are a variety of means of techniques available for estimating it [5, 6]. In real applications, the right-hand side error-free vector b~ is not accessible. Instead, the vector (3)b=b~+ηis known, where the vector η represents the additive noise. That is, the observed image is not only blurred but also contaminated with noise. Commonly, η is assumed to be the white Gaussian noise, and its Euclidean vector norm is considered to be a priori but the noise vector itself is not.In this work, we aim to obtain an approximation of the original imagex by computing a solution of the linear system of equations (4)Ax=b,A∈Rn2×n2,x,b∈Rn2.If the observed image array has dimension n×n, then f and g are vectors of length n2, and A is an n2×n2 matrix. Typical values of n are 256, 512, and 1024, so the dimensions of the matrix A can be extremely large [7]. Then the computations with A can be very expensive.Fortunately, the matrixA has a special structure when an appropriate boundary condition is imposed. Then the computational cost of matrix-vector multiplication can be alleviated to some extent. For large-scale problems, such as image restoration problems, the direct regularization method cannot always obtain good solutions, but the iterative method is a better choice. Krylov subspace iterative methods are the most commonly used approaches that can be employed for solving (4). In [8], the authors proposed to employ the well known BiCG and QMR methods for image restoration. They also considered using a popular iterative method GMRES which was first proposed by Saad and Schultz for image restoration in [9].Equation (4) can be replaced by new ones involving matrix equations, if the matrix A can be decomposed as Kronecker products, and then the computations with A can be reduced. In [10], the authors first proposed the global Krylov subspace methods to solve the matrix equations. The methods were proved to be very effective for large-scale matrix equations. Later in [11], Bouhamidi and Jbilou applied the global GMRES method to image restoration problems. Their numerical tests demonstrated that the global GMRES method was better than the GMRES method.Due to the error in the right-hand side and the severe ill-conditioning property of the matrixA, the straightforward solution of (4) typically does not yield a meaningful approximation [4, 12, 13]. Therefore, instead of solving the system (4) directly, we replace it by a nearby linear system with a less ill-conditioned matrix and solve the corresponding new linear system. This replacement is commonly referred to as regularization [9]. Probably the most renowned regularization approach to overcome ill-conditioning dates back to Tikhonov and Arsenin [14].In this paper, we consider the implementation of the global quasi-minimal residual (QMR) method for image restoration problems. The approach discussed here can be considered as an extension and a specific real application of the method introduced in [15] where the authors applied this method to solve the general Sylvester equation. This approach is motivated by the work of Bouhamidi and Jbilou in [11]. The numerical experiments show that the global QMR method is very effective compared with the global-GMRES method and the classic conjugate gradient method for least square (CGLS) problem.The outline of this paper is as follows. In the next section, we give some notations and definitions that will be used throughout this paper. Section3 introduces the global QMR method for image restoration problems. We present some numerical experiments to show the efficiency of the global QMR method in Section 4. Finally concluding remarks can be found in Section 5.


## 2. Preliminaries

As shown in [7], some blurring operators (e.g., Gaussian) are separable and therefore can be factored as Kronecker product of two matrices. The computation for the solution of (4) can be reduced if the Kronecker product approximation of A is employed.Suppose thatP∈Rn×n is the discretized PSF. If the PSF is separable, that is, P can be decomposed as (5)P=abT,where a and b are n×1 vectors, the matrix A constructed from P has block structure of the form (6)A=Ar⊗Ac=a11rAca12rAc⋯a1nrAca21rAca22rAc⋯a2nrAc⋮⋮⋮an1rAcan2rAc⋯annrAc,where the matrices Ar and Ac have parameters a and b, respectively, with the specific structures depending on the imposed boundary condition [6, 16], and “⊗” denotes Kronecker product. We refer the readers to [17] for details about the properties of Kronecker product.If the PSFP is inseparable, then the corresponding matrix A is inseparable. However, we can find the Kronecker product approximation of A by using SVD technique so that A can be approximately decomposed as the following form: (7)A≈∑i=1r‍Aai⊗Abi,where P=∑i=1r‍biaiT with a given integer r≤rank(P). In particular, the authors in [2, 16] pointed out that Aa1⊗Ab1 is the best (as measured by the Frobenius norm) Kronecker approximation of A.According to the properties of Kronecker product, (4) can be rewritten as (8)AcXArT=B,where x=vec(X) and b=vec(B). Note that z=vec(Z) with Z∈Rn×n is the n2×1 vector obtained by stacking n columns of the matrix Z. Define an operator A:X∈Rn×n→AcXArT and AT:X∈Rn×n→ArXAcT; then (4) can be rewritten as (9)AX=B.We use the notation (10)KmA,V=spanV,AV,A2V,…,Am-1V,for the global Krylov subspace of Rn×n generated by the matrix V∈Rn×n and the operator A. Note that (11)Z∈KmA,V⟺Z=∑i=1m‍αiAi-1V,αi∈R.LetX,Y∈Rn×n; we define the inner matrix product 〈X,Y〉F=tr(XTY), where tr(Z) denotes the trace of the square matrix Z and XT the transpose of the matrix X. The associated norm is the Frobenius norm ·F. The matrices X,Y are said to be F-orthonormal if tr(XTY)=0.In the following, we will introduce an algorithm of the global Lanczos biorthogonal process, which has been elaborately discussed in [15, 18]. This process is used to construct a pair of biorthogonal basis V1,V2,…,Vm and W1,W2,…,Wm of the two Krylov subspaces Km(A,V1) and Km(AT,W1), respectively, such that(12)Vi,Wj=trViTWj=δij=1,i=j,0,i≠j.fori,j=1,2,…,m.The construction process can be summarized as in Algorithm1.Algorithm 1:The global Lanczos biorthogonal process.
(1) GivenA,V1 and W1 such that 〈V1,W1〉F=1;
(2) Setβ1=δ1=0 and V0=W0=On×n;
(3) Fori=1,2,…,m
αj=tr⁡(WjTAVj);
V^j+1=AVj-αjVj-βjVj-1;
W^j+1=ATWj-αjWj-δjWj-1;
δj+1=|tr⁡(V^j+1TW^j+1)|1/2, if δj+1=0, stop;
βj+1=tr⁡(V^j+1TW^j+1)/δj+1;
Vj+1=V^j+1/δj+1;
Wj+1=W^j+1/βj+1;
EndFor convenience, we denote byVm and Wm the n×nm block matrix, that is, Vm=[V1,V2,…,Vm] and Wm=[W1,W2,…,Wm]; two matrices are of dimension n×mn. Suppose that the tridiagonal matrix Tm is denoted by (13)Tm=α1β2δ2α2⋱⋱⋱βm-1δmαm,where αi,βi, and δi (i=1,2,…,m) are the scalars defined in Algorithm 1.To derive the relation betweenAVm and Tm, we define the matrix T~m=Tmδm+1emT, where em=(0,…,0,1)T∈Rm.Recall the notation∗ in [10]:(14)Vm∗y=∑i=1m‍yiVi=Vmy⊗In,where y=(y1,y2,…,ym)T is a vector of Rm, In is the n×n identity matrix, and(15)Vm∗Tm=Vm∗T.,1,…,Vm∗T.,m,where T.,i denotes the ith column of the matrix Tm. Then, for j=1,2,…,m-1, we have the following relations: (16)Vm∗T.,j=αjVj+δj+1Vj+1+βjVj-1=AVj,Vm∗T.,m=βmVm-1+αmVm=AVm+δm+1Vm+1.Then we get (17)AVm=Vm∗Tm+δm+1On×n,…,On×n,Vm+1=Vm∗Tm+δm+1emTVm+1=Vm,Vm+1∗T~m.That is, by the global Lanczos biorthogonal process, we can obtain (18)AVm=Vm+1∗T~m.It was pointed out in [15] that the global Lanczos algorithm had significant advantages over the Arnoldi method for its fewer matrices of storage.


## 3. The Global QMR Method for Image Restoration

The quasi-minimal residual (QMR) method was first introduced by Freund and Nachtigal [19] to solve the linear equation Ax=b. The main idea of this algorithm is to solve the reduced tridiagonal system in a least squares sense. Additionally, the QMR method uses the look-ahead technique to avoid breakdowns in the underlying Lanczos process, which makes it more robust than the BiConjugate Gradient method (BiCG) [20], and when BiCG makes no progress at all, QMR may still show slow convergence. Since the linear system is usually of large scale in applications such as image restoration, it needs enormous computation.Fortunately, by applying the Kronecker product approximation of the matrixA [16, 17], the large-scale problems such as image restoration could be simplified intensively. In [10], the authors first introduced a global approach for solving matrix equations and derived the global FOM and the global GMRES methods. These methods are generalizations of the global MINRES method proposed by Saad [20]. The authors proved that these methods were effective when applied for matrix equations of large scale and low rank [21]. More recently, Wang and Gu [15] applied the global QMR method to solve the Sylvester equations. In this work, we will focus on the global QMR method for image restoration.Suppose that the operatorA=Ar⊗Ac is a good approximation of A, x=vec(X) and b=vec(B). In the following, we give details of the global QMR method for image restoration. Let X0∈Rn×n be the initial solution of (9) and let R0=B-AX0 be the corresponding residual. Usually we set the black image to be the initialization. By using Algorithm 1 for (9), the iterate Xm at step m satisfies that (19)Xm-X0∈KmA,R0.DefineV1=R0/ρ,ρ=R0F and W1 such that, 〈V1,W1〉F=1. Suppose that the matrix Krylov subspaces Km(A,R0) and Km(AT,W1) are generated by the sets of matrices {V1,V2,…,Vm} and {W1,W2,…,Wm} constructed by Algorithm 1. Then according to (19), we can obtain an approximate solution of (9): (20)Xm=X0+Vm∗ym,where ym=(y1,y2,…,ym)T. Consequently, we can get the associated residual matrix (21)Rm=B-AXm=B-AX0+Vm∗ym=R0-AVm∗ym=ρV1-Vm+1∗T~mym=Vm+1∗ρe1-T~mym,where e1=(1,0,…,0)T. Hence, the norm of the residual matrix is (22)RmF=Vm+1∗ρe1-T~mymF.An approximate solution of (9) can be obtained by computing the minimizer from (22) with respect to ym. Generally, the Vi’s obtained by Algorithm 1 are not F-orthonormal. However, as shown in [20], it is still reasonable to obtain that (23)y^m=argminym∈Rmρe1-T~mymF.What has been shown above is the key idea of the global QMR method; hence the approximated solution by the global QMR method can be given as(24)Xm=X0+Vm∗y^m,where y^m=argminym∈Rmρe1-T~mymF. We refer the readers to see [15, 20] for the details on how to compute y^m.To sum up, the algorithm for obtaining the approximate solution of (9) arising from image restoration can be described as in Algorithm 2.Algorithm 2:The global QMR method for image restoration.
(1) GivenA, and B,X0;
(2) SetR0=B^-AX0,ρ=‖R0‖F,W1=V1=R0/ρ;
(3) Form=1,2,…
Compute T~m,Vm,Wm using Algorithm 1;
y^m=argminym∈Rm⁡‖ρe1-T~mym‖2;
Xm=X0+Vm*y^m;
EndWe note that the discrepancy principle can be used as the stopping criterion in Algorithm2; that is, the computations will be terminated if the associated residual error corresponds to the approximate solution Xm(25)AXm-B≤μϵ,AXm-1-B>μϵ,where ϵ is the noise’s Frobenius Norm which is supposed to be a priori and μ≥1 is a fixed constant. For details on the discrepancy principle, we refer to [22] and references therein for more details.Since the matrixA is usually ill-conditioned, the solution is sensitive to the noise in the observed image. In the following, in order to improve the accuracy of the solution, we consider combining the global QMR method with Tikhonov regularization technique. Motivated by the work in [11], we can obtain the following algorithm which is named as the global Tik-QMR method.In Algorithm3, the regularization parameter λ can be determined by the L-curve criterion or the GCV method. We choose the latter here. Note that the regularization step in our work is different from the work in [11], since we adopt the regularization after Lanczos process while the authors in [11] used the regularization before the Lanczos process. Then our method needs fewer computations than theirs.Algorithm 3:The global Tik-QMR method for image restoration.
(1) GivenA,B^,X0,μ and ϵ;
(2) SetR0=B^-AX0,ρ=‖R0‖F,W1=V1=R0/ρ;
(3) Form=1,2,…
Compute T~m,Vm,Wm using Algorithm 1;
y^m=argminym∈Rm⁡‖ρe1-T~mym‖22+λ2‖ym‖22;
Xm=X0+Vm*y^m;
Rm=Vm+1*(ρe1-T~mym);
if ‖Rm‖F satisfies the discrepancy principle, stop
else
X0=Xm,ρ=‖Rm‖F,V1=Rm/ρ;
End


## 4. Numerical Experiments

In this section, we report some numerical examples to illustrate the performance of the global QMR method for image restoration problems. The results show that the quality of images restored by the global QMR method is better than those obtained by other methods of the same kind, such as the classic CGLS method [23] and the global GMRES method proposed by Jbilou et al. [10]. The experiments are carried out in Matlab 7.0 on a PC equipped with a 2.93 GHz Intel Core Duo CPU, with 2 GB of RAM, under Microsoft Windows XP.Example 1.
Our first example is to show the practical efficiency of the global QMR method. The original imageX of size 256×256 is shown in Figure 1(a), which can be obtained from theTelescope Science Institute, and intended to simulate a star cluster image taken by the Hubble space telescope before its defective mirror was replaced [24]. Let X denote the exact star cluster. The PSF used in this example is the so-called Moffat function [6]. The PSF is given by(26)Pij=i-162+j-1622-5if1≤i,j≤30,0,otherwise.
This PSFP is nonsymmetric and unseparable, so the blurring matrix A constructed from P is nonsymmetric and unseparable. If the zero boundary condition is imposed, the matrix A can be represented as the Kronecker product approximation of Toeplitz matrices Ar and Ac; that is, A=Ar⊗Ac. We add 1% white Gaussian noise to the blurred image to simulate the observed image (Figure 1(b)). The PSNR of the observed image is 30.82 dB. We set the parameter μ=1.05 for the discrepancy principle. Using the global QMR method, we obtained the estimated image after 4 iterations when the discrepancy principle of its associated residual is satisfied. The numerical results in terms of PSNR are reported in Table 1. From the table, we see that the PSNR of the restored image by the global QMR method is 41.16 dB and the consumed CPU time is 0.25 s. The restored image is shown in Figure 1(c).Table 1
Numerical results for the experiments, in terms of PSNR (dB) and CPU time (second).
Test problem
Method
restored PSNR
CPU time
observed PSNR
starcluster
global QMR
41.16
0.25 s
30.82
MRI
global QMR
20.61
0.156 s
17.95
global Tik-QMR
22.16
0.344 s
Indian man
global Tik-QMR
25.69
8.078 s
23.51
CGLS
25.05
1.7 s
bridge
global Tik-QMR
19.91
49.42 s
15.66
global GMRES
20.05
67.95 sFigure 1
Example1: (a) original image. (b) Observed image contaminated with blur and noise. (c) Image restored by the global QMR method.
(a)
Exact image
(b)
Observed image
(c)
Restored imageExample 2.
In order to suppress the sensitivity of solution to noise, we employ Tikhonov regularization technique to get a more accurate solution. In this example, we compared the performance of the global Tik-QMR method and the global QMR method. We consider the problem of restoring the image of the MRI data from Matlab (Figure2(a)). The data size is 128×128. The blurred and noisy image is shown in Figure 2(b). The PSF for blurring in this test is the truncated separable Gaussian function, and the variance of the Gaussian blur is 3 and 1% white Gaussian noise is added.
The restored images obtained by the global QMR method and the global Tik-QMR method are shown in Figure3, respectively. From Figure 3, it is easy to see that the image restored by the global Tik-QMR method has higher visual quality than that by the global QMR method. The numerical results are shown in Table 1, and it is not difficult to see that the global Tik-QMR method outperforms the global QMR method.Figure 2
Example2: exact and observed images.
(a)
Exact image
(b)
Observed imageFigure 3
Example2: restored images by the global QMR method and the global Tik-QMR method.
(a)
Global QMR image
(b)
Global Tik-QMR imageExample 3.
The third example consists in restoring the image of512×512 “Indian man” degraded by the Gaussian blur and 0.1% additive noise. The true image and degraded image are shown in Figure 4. We compare the global Tik-QMR method with the classic CGLS. The PSNR of the restored images by the two methods and computational CPU time are given in Table 1. The restored images are shown in Figure 5.
From Figure5 and Table 1, we see that the global Tik-QMR method is quite competitive with the CGLS method.Figure 4
Example3: exact and observed images.
(a)
Exact image
(b)
Observed imageFigure 5
Example3: exact and observed images.
(a)
Global Tik-QMR image
(b)
CGLS imageExample 4.
In the last experiment, the256×256 bridge image has been contaminated by a nonsymmetric wavefront blur [25] and 0.1% additive noise. The true image, the wavefront PSF, and the degraded image are shown in Figure 6. The PSF is also unseparable. Then the corresponding blurring matrix A is approximated by the Kronecker product of two small matrices. We compared the behavior of the global Tik-QMR method and the global GMRES method [11] in this experiment.
The numerical results are given in Table1. From the table, we see that the PSNR of the restored image by the global GMRES method is slightly higher than the global Tik-QMR method, but the CPU time by using the global Tik-QMR method is much less than the global GMRES method. The visual quality of restored images is very close. The restored images by using the two methods are displayed in Figure 7.
At the end of this section, a general comment about the presented numerical experiments is worth mentioning. The first example illustrates efficiency of the proposed method for image restoration problems. In general, the global Tik-QMR method behaves better than the classic CGLS method and the global GMRES method.Figure 6
Example4: exact image, wavefront PSF, and observed image.
(a)
Exact image
(b)
Wavefront PSF
(c)
Observed imageFigure 7
Example4: restored images using the global GMRES method and the global Tik-QMR method.
(a)
Global GMRES image
(b)
Global Tik-QMR image


## 5. Conclusion

In [10], Jbilou et al. first introduced the global methods. In this paper, we take the advantage of the global QMR method for image restoration and compare it with other popular methods. Numerical results show that the global QMR method is very efficient and is competitive with the classic CGLS method and the global GMRES method in [11]. In addition, when combining with the Tikhonov regularization, the global QMR method can behave much better.



---

*Source: 943072-2015-06-23.xml*