Hazardous Continuation Backward in Time in Nonlinear Parabolic Equations, and an Experiment in Deblurring Nonlinearly Blurred Imagery

Identifying sources of ground water pollution, and deblurring nanoscale imagery as well as astronomical galaxy images, are two important applications involving numerical computation of parabolic equations backward in time. Surprisingly, very little is known about backward continuation in nonlinear parabolic equations. In this paper, an iterative procedure originating in spectroscopy in the 1930’s, is adapted into a useful tool for solving a wide class of 2D nonlinear backward parabolic equations. In addition, previously unsuspected difficulties are uncovered that may preclude useful backward continuation in parabolic equations deviating too strongly from the linear, autonomous, self adjoint, canonical model. This paper explores backward continuation in selected 2D nonlinear equations, by creating fictitious blurred images obtained by using several sharp images as initial data in these equations, and capturing the corresponding solutions at some positive time T. Successful backward continuation from t=T to t = 0, would recover the original sharp image. Visual recognition provides meaningful evaluation of the degree of success or failure in the reconstructed solutions. Instructive examples are developed, illustrating the unexpected influence of certain types of nonlinearities. Visually and statistically indistinguishable blurred images are presented, with vastly different deblurring results. These examples indicate that how an image is nonlinearly blurred is critical, in addition to the amount of blur. The equations studied represent nonlinear generalizations of Brownian motion, and the blurred images may be interpreted as visually expressing the results of novel stochastic processes.


Introduction
This paper presents an effective iterative procedure that can be used to solve a wide class of 2D nonlinear parabolic equations backward in time. However, previously unsuspected difficulties are also uncovered that may preclude useful backward continuation in parabolic equations deviating too strongly from the linear, autonomous, self adjoint, canonical model. Instructive 1D examples of ill-behaved continuations were previously reported in [6].
Continuation backward in time in parabolic equations is a notoriously ill-posed problem with some intriguing applications. Of major current interest are hydrologic inversion and image deblurring. Hydrologic inversion seeks to identify sources of groundwater pollution by backtracking contaminant plumes, [2][3][4][5]15,18]. This involves solving the advection dispersion equation (ADE) backward in time, given the contaminant spatial distribution g(x,y) at the current time T. In image science, images blurred by Let Ω be a bounded domain in R n with smooth boundary ∂Ω . Let L be a linear or nonlinear elliptic operator in Ω , acting on smoothly differentiable functions satisfying homogeneous Dirichlet or Neumannn conditions on ∂Ω . Let L be such that the forward initial value problem = , > 0, (0) = , t w Lw t w f is wellposed in 2 ( ) L Ω . Let 1 ( , ) w x t and 2 ( , ) w x t be any two solutions, and let Using logarithmic convexity techniques, [1,12,17], the folowing inequality can be established for a wide class of parabolic equations = t w Lw ,

Backward Continuity and the Hölder Exponent ( ) t µ
In many engineering or applied science contexts, only educated guesses would generally be available to estimate δ and M, rather than exact values. Typically, the L 2 relative error might be expected to be on the order of 1% or thereabouts. Since the given data ( ) f x may simultaneously approximate several distinct solutions ( , ) p w x t of Eq. (2) at time T, there are, in general, infinitely many possible solutions of Eqs. (2) and (3). If δ is small, it is generally assumed that any two such solutions would differ only slightly. The extent to which this expectation is justified depends on the decay behavior in the Hölder exponent ( ) t µ as illustrated in Fig. 1. In the best possible case, that of a linear self adjoint elliptic operator L with time-independent coefficients, we have ( ) = / t t T µ , so that ( ) t µ decays linearly to zero as continuation progresses from = t T to = 0 t . At = /2 t T , we have ( /2) = 1/2 T µ , and This indicates a loss of acccuracy from ( ) O δ to ( ) O δ , while still only half way to = 0 t . More typically, ( ) t µ is sublinear in t, possibly with rapid exponential decay. This can lead to much more severe loss of accuracy as reconstruction progresses to = 0 t . Such rapid decay of µ to zero can be brought about by various factors, including nonlinearity, non self adjointness, diffusion coefficients that grow rapidly with time, or adverse spectral properties in the elliptic operator L. In all cases, Eq. (4) does not guarantee any accuracy at = 0 t , but only provides the redundant information

Exponentially Decaying Hölder Exponent
The following is a simple example of a parabolic equation with exponentially decaying ( ) t With fixed > 0, In this linear self adjoint problem with growing time-dependent diffusion coefficient, ( ) t

