Faster gradient descent and the efficient recovery of images

Introduction

The tasks of deblurring and denoising are fundamental in image restoration and have received a lot of attention in recent years; see, e.g., vogelbook ; chsh ; of ; dafrmo ; mallat and references therein. These are inverse problems, and they can each be formulated as the recovery of a 2D surface model $m$ from observed (or given) 2D data $b$ based on the equation

Here $F(m)$ is the predicted data, which is a linear function of the sought model $m$, and $\epsilon$ is additive noise. Both $m$ and $b$ are defined on a rectangular pixel grid, and we assume without loss of generality that this grid is square, discretizing the unit square $\Omega$ with $n^{2}$ square cells of length $h=1/n$ each. Thus, $m=\{m_{i,j}\}_{i,j=0}^{n}$.

where the sensitivity matrix $J=\frac{\partial F}{\partial m}$ is constant. Note that $N$ can easily exceed $1,000,000$ in applications; however, there are fast matrix-vector multiplication algorithms available to calculate $F$ for both deblurring vogelbook and (trivially) denoising.

The second form is really an extension of $m$ into a piecewise smooth function $m(x,y)$ defined on $\Omega$. This allows us to talk about integral and differential terms in $m$, as we proceed to do below, with the understanding that these are to be discretized on the same grid to attain their true meaning.

Both inverse problems are ill-posed, the deblurring one being so even without the presence of noise chsh . Some regularization is therefore required. Tikhonov-type regularization is a classical way to handle this, leading to the optimization problem

where $R(m)$ is the regularization operator and $\beta>0$ is the regularization parameter ehn1 ; ta . It is important to realize that the determination of $\beta$ is part of the regularization process. The least-squares norm $\|\cdot\|$ is used in the data-fitting term for simplicity, and it corresponds to the assumption that the noise is Gaussian. The necessary condition for optimality in (3) yields the algebraic system

The choice of $R(m)$ should incorporate a priori information, such as piecewise smoothness of the model to be recovered, which is the vehicle for noise removal. Consider the one-parameter family of Huber switching functions, whereby $R$ is defined as the same-grid discretization of

The corresponding necessary conditions (4) are the same-grid discretization of an elliptic PDE with the leading differential term given by

Thus, for $\gamma$ large enough so that always $\max\{\gamma,|{\boldsymbol{\nabla}}m|\}=\gamma$ the objective function is a convex quadratic and (4) is a linear symmetric positive definite system. However, this choice smears out image edges. At the other extreme, $\gamma=0$ yields the total variation (TV) regularization clmu ; rof ; vogelbook ; fo . But this requires modification when $m$ is flat to avoid blowup in $L$. For most examples reported here we employ the adaptive choice

proposed in ahh , which basically sets a resolution-dependent switch, modifying TV.

For the approximate solution of the optimization problem (3), consider combining two iterative techniques. The first iteration applies for the case where $\gamma$ is small enough so that $R_{m}$ is nonlinear in $m$, as is the case when using (7). In this case a fixed point iteration called lagged diffusivity (or IRLS) is employed vogelbook , where $m^{0}$ is an initial guess and $m^{k+1}$ is defined as the solution of the linear problem

for $k=0,1,2,\ldots$. See ahh ; chmu for a proof of global convergence of this iteration. In practice a rapid convergence rate is observed at first, typically slowing down only in a regime where the iteration process would be cut off anyway.

Next, for the solution of (3) or the potentially large linear system (8) consider the iterative method of gradient descent given by

Such a method was advocated for the denoising problem already in rof ; pema , and it corresponds to a forward Euler discretization of the embedding of the elliptic PDE in a parabolic PDE with artificial time $t$, written as

