Global Quasi-Minimal Residual Method for Image Restoration

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 QMRmethod 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 QMRmethod outperforms the classic CGLS method and the global GMRES method.


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 ∬ Ω ℎ (, ; , )  (, )   =  (, ) , where the functions ,  represent the original and blurred images, respectively.The kernel ℎ is a point 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: where  is a vector representing the true image and b is a vector representing the blurred image, which are the discretized versions of  and  in (1), respectively.The matrix  is the blurring matrix constructed from the discretized version of the PSF ℎ.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 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 image  by computing a solution of the linear system of equations If the observed image array has dimension  × , then  and  are vectors of length  2 , and  is an  2 ×  2 matrix.Typical values of  are 256, 512, and 1024, so the dimensions of the matrix  can be extremely large [7].Then the computations with  can be very expensive.

Mathematical Problems in Engineering
Fortunately, the matrix  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  can be decomposed as Kronecker products, and then the computations with  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 matrix , 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.Section 3 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.

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  is employed.
Suppose that  ∈ R × is the discretized PSF.If the PSF is separable, that is,  can be decomposed as where a and b are  × 1 vectors, the matrix  constructed from  has block structure of the form where the matrices   and   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 PSF  is inseparable, then the corresponding matrix  is inseparable.However, we can find the Kronecker product approximation of  by using SVD technique so that  can be approximately decomposed as the following form: where  = ∑  =1 b  a   with a given integer  ≤ rank().In particular, the authors in [2,16] pointed out that  a 1 ⊗  b 1 is the best (as measured by the Frobenius norm) Kronecker approximation of .
According to the properties of Kronecker product, (4) can be rewritten as where  = vec() and  = vec().Note that  = vec() with  ∈ R × is the  2 ×1 vector obtained by stacking  columns of the matrix .Define an operator A :  ∈ R × →      and A  :  ∈ R × →      ; then (4) can be rewritten as We use the notation for the global Krylov subspace of R × generated by the matrix  ∈ R × and the operator A. Note that Let ,  ∈ R × ; we define the inner matrix product ⟨, ⟩  = tr(  ), where tr() denotes the trace of the square matrix  and   the transpose of the matrix .The associated norm is the Frobenius norm ‖⋅‖  .The matrices ,  are said to be F-orthonormal if tr(  ) = 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  1 ,  2 , . . .,   and  1 ,  2 , . . .,   of the two Krylov subspaces K  (A,  1 ) and K  (A  ,  1 ), respectively, such that The construction process can be summarized as in Algorithm 1.
To derive the relation between AV  and   , we define the matrix T = ( ), where   = (0, . . ., 0, 1)  ∈ R  .Recall the notation * in [10]: where  = ( 1 ,  2 , . . .,   )  is a vector of R  ,   is the  ×  identity matrix, and where  .,denotes the th column of the matrix   .Then, for  = 1, 2, . . .,  − 1, we have the following relations: Then we get That is, by the global Lanczos biorthogonal process, we can obtain It was pointed out in [15] that the global Lanczos algorithm had significant advantages over the Arnoldi method for its fewer matrices of storage.

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  = .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 matrix  [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 operator A =   ⊗   is a good approximation of ,  = vec() and  = vec().In the following, we give details of the global QMR method for image restoration.Let  0 ∈ R × be the initial solution of (9) and let  0 =  − A 0 be the corresponding residual.Usually we set the black image to be the initialization.By using Algorithm 1 for (9), the iterate   at step  satisfies that Define  1 =  0 /,  =      0     and  1 such that, ⟨ 1 ,  1 ⟩  = 1.Suppose that the matrix Krylov subspaces K  (A,  0 ) and K  (A  ,  1 ) are generated by the sets of matrices { 1 ,  2 , . . .,   } and { 1 ,  2 , . . .,   } constructed by Algorithm 1. Then according to (19), we can obtain an approximate solution of (9): where y  = ( 1 ,  2 , . . .,   )  .Consequently, we can get the associated residual matrix where  1 = (1, 0, . . ., 0)  .Hence, the norm of the residual matrix is An approximate solution of ( 9) can be obtained by computing the minimizer from ( 22) with respect to y  .Generally, the   's obtained by Algorithm 1 are not F-orthonormal.However, as shown in [20], it is still reasonable to obtain that ŷ = arg min 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 where ŷ = arg min   ∈R        1 − T y       .We refer the readers to see [15,20] for the details on how to compute ŷ .To sum up, the algorithm for obtaining the approximate solution of ( 9) arising from image restoration can be described as in Algorithm 2.
We note that the discrepancy principle can be used as the stopping criterion in Algorithm 2; that is, the computations will be terminated if the associated residual error corresponds to the approximate solution 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 matrix  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 Algorithm 3, 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.

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 image  of size 256 × 256 is shown in Figure 1(a), which can be obtained from the Telescope Science Institute, and intended to simulate a star cluster image taken by the Hubble space telescope before its defective mirror was replaced [24].Let  denote the exact star cluster.The PSF used in this example is the socalled Moffat function [6].The PSF is given by This PSF  is nonsymmetric and unseparable, so the blurring matrix  constructed from  is nonsymmetric and unseparable.If the zero boundary condition is imposed, the matrix  can be represented as the Kronecker product approximation of Toeplitz matrices   and   ; that is,  =   ⊗  .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.82dB.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).
Example 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 (Figure 2(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 Figure 3, 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.
Example 3. The third example consists in restoring the image of 512 × 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 Figure 5 and Table 1, we see that the global Tik-QMR method is quite competitive with the CGLS method.
Example 4. In the last experiment, the 256×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  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 Table 1.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.

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.

Figure 1 :
Figure 1: Example 1: (a) original image.(b) Observed image contaminated with blur and noise.(c) Image restored by the global QMR method.

Figure 3 :
Figure 3: Example 2: restored images by the global QMR method and the global Tik-QMR method.

Figure 7 :
Figure 7: Example 4: restored images using the global GMRES method and the global Tik-QMR method.

Table 1 :
Numerical results for the experiments, in terms of PSNR (dB) and CPU time (second).