Effective Backward Non-Uniqueness in Non Self Adjoint Problem
Reconstructing the correct backward solution from reasonably accurate data at some > 0 T , can be a major challenge even with slowly varying diffusion coefficients. The following counterexample was discovered computationally, using the parabolic solver methodology described in [6]. It involves a linear non self adjoint equation with variable coefficients, non-negative initial values, and non-negative solution.
Evidently, quite distinct initial values at = 0 t can produce almost identical solutions at = 1 t . This is a good example of effective backward non uniqueness. Indeed, if 0 ( ) red w x is the true solution in this example, the false solution 0 ( ) green w x would seem to be the more likely initial value, given the black = 1 t trace in Fig. 2. In ill-posed inverse problem computations, smoothness and non negativity of solutions are considered beneficial regularizing constraints. Here, both traces are smooth and non negative, satisfy the reasonable 2 L bounds in Eq. (10), and yet the ambiguity remains.  x .

2D Nonlinear Parabolic Equations and the Solution Operator T Λ
The computational experiments discussed below involve four images and two parabolic equations. Numerous other equations can be considered, and a large variety of unexpected phenomena are yet to be uncovered. Let Ω be the unit square 0 < , < 1 x y in the ( , ) x y plane. With fixed > 0 T , and homogeneous Neumann boundary conditions on ∂Ω , the following initial value problem will be studied, is well-defined on 2 ( ) L Ω . The nonlinear operator T Λ is not known explicitly. Rather, ( , , 0) T w x y Λ must be found by solving the appropriate initial value problem Eq. (11), or Eq. (13), and obtaining the corresponding solution at time T. Note that ( , , ) w x y T necessarily belongs to a very restricted class of smooth functions.
In the image deblurring experiments in Sec. 6, Eq. (11) will be used to blur the sharp MRI brain image (image A), and the sharp Marylin Monroe image (image D), by using these images as the initial data ( , ) g x y . The sharp USS Eisenhower image (image G), and the sharp Sydney Opera House image (image J), will be blurred using Eq. (13).

Continuation Backward in Time and the Van Cittert Iteration
In its original formulation, given the data ( ) f x and the explicitly known 1D linear convolution integral operator S with Fourier transform ˆ( ) > 0 S ω , the Van Cittert method solves = Sg f for the unknown ( ) g x by means of the iterative procedure Here, > 0 λ is a fixed relaxation parameter chosen so that , and the expectation is that m h g → . In fact, in spectroscopy and image processing applications [11,14], the Van Cittert method generally produces useful results after finitely many iterations, even though it may not converge.
We consider using this in the present parabolic context to recover ( , ,0) = ( , ) w x y g x y in Eqs. (11) with some fixed λ such that 0 < < 1 λ , and 1 ( , ) = N h x y is found to be a useful approximation to ( , , 0) w x y . As noted in Sec. 3, backward continuation in certain classes of parabolic equations can be especially challenging. Accordingly, interesting continuation problems may exist where the procedure in Eq. (17) cannot produce useful results.

Explicit Finite Difference Scheme for Computing
( , , 0) T w x y Λ A convenient and effective numerical procedure for solving the nonlinear initial value problems in Eqs. (11) and (13), is based on finite differences, using explicit time diffencing and centered space differencing.
This leads to modest ( ) ( , ) = f x y W . We stress that the blurred image ( , ) f x y so obtained is only an approximation to the true solution ( , , ) w x y T in Eq. (11) or Eq. (13). Image diagnostic statistical information will use the discrete 1 2 , L L , and total variation norms, defined by .
In addition, the peak signal to noise ratio ( ) PSNR , will be used as an image quality metric. If ( , ) g x y is the original sharp image, and ( , ) f x y is any degraded version of ( , ) g x y , this is defined by

Nonlinear Blurring and Deblurring Experiments
Very little seems to be currently known regarding backward in time continuation in multidimensional nonlinear parabolic equations, and the experiments described below, involving relatively simple nonlinearities, already represent uncharted waters. An important advantage of the Van Cittert method is the 'self regularizing' property of the iterative process, whereby low frequency information is reconstructed in the first few iterations, while many more iterations are needed to acquire high frequency information. Several other iterative restoration methods have this property. As a consequence, useful information can often be retrieved by stopping the iterative process before amplification of high frequency noise overwhelms the reconstruction.
The results developed in Sec. 3, concerning backward stability and the Hölder exponent ( ) t µ , will inform the subsequent discussion. While backward uniqueness is characteristic of a large class of linear and nonlinear parabolic equations, the major practical difficulty lies in recovering the correct solution from the limited precision available in the given continuation data. Deblurring nonlinearly blurred imagery involves the recovery of fairly complex initial data at time = 0 t , by nonlinear backward continuation of imprecise data at some > 0 T . In fact, as is typically the case in applications, the accuracy in the blurred image data 400 ( , ) = f x y W in Eq. (18), as an approximation to the elusive true solution ( , , ) w x y T in Eq. (11) or Eq. (13), is actually unknown.
In addition, Eqs. (11) and (13) are strongly nonlinear through the functions ( ) r w and ( ) s w , with ( , , ) w x y t ranging between 0 and 255. Moreover, there is the space and time dependent function ( , , ) q x y t , and the terms in x ww and y ww . Such equations deviate strongly from the autonomous, linear, self adjoint case, for which substantial computational experience exists. While a stability inequality such as Eq. (4) can be derived for each of Eqs. (11) and (13), the resulting functional form for the Hölder exponent ( ) t µ is unlikely to be precise. In summary, neither δ nor ( ) t µ are likely to be known in the stability estimate for either of Eqs. (11) or (13). The results in Figs. 1 and 2, together with the examples in [6], indicate that only a modest degree of success can be expected in backward continuation in Eqs. (11) and (13). In the present paper, knowledge of the original sharp image can be used to gauge the usefulness of the deblurred image produced by backward continuation. However, in applications unrelated to imaging, using field data of unknown precision, the degree of success or failure in nonlinear backward continuation may not be as easily ascertained. As shown in [6], there is the possibility of producing a smooth, physically plausible, yet false reconstruction. http://dx.doi.org/10.6028/jres.118.010