The step size $\tau_{k}$ in (9) had traditionally been determined for the quadratic case by exact line search, yielding the steepest descent (SD) method. But this generally results in slow convergence, as slow as when using the best uniform step size, unless the condition number of $J^{T}J+\beta L$ is small akaike ; nosazh ; asdohusv . The convergence rate then is far slower than that of the method of conjugate gradients (CG). In recent years much attention has been paid to gradient descent methods where the steepest descent step size value from the previous rather than the current iteration is used babo , or where it gets updated only once every second iteration rasv ; fmmr . These step selection strategies yield in practice a much faster convergence for the gradient descent method applied to convex quadratic optimization, although they are still slower than CG and their theoretical properties are both poorer and more mysterious asdohusv . Let us refer to these variants as faster gradient descent methods. These methods automatically combine an occasionally large step that may severely violate the forward Euler absolute stability restriction but yields rapid convergence with small steps that restore stability and smoothness. The resulting dynamical system is chaotic doas3 .

Faster gradient descent methods have seen practical use in the optimization context of quadratic objective functions subject to box constraints dafl , especially when applied to compressed sensing and other image processing problems finowr ; befr . For unconstrained optimization they are generally majorized by CG, and the same holds true for their preconditioned versions. Here, however, the situation is more special. For one thing, the PDE (10) has in the case of denoising, where $J$ is the identity, the interpretation of modeling the physical process of anisotropic diffusion (isotropic for sufficiently large $\gamma$). This has a smoothing effect that is desirable for the denoising problem. CG is also a smoother hansen , but it does not have the same physical interpretation. Moreover, CG is more susceptible to perturbation caused by the lagged diffusivity or IRLS method, where the quadratic problem solved varies slightly from one iteration to the next (see also doas3 ). Related and important practically, the iteration process (9) need not be applied all the way to strict convergence, because a desirable regularization effect is obtained already after relatively few iterations. To understand this, recall that the parameter $\beta$ in (3) is still to be determined. Its value affects the Tikhonov filter function and relates to the amount of noise present vogelbook . But it is well known that an alternative to the Tikhonov filter function is the exponential filter function—see, e.g., vogelbook ; care —and the latter is approximately obtained by what corresponds to integrating (10) up to only a finite time; see ashudo ; gigo and references therein. In fact, upon integration to finite time one may optionally set $\beta=0$, replacing the Tikhonov regularization altogether. The question now is, in the current problem setting, do faster gradient descent methods perform better than steepest descent, taking larger steps and yet maintaining the desired regularization effect at the same time? More specifically, does the more rapid convergence provided by the large steps and the regularization (i.e., piecewise smoothing) effect provided by the small steps automatically combine to form a method that achieves the desired effect in much fewer steps?

In asdohusv we started to answer this question for the denoising problem. Here we continue and extend that line of investigation. In Section 2 we complete the results presented in asdohusv and propose a new, efficient, explicit-implicit denoising scheme with an edge sharpening option.

In Section 3, the main section of this article, we apply the methods described above to the much harder deblurring problem, and also propose a new method for the case where there are both blur and significant noise present.

The two step size selections that are being compared for the deblurring problem can be written as

All presented numerical examples were run on an Intel Pentium 4 CPU 3.2 GHz machine with 512MB RAM, and CPU times are reported in seconds. Conclusions are offered in Section 4.

Denoising

For the denoising problem we have $F(m)=m$, so the sensitivity matrix $J$ is the identity, and it is trivially sparse and well-conditioned. In asdohusv we have considered the well-known gradient descent algorithm

obtained as a special case of (9) upon rescaling the artificial time by $\beta$ and then letting $\beta\rightarrow\infty$. This gets rid of the annoying need to determine $\beta$. The influence of the data is only through the initial conditions, i.e., it is a pure diffusion simulation. Correspondingly, in the step size definitions (11), $G$ is replaced by $R$ and $J^{T}J+\beta L$ is replaced by $L$.

The results in asdohusv clearly indicate not only the superior performance of the parameter selection (7) over using large $\gamma$ but also the improved speed of convergence using LSD over SD, occasionally by a significant factor. For instance, cleaning the Cameraman image used below in Fig. 3, which is $256\times 256$ and corrupted by $20\%$ Gaussian white noise, agreeable results are obtained after $382$ iterations using SD as compared to only $117$ iterations using LSD. These numbers arise upon using the same relative error norm, defined by

