64 kB PostScript reprint
N. Wu
Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218
The IRAF mem0 package, version B, was released to the public in
October 1992, and was upgraded to version C in December 1993. Version C
resides in the anonymous ftp site at NOAO ( iraf.noao.edu, directory
/contrib/).
Now, the most important MEM task of the mem0 package has been
upgraded to version D; the other four tasks remain unchanged.
This task is included in the subpackage stsdas.analysis.restore
with the name mem. Its executable code on SUN/UNIX has also been put
into the ftp site at NOAO with the task name dmem. The two relevant
files are dmem.notes and dmem.tar.Z.
The task mem, version D, features the following sophisticated functions and techniques:
CCD images from the Wide Field and Planetary Camera (WF/PC) of
Hubble Space Telescope ( HST) have both readout noise, of Gaussian type, and
signal noise, of Poisson type. In this case the total noise variance is
calculated using the parameters noise, in electrons, for the former,
and adu (analogue-to-digital unit conversion constant or gain),
in electrons/DN, for the latter. Then this variance is used in the
Gaussian likelihood function or, equivalently, in the expression of
.
On the other hand, images from the photon counting detectors of the Faint Object Camera (FOC) of HST contain only signal noise of Poisson type. Not only should the parameter noise be set to zero, but, more importantly, the Poisson likelihood function must be used in calculation. Failure to do so will lead to wrong and unacceptable results.
Earlier versions of the task mem could not handle correctly the case of Poisson noise only. In version D, setting the parameter poisson to yes selects the correct expression for the likelihood function, automatically sets the parameter noise to zero, and takes other actions especially designed for the case of Poisson noise only.
It has been shown that the subpixelization (subsampling) technique may improve resolution in restored images, or at least result in more pleasant appearances for objects (Weir & Djorgovski 1990).
Subpixelization means to restore an image on a grid finer than that on which the input data image is defined. In this process one ``normal size'' pixel of the input image is ``split'' to several ``subpixels'' of the restored image. This mechanism must be built into the program, but cannot be accomplished simply by replicating each pixel of the input image before restoration.
Subpixelization is activated by setting the parameter nsub
to a value greater than one. In this case all input images, including
the point spread function (PSF) but excluding the data image, must be
subpixelized by the user by a factor of nsub in each dimension.
The required core memory and computational time are approximately
proportional to nsub
.
This technique was described in detail in Wu (1994).
The MEM program uses a double iteration scheme: the values of the
Lagrange multipliers 
and 
, respectively for
the data constraint and the total power constraint, are revised
in the outer iteration, while the inner iteration is
for finding the ME solution for the particular 
and
of each outer iteration.
The basic idea behind the model updating technique is to use the ME image
converged for particular 
and 
in each outer iteration as
the model to start the next iteration. In this way in maximization,
the approximate solution to the linear equations used in the
preconditioned conjugate method is more accurate, or the accurate solution
to the linear equations required in the accurate Newton method is
easier to be found. Therefore, the total number of iterations is
considerably reduced and much computational time is saved.
The restored image also has improved photometric linearity.
In previous versions of the task mem, only the zeroth-order approximate Newton method of maximization was available. Hence the name of the package: mem0. Zeroth-order approximation means that in the solution of a large set of linear equations (or equivalently the inversion of a large matrix) non-diagonal elements are ignored, under the assumption that the diagonal ones dominate. In this way, solving the equations becomes a simple operation. However, this simple method may result in very slow convergence.
Now, in version D, much more sophisticated methods are used to calculate
the change in the iteration, i.e., to determine the search direction in
maximization.
The first method is the preconditioned conjugate method (as opposed to
the ``standard'' conjugate method commonly used), or the conjugate
method based on the approximate Newton method described in the above.
More specifically, in each outer iteration for particular 
and
, the approximate Newton method is used to calculate
the search direction for the first inner iteration. Thereafter,
in each inner iteration, a direction is calculated using the
approximate Newton method. The component orthogonal or conjugate to
the search direction used in the previous inner iteration is calculated,
and used as the new search direction for the maximum point.
By careful programming, the core memory requirement and the number of FFTs
are the same as in the approximate Newton method.
The second method is the accurate Newton method (as opposed to the approximate Newton method). Here the accurate solution to a set of linear equations of large size is calculated by iteration, each of which requires two convolutions, i.e., four FFTs.
The accurate Newton method is the most efficient in the sense that the fewest total number of (inner) iterations is needed because the search direction is determined accurately. However, many FFTs may well be required in each (inner) iteration to solve the linear equations. In contrast, the preconditioned conjugate method requires a greater total number of (inner) iterations, but no extra FFTs in each (inner) iteration are needed to calculate the conjugate direction.
By default, in the case of Poisson noise only ( poisson=yes) the accurate Newton method is used, otherwise ( poisson=no) the preconditioned conjugate method is used. This is the best choice. In the case of Poisson noise only, the Hessian of the objective function, which is a measure of the curvature of the image space, changes rapidly from (inner) iteration to iteration. Consequently, the preconditioned conjugate method, as a method taking advantage of memory in the iteration, is not effective in determining the search direction.
After determining the search direction, an optimal step (length) in this direction should be calculated. In earlier versions of the task mem, quadratic extrapolation and cubic interpolation are used for this purpose. They are both approximate methods for calculating the optimal step. Now, in version D, the accurate one-dimensional (1-D) search is available for interpolation. Specifically, the approximate maximum point found by cubic interpolation is used as the initial guess to start a search for the accurate maximum point, using the Newton method in a single variable. In such a way, the maximum point is found with little extra effort but much higher accuracy.
The (approximate) quadratic extrapolation remains. In most cases, especially at the late stage of iteration, interpolation but not extrapolation is desirable in the 1-D search.
Apart from employing the methods described in the previous subsections to enhance the task's function and to speed up convergence, having a good user interface is also important. Every effort has been made to create a user friendly interface. The number of positional parameters is kept to a minimum. The default values of hidden parameters are carefully chosen. Automatic schemes to adjust some variables and to deal with some predicted ill-conditioned cases are built in the program. Options and parameters used only for testing are hidden from the user. The diagnostic messages are informative and grouped logically, and can be output at three different levels of verboseness to meet the user's needs. The most detailed level ( message=3) is primarily for debugging purposes. The help file is well written and should be read before the first attempt to run the task.
The current version of the task mem has its limitations: Like other MEM programs, it can only handle the case of space-invariant PSF. It cannot be used for multi-channel and multi-data set restoration like the package MEM/MemSys5 (Weir 1991). The algorithm used in the accurate Newton method to find the accurate solution of a large set of linear equations should be improved, e.g., using the preconditioned conjugate method. Finally, criteria for convergence should be investigated, especially in the case of Poisson noise only, and more reasonable ones should be adopted. For this version of the task, the user's judgment is very important in obtaining satisfactory results.
Weir, N. 1991, in Proceedings of the 3rd ESO/ST-ECF Data Analysis Workshop, ed. P. J. Grosbøl & R. H. Warmels (Garching, ST-ECF), p. 115
Wu, N. 1994, in The Restoration of HST Images and Spectra II, ed. R. J. Hanisch & R. L. White (Baltimore, Space Telescope Science Institute), p. 58
64 kB PostScript reprint