MRI Brain Image
In Fig. 3, the original sharp MRI brain image (A) is blurred to form image (B), by applying the finite difference scheme in Eq. (18) to the parabolic equation Eq. (11), with coefficients = = 0 a b . A different blurred image is then obtained, image (C), by repeating this process with coefficients = 1.25, = 0.6 a b . Images (B) and (C) appear very similar in quality, and, from Table 1, both these images have almost the  same values for   1 2 In particular, 1 f ∇   has been reduced by almost a factor of two from its original value in image (A), reflecting substantial blurring. The PSNR value in image (C) is noticeably smaller than in image (B), indicating greater degradation in image (C). However, since the PSNR metric requires knowledge of the original sharp image, in practice, such increased degradation in image (C) would not be known to a user. In fact, both images (B) and (C) appear to have been blurred, more or less equally, by convolution with a type of Gaussian-like point spread function.    Fig. 3.   The results of deblurring image (C), shown in in Fig. 4, are sharply different, and unexpected. After 10 Van Cittert iterations, using = 1.25, = 0.6 a b , in Eq. (18), most useful information has been lost in the deblurred image, and this without explosive noise amplification. Indeed, in Table 2, the values for 1 f   and 2 f   after 10 iterations are about 12 % larger than their true values in the original image (A), well within the range of what might have been prescribed to stabilize continuation. One possible explanation is that the inclusion of the terms in x ww and y ww in Eq. (11) renders backward stability more precarious, (see Fig. 1), and the accuracy in the continuation data represented by image (C) is no longer sufficient to recover the sharp image.

Marilyn Monroe Image
In view of the unexpected failure in deblurring image (C), the experiments in Figs. 5 and 6 aim at elucidating the influence of the nonlinear terms involving  Table 3, both these images have almost the same values for 1 2 , , f f     and 1 f ∇   . However, image (F) has a smaller PSNR value. Evidently, the x ww term in Eq. (11) is responsible for the increased degradation in image (F).    Figure 6 displays the results of deblurring these two images. Again, remarkably, despite the strongly nonlinear blurring in image (E) through the function ( ) r w , and the inclusion of the Figure 8 displays the results of deblurring these two images. Despite the strongly nonlinear blurring in image (H), through the function ( ) s w and the inclusion of both the | | x w w and y ww terms in Eq. (13), reasonably good results are obtained after 100 Van Cittert iterations. The carrier's command 'island' has been recovered, along with the two rows of planes on deck. As shown in Table 6, the values of  Table 6. As was the case in the deblurred Marilyn Monroe image (F), substantial sharpening has occurred in the deblurred image (I), but the sharpened image is seriously marred by artifacts. Because of the moderating effect of the factor 2 cos w , it was not anticipated that the term

Concluding Remarks
The successful deblurring of images (B), (E), (H), and (K), indicate the Van Cittert iterative procedure to be a useful tool for backward in time continuation in an important class of 2D nonlinear parabolic equations. A wide variety of nonlinear problems remains to be explored. Surprisingly, the simple nonlinearities in Eqs. (11) and (13) involving the terms in ww x and ww y , were found to be potentially destabilizing, and capable of preventing useful continuation. The practical difficulty of reconstructing the correct backward solution, using data of limited and unknown precision, was stressed. Other limitations include the fact that the fundamental stability inequality governing a particular continuation, Eq. (4), can seldom be obtained with sufficient precision. In particular, the rate at which the Hölder exponent ( ) t µ tends to zero as 0 t ↓ , which is of vital interest, is typically unknown. The unsuccessful deblurring in images (C), (F), (I), and (L), is of major interest. Visually, the amount of blurring in each of these images is no greater than in the successfully deblurred companion image, and the corresponding values of 1 f ∇   are almost equal in every case. Evidently, how the image was blurred is critical, not just the amount of blur. These failed continuations suggest that the presence of ww x and ww y terms in Eqs. (11) and (13), with relatively large coefficients, unexpectedly leads to faster decaying Hölder exponents ( ) t µ , such as is shown in Fig. 1, and the accuracy in the blurred image data becomes insufficient for useful continuation.