to stop the iteration in both cases. Unless otherwise noted we have used a sufficiently strict tolerance to ensure that the resulting images are indistinguishable. Similar results are obtained for other test images.

Still, these methods can often be further improved. With LSD, occasionally the larger step sizes could produce a slightly rougher image than desired (compare Figs. 4(f) and 4(e) in asdohusv ), whereas SD occasionally simply takes too long. We may therefore wish to switch to solving the Tikhonov equations (4), which here read

However, this brings back the question of effectively selecting the parameter $\beta$.

Fortunately, there is a fast way to determine $\beta$, at least for the case where the noise is Gaussian. Using (5) with (7), the reconstructed model becomes much closer to the true image than to the noisy one already after a few LSD or even SD iterations, see asdohusv . The computable misfit, defined by

therefore provides a very good, cheaply obtained approximation for the noise level. Thus, denoting the roughly denoised image by $\bar{m}$, we have $\int_{\Omega}|\bar{m}-b|^{2}\approx\eta^{2}$. According to the discrepancy principle, we wish our intended reconstruction $m$ to maintain this constraint invariant in “time”, that is

This therefore leads us to determine $\beta$ as

Now that we have $\beta$ (and $\bar{m}$ for an initial guess) we solve the equations (14) with (5) and (7), which corresponds to an implicit artificial-time integration step, using a combination of lagged diffusivity (IRLS) and CG with a multigrid preconditioner; see ahh for details.

In Figs. 1(d) and 1(e), we compare the result of this hybrid explicit-implicit scheme with that of the pure explicit scheme (12). The pre-denoised image after 17 SD steps, shown in Fig. 1(c), acts as a warm start and produces a good regularization upon setting $\beta=0.083$ by (16) for the following implicit process. Then, after only 3 IRLS iterations, the denoised image in Fig. 1(e) looks already quite comparable with—even slightly smoother than—the one in Fig. 1(d), which is continually denoised using SD for a stricter relative error tolerance. The processing time of the hybrid scheme is only a small fraction of that required by the explicit scheme, even with the faster LSD step size. Figs. 4(e) in asdohusv and 3(d) in the present article tell us essentially the same story for a different example.

Besides employing the implicit method to improve the quality of pre-denoised images, we may wish to sharpen the reconstructed image in a way that TV cannot provide sapiro . One possibility is to sharpen them by making $R(m)$ in (12) gradually more and more non-convex pema ; Black ; sapiro ; Dur and so reducing penalty on large jumps. This can be done by replacing the Huber switching function depicted in Fig. 2(a) with the Tukey function depicted in Fig. 2(b).

The Tukey function, scaled similarly to the Huber function, is defined by

To ensure that the Tukey and Huber functions start rejecting outliers at the same value we set

where $\gamma$ for the Huber function is still defined by (7).

Thus, the influence of the Tukey function decreases all the way to zero. Comparing the two functions $\phi(\sigma)=\sigma g(\sigma)$ in the horizontal center of Fig. 2, the Huber function gives all outliers a constant weight of one whereas the Tukey function gives zero weight to outliers whose magnitude is above a certain value. From such shapes of $\phi$ we can correctly predict that smoothing with Tukey produces sharper boundaries than smoothing with Huber. We can also see how the choice of edge-stopping function acts to hold off excessive smoothing: given a piecewise constant image where all discontinuities are above a certain threshold, Tukey will leave the image unchanged whereas Huber will not. Results in Figs. 1(f) and 3(e) confirm our predictions. After a quick pre-denoising with Huber, only 10 SD steps with Tukey result in obviously sharper discontinuities, i.e., image edges, than those without switching the regularization operator; see Figs. 1 and 3.

It is important to note that a rough pre-denoising is necessary here. It helps us avoid strengthening undesirable effects of heavy noise by using the Tukey function.

