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.
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 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 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 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 , and is an 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.
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.
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 is employed.
Suppose that is the discretized PSF. If the PSF is separable, that is, can be decomposed as where and are vectors, the matrix constructed from has block structure of the form where the matrices and have parameters and , 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 with a given integer . In particular, the authors in [2, 16] pointed out that 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 and . Note that with is the vector obtained by stacking columns of the matrix . Define an operator and ; then (4) can be rewritten as We use the notation for the global Krylov subspace of generated by the matrix and the operator . Note that
Let ; we define the inner matrix product , where 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 .
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 and of the two Krylov subspaces and , respectively, such that
The construction process can be summarized as in Algorithm 1.
| ||||||||||||||||||||||
For convenience, we denote by and the block matrix, that is, and ; two matrices are of dimension . Suppose that the tridiagonal matrix is denoted by where , and () are the scalars defined in Algorithm 1.
To derive the relation between and , we define the matrix , where .
Recall the notation in [10]:where is a vector of , is the identity matrix, andwhere denotes the th column of the matrix . Then, for , 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.
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 . 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 is a good approximation of , and . In the following, we give details of the global QMR method for image restoration. Let be the initial solution of (9) and let 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 and such that, . Suppose that the matrix Krylov subspaces and are generated by the sets of matrices and constructed by Algorithm 1. Then according to (19), we can obtain an approximate solution of (9): where . Consequently, we can get the associated residual matrix where . 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 . Generally, the ’s obtained by Algorithm 1 are not F-orthonormal. However, as shown in [20], it is still reasonable to obtain that
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 . 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 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.
| ||||||||||||||||||||||
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 image of size 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 so-called Moffat function [6]. The PSF is given byThis 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 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 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).
(a) Exact image
(b) Observed image
(c) Restored image
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 . 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 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.
(a) Exact image
(b) Observed image
(a) Global QMR image
(b) Global Tik-QMR image
Example 3. The third example consists in restoring the image of “Indian man” degraded by the Gaussian blur and 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.
(a) Exact image
(b) Observed image
(a) Global Tik-QMR image
(b) CGLS image
Example 4. In the last experiment, the bridge image has been contaminated by a nonsymmetric wavefront blur [25] and 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.
(a) Exact image
(b) Wavefront PSF
(c) Observed image
(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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is supported by NSFC (61370147, 61170311), 973 Program (2013CB329404), and Sichuan Province Sci. & Tech. Research Project (2012GZX0080).