From an experimental point of view, we recommend to apply the explicit-implicit LSD scheme if a smooth image is desired; otherwise, employ the explicit Tukey regularization, starting from the result of the rough explicit Huber regularization, to get a sharper version.

Deblurring

Here we consider the deblurring problem discussed in vogelbook ; chsh ; Hardy . The blurring of an image can be caused by many factors: (i) movement during the image capture process, by the camera or, when long exposure times are used, by the subject; (ii) out-of-focus optics, use of a wide-angle lens, atmospheric turbulence, or a short exposure time, which reduces the number of photons captured; (iii) scattered light distortion in confocal microscopy. Mathematically, in most cases blurring can be linearly modeled to be shift-invariant with a point spread function (PSF), denoted by $f(x,y)$. Further, it is well known in signal processing and systems theory OppSch ; OppWil that a shift-invariant linear operator must be in the form of convolution, written as

So, an observed blurred image $b$ is related to the ideal sharp image $m(x,y)$ by

where $*$ denotes convolution product and the point spread function $f(x,y)$ may vary in space. Thus, the matrix-vector multiplication $Jm$ represents a discrete model of the distortion operator (18) convolved by the PSF.

In the spatial domain, the PSF describes the degree to which an optical system blurs or spreads a point of light. The PSF is the inverse Fourier transform of the optical transfer function (OTF). In the frequency domain, the OTF describes the response of a linear, position-invariant system to an impulse. The distortion is created by convolving the PSF with the original true image, see vogelbook ; chsh . Note that distortion caused by a PSF is just one type of data degradation, and the clear image $m_{true}$ generally does not exist in reality. This image represents the result of perfect image acquisition conditions. Nonetheless, in our numerical experiments we have a “ground truth” model $m_{true}$ which is used to synthesize data and judge the quality of reconstructions.

In the implementation we apply the same discretization as in Chapter 5 of vogelbook , and then the convolution (18) is discretized into a matrix-vector multiplication, yielding a problem of the form (1), (2), where $J$ is an $N\times N$ symmetric, doubly block Toeplitz matrix. Such a blurring matrix $J$ can be constructed by a Kronecker product. However, $J$ is now a full, large matrix, and avoiding its explicit construction and storage is therefore desirable. Using a gradient descent algorithm we actually only need two matrix-vector products to form $J^{T}(Jm)$, and there is no reason to construct or store the matrix $J$ itself. Moreover, it is possible to use a 2D fast Fourier transform (FFT) algorithm to reduce the computational cost of the relevant matrix-vector multiplication from $\mathcal{O}(N^{2})$ to $\mathcal{O}(N\log N)$. Specifically, after discretizing the integral operator (18) in the form of convolution, we have a fully discrete model

where $\epsilon_{i,j}$ denotes random noise at the grid location $(i,j)$. In general, the discrete PSF $\{f_{i,j}\}_{i,j=0}^{n}$ is 2D-periodic, defined as

where $\mathcal{F}$ is the discrete Fourier transform and $.*$ denotes component-wise multiplication.

Thus, we consider the gradient descent algorithm (9), starting from the data $m^{0}=b$, and compare the step size choices (11a) vs. (11b). As before, strictly speaking, these would be steepest descent and lagged steepest descent only if $L$ defined in (6) were constant, i.e., using least-squares regularization. But we proceed to freeze $L$ for this purpose anyway, which amounts to a lagged diffusivity approach.

One difference from the denoising problem is that here we do not have an easy tool for determining the regularization parameter $\beta$, and it is determined experimentally instead. However, for the problems discussed below this turns out not to be a daunting task.

To illustrate and analyze our deblurring algorithm, we generate some degraded data at first. We use the Matlab function fspecial to create a variety of correlation kernels, i.e., PSFs, and then deliberately blur clear images by convolving them with these different PSFs. The function fspecial(type,parameters) accepts a filter type plus additional modifying parameters particular to the type of filter chosen. Thus, fspecial(‘motion’,len,theta) returns a filter to approximate, once convolved with an image, the linear motion of a camera by len pixels with an angle of theta degrees in a counterclockwise direction, which therefore becomes a vector for horizontal and vertical motions (see Fig. 4); fspecial(‘log’,hsize,sigma) returns a rotationally symmetric Laplacian of Gaussian filter of size hsize with standard deviation sigma (see Fig. 6); fspecial(‘disk’,radius) returns a circular mean filter within the square matrix of side 2 radius+1 (see Fig. 7); fspecial(‘unsharp’,alpha) returns a $3\times 3$ unsharp contrast enhancement filter, which enhances edges and other high frequency components by subtracting a smoothed unsharp version of an image from the original image, and the shape of which is controlled by the parameter alpha (see Fig. 8); fspecial(‘gaussian’,hsize,sigma) returns a rotationally symmetric Gaussian low-pass filter of size hsize with standard deviation sigma (see Fig. 9); fspecial(‘laplacian’,alpha) returns a $3\times 3$ filter approximating the shape of the two-dimensional Laplacian operator and the parameter alpha controls the shape of the Laplacian (see Fig. 10).

All six images presented and used in this section are $256\times 256$. For the first four experiments, we only add a small amount of random noise into blurred images, say 1% ($\eta=1$), and stop deblurring when the relative error norm (13) is below $10^{-4}$. For a given PSF, the only parameter required is the regularization parameter $\beta$. From Fig. 4 we can clearly see that the smaller $\beta$ is, the sharper the restored solution is, including both image and noise. So the Boat image reconstructed with $\beta=10^{-3}$ in Fig. 4(c) still looks blurry and that cannot be improved by running more iterations. The Boat image deblurred with $\beta=10^{-5}$ in Fig. 4(e) becomes much clearer; however, unfortunately, such a small $\beta$ also brings the undesirable effect of noise amplification. The setting $\beta=10^{-4}$ seems to generate the best approximation of the original scene, and so it does in the following three deblurring experiments.

As in the case of denoising, the LSD step size selection (11b) usually yields faster convergence than SD. When $\beta$ is small, e.g., $10^{-4}$ or less, the step size (11a) is very close to the strict steepest descent selection obtained for a constant $L$, and so the famous two-periodic cycle of akaike (see also asdohusv ) appears in the step sequence, resulting in a rather slow convergence; see Fig. 5(b). The lagged step size (11b) breaks this cycling pattern (see Fig. 5(d)), providing a much faster convergence for the same error tolerance. Moreover, with the same parameter $\beta$, the reconstructed images using both SD and LSD are quite comparable, and it is difficult to tell any differences between them by the naked eye. For $\beta=10^{-5}$, 335 SD steps are required to reach a result comparable to Fig. 4(e) and the corresponding CPU time is 441.6 sec. Using LSD we only need 48.7 sec, reflecting the fact that the CPU time is roughly proportional to the number of steps required.

Figs. 6 – 8 reinforce our previous observations that the gradient descent deblurring algorithm with LSD step selection (11b) works very well, and that the improvement of LSD over SD is even more significant here than in the case of denoising. This is especially pronounced when a small value of $\beta$ must be chosen and when accuracy considerations require more than 20 or so steepest descent steps.

In asdohusv we have discussed other faster gradient descent methods. The half-lagged steepest descent (HLSD) method rasv ; fmmr was generally found there to consistently be at par with LSD, while other variants performed somewhat worse. In the present article we have applied HLSD in place of LSD for the above four examples, where the step size selection counts most. With HLSD the formula (11a) is applied only at each even-numbered step and then the same step size gets reused in the following odd-numbered one. The results were found to be again comparable to those using LSD, both in terms of efficiency and in terms of quality.

We have also experimented with a version of CG where we locally pretend, as in (11), that the minimization problem is quadratic. However, the CG method is well-known to be more sensitive than gradient descent to violations of its premises, and its iteration counts when applied to each of the examples in Figs. 4 – 7 were consistently about 20% higher than the better of LSD and HLSD.

The operation of deblurring essentially sharpens the image, whereas denoising essentially smooths it. Thus, trouble awaits any algorithm when both significant blur and significant noise are present in the given data set.

In our present setting, if we add more noise to the blurred images when synthesizing the data, say $\eta\geq 5$, then directly running the deblurring gradient descent algorithm as above may fail due to the more severe effect of noise amplification; see, e.g., Figs. 9(c) and 10(c). After a few iterations, the restored image can have a speckled appearance, especially for a smooth object observed at low signal-to-noise ratios. These speckles do not represent any real texture in the image, but are artifacts of fitting the noise in the given image too closely. Noise amplification can be reduced by increasing the value of $\beta$, but this may result in a still blurry image, as in Fig. 4(c). Since the effects of noise and blur are opposite, a quick and to some extent effective remedy is splitting, described next.

At first, we only employ denoising up to a coarser tolerance, e.g., $10^{-4}$. This can be carried out in just a few steps of (12) with either the SD or LSD step size selection. Starting with the lightly denoised image, as in Figs. 9(d) and 10(d), we next apply the deblurring algorithm (9) with the LSD step sizes to correct PSF distortion. Since now the noise level becomes higher, we slightly increase the value of the regularization parameter and apply $\beta=5\times 10^{-4}$ for the Toy example and $\beta=10^{-3}$ for the Pepper example. Observe that, even though some speckles still appear on the image cartoon components, the results presented in Figs. 9(e) and 10(e) are much more acceptable than those in Figs. 9(c) and 10(c), which were deblurred by the same number of iterations without pre-denoising. Finally, we can use the Tukey regularization in (12) to further improve the reconstruction, carefully yet rapidly removing unsuitable speckles and enhancing the contrast, i.e., sharpening. The results are demonstrated in Figs. 9(f) and 10(f). The total CPU time, given in the captions of Figs. 9 and 10, clearly shows the efficiency of this hybrid deblurring-denoising scheme.

Conclusions

In this paper we have examined the effect of replacing steepest descent (SD) by a faster gradient descent algorithm, specifically, lagged steepest descent (LSD), in the practical context of image deblurring and denoising tasks. We have also proposed several highly efficient schemes for carrying out these tasks, independently of the step size selection.

Our general conclusion is that in situations where many (say, over 20) steepest descent steps are required, thus building slowness into the solution procedure, the faster gradient descent method offers substantial advantages.

Specifically, four scenarios have been considered. The first is a straightforward denoising process using anisotropic diffusion asdohusv . Here the LSD step selection offers an efficiency improvement by a factor of roughly $3$.

In contrast, the second denoising scenario does not allow slowness buildup by SD because after a quick rough denoising we switch to an implicit method, with a good estimate for $\beta$ at hand, or to a sharpening phase using (17). The resulting method is new and effective, although not because of a dose of LSD.

Switching to the more interesting and challenging deblurring problem, the third scenario envisions the presence of little additional noise, so we directly employ the gradient descent method (9) with $J$ as described in Section 3, $R$ given by (5) and (7), and $\beta=10^{-4}$. This allows for slowness buildup when using SD step sizes, and the faster LSD variant then excels, becoming up to 10 times more efficient. The HLSD variant is overall as effective as LSD.

Finally, in the presence of significant noise effective deblurring becomes a harder task. We propose a splitting approach whereby we switch between the previously developed denoising and deblurring algorithms. This again creates a situation where LSD does not contribute much improvement over SD. The splitting approach has been demonstrated to be relatively effective, although none of the reconstructions in Fig. 10, for instance, is amazingly good. The problem itself can become very hard to solve satisfactorily, unless some rather specific knowledge about the noise is available and can be used to carefully remove it before deblurring begins.

A future problem to be considered concerns the case where the PSF causing blurring is not known.

References