# Applications of a nonnegatively constrained iterative method with statistically based stopping rules to CT, PET, and SPECT imaging.

1. Introductory material. The medical imaging modalities of computed tomography (CT), positron emission tomography (PET), and single photon emission computed tomography (SPECT) are widely used in the medical professions. Because of the unique set of strengths and weakness of each of these methods, they are utilized in different settings. However, mathematically they are closely related.

In this introductory section, we set the stage for the algorithmic discussion that constitutes the main result of the paper by briefly introducing the mathematical models for CT, PET, and SPECT imaging, their numerical discretization, and the associated statistical models. The mathematical discussion is included because it is accessible and will be of interest to the unfamiliar reader, while the statistical models are integral to the development of the computational methods that are our focus.

1.1. Mathematical models. CT is the most widely used of the three methods and has the simplest mathematical model. In the two-dimensional case, which is our focus in this paper, CT involves the reconstruction of the mass absorption function [mu] of a body from one-dimensional projections of that body. A particular one-dimensional projection is obtained by integrating [mu] along all parallel lines making a given angle w with an axis in a fixed coordinate system. Each line L can be uniquely represented in this coordinate system by [omega] together with its perpendicular distance y to the origin.

Suppose L([omega],y) = {x(s)| 0 [less than or equal to] s [less than or equal to] S}, with an X-ray source located at s = 0 and a sensor at s = S. The standard assumption is that the intensity I of the X-ray along a line segment ds is attenuated via the model [6]

dI = -[mu](x(s))I ds.

The resulting ordinary differential equation can be solved using the method of separation of variables to obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where I(0) is the intensity at the source and I(S) is the intensity at the receiver. Setting z = - ln(I(S)/I(0)), we obtain the Radon transform model for CT:

(1.1) z([omega], y) = [[integral].sub.L([omega],y)] [mu](x(s))ds.

A discretized version of (1.1) is what is solved in the CT inverse problem, where z corresponds to collected data, and [mu] is the unknown.

In both PET and SPECT, a somewhat different problem is solved. In both cases, a radioactive tracer element is injected into the body. The tracer then exhibits radioactive decay, resulting in photon emission. The emitted photons that leave the body are, in theory, all recorded by a photon detector, which also determines the line of response (LOR) L([omega], y). In PET, two photons are emitted, and if they reach the detector ring at the same time (i.e., within a few nanoseconds of each other), an event is registered along the connecting line, which is also the LOR. In SPECT, single photons are emitted and detected. The LOR is determined using a method known as columnation. In both cases, the task is to reconstruct the tracer density distribution u within the subject given the collected photon count data.

In both PET and SPECT, the data z([omega], y) for the line L([omega], y) correspond to the number of detected incidents along that line. The model relating the tracer density u to the data is similar to (1.1) and is given by

z([omega], y) = [[integral].sub.L([omega],y)] [g.sub. [omega],y](x(s))u(x(s))ds,

where the impulse response function [g.sub. [omega],y](x(r)) can be viewed as the probability that an emission event located at x(r) along L([omega], y) is detected.

For SPECT, we assume that the detector is located at x(S). Then, making assumptions for probabilities analogous to those made for intensities in the derivation of the CT model above, and assuming probability of 1 at x(r), one can obtain [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], so that the full model becomes

(1.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where the interior integral is along the line L([omega], y).

For PET, since there are two photons that have to reach the respective detectors at x(0) and x(S), the impulse response is the product of the probabilities,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

which does not depend on r and hence we have the somewhat simpler mathematical model

(1.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

A discussion of the PET mathematical and statistical models is given in [22].

Note the appearance of the absorption density [mu] in both (1.2) and (1.3), which must be known beforehand in order to solve the PET and SPECT inverse problems. Estimates of [mu] can be obtained using, for example, CT. Note also that brief derivations of all of these models can be found in [14], whereas [6] focuses on CT, [15] on PET, and [4] on SPECT.

1.2. Numerical discretization. After discretization, each of (1.1), (1.2), and (1.3) can be written as a system of linear equations. The discretization occurs both in the spatial domain where [mu] and u are defined, as well as in the Radon transform domain where the independent variables are [omega] and y. We will use a uniform n x n spatial grid, and a grid for the transform domain with [n.sub.w] angles and [n.sub.s] sensors, both uniformly spaced. Then, after column-stacking the resulting two-dimensional grids, we obtain a matrix-vector system

(1.4) z = Au,

with the data vector z [member of] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], the unknown vector u [member of] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and the forward model matrix A [member of] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

In the case of CT, for the i-th line [L.sub.i], we have the discrete model

[z.sub.i] = [[n.sup.2].summation over (j=1)][a.sup.Radon.sub.ij] [[mu].sub.j],

where [a.sup.Radon.sub.ij] is the intersection length of [L.sub.i] with pixel j. Written as (1.4), we have [[A].sub.ij] = [a.sup.Radon.sub.ij] and [u.sub.j] = [[mu].sub.j].

For both PET and SPECT, the system of linear equations instead has the form

where [g.sub.ij] is the discrete impulse response function, which for PET is

(note no dependence on j) and for SPECT is

[g.sub.ij] = exp (-[[n.sup.2].summation over (k=j)][a.sup.Radon.sub.ik][[mu].sub.k]).

Written as (1.4), we have [[A].sub.ij] = exp [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for PET and [[A].sub.ij] = exp [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for SPECT.

1.3. Statistical models. The character of the noise in the data z is an important consideration, particularly in the cases of PET and SPECT.

In CT, because intensities are typically high, a statistical model of the form

(1.5) z = Au + n,

where n is a zero mean, independent and identically distributed normal distribution with variance [[sigma].sup.2], is not uncommon. The resulting likelihood function takes the form

(1.6) [p.sub.z](z; u) [varies] exp (-1/2[[sigma].sup.2][[parallel]AU - 2[parallel].sup.2]).

PET and SPECT data are typically much more noisy, and hence it is important to accurately model noise statistics in these cases. We follow [4, 15], and use a Poisson statistical model

(1.7) z = Poiss(Au + [gamma]).

In this case, the likelihood function takes the form

(1.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Given image data z arising from model (1.5) or (1.7), the maximum likelihood estimate of u is obtained by maximizing [p.sub.z](z; u) with respect to u [greater than or equal to] 0, or equivalently by solving

(1.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where for likelihood (1.6),

(1.10) T(u; z) = 1/2[[sigma].sup.2] [[parallel]Au - z[parallel].sup.2],

whereas for likelihood (1.8),

(1.11) T(u; z) = [n.summation over (i=1)]{[([Au].sub.i] + [[gamma].sub.i]) - [z.sub.i]ln([[Au].sub.i] + [[gamma].sub.i])}.

Solutions of (1.9) tend to be noise-corrupted due both to random errors in z and ill-conditioning of the matrix A, which in all cases is the discretization of a compact operator [13]. Thus, some form of regularization is needed. One way that this can be accomplished is by truncating an iterative method applied to (1.9) [5]. When this approach is taken, the choice of stopping the iteration becomes extremely important and is akin to the choice of a regularization parameter in the standard Tikhonov approach to regularization [5, 23].

In the context of PET imaging, the Richardson-Lucy algorithm (RL) is widely used. RL is an iterative method for solving (1.9), (1.11) [10, 19], and stopping rules for its iterations have been given [9, 16, 17, 18]. Much effort has gone into the development of efficient iterative methods for PET; see, e.g., the review paper [15] and the references therein and also [8]. In this paper, we apply the iterative method and stopping rules of [1, 3], developed in the context of astronomical imaging, to (1.9) for both (1.10) and (1.11). The iterative method, called modified residual norm steepest descent (MRNSD) in [12], was originally introduced in [7] in the context of PET.

The remainder of the paper is organized as follows. In Section 2, we present the iterative method and stopping rules that we will use for the image reconstruction step, and then in Section 3 we test our methods on some synthetically generated examples. Finally, we end with conclusions in Section 4.

2. An iterative method and stopping rules. It is now time to present our numerical method and iteration stopping rules for approximately solving (1.9). The methods require a least-squares formulation of the nonnegatively constrained optimization problem (1.9). For noise model (1.5) and T defined by (1.10), this is already the case. However, for noise model (1.7) an approximation of (1.11) must be computed. For this, we follow [1, 2], where a Taylor series argument is used to motivate the following quadratic approximation of T:

T(u; z) [approximately equal to] T([u.sub.e]; z) + [T.sub.wls](u; z),

where

(2.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

with C = diag(z). Here we assume that z > 0. This leads to the following approximation of(1.9):

(2.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Before continuing, we note that T defined by (1.10) also can be written as (2.1), with C = [[sigma].sup.2]I and [gamma] = 0. In what follows we will use this convention in order to unify notation. However, we note that this implies that the user knows the noise level [[parallel]n[parallel].sup.2] in (1.5) a priori.

2.1. The iterative method. The algorithm of [3] for the numerical solution of (2.2) has the form

(2.3) [u.sub.k+1] = [u.sub.k] - [[tau].sub.k][u.sub.k] [dot encircle] [A.sup.T][C.sup.-1]([Au.sub.k] - (z - [gamma])),

where "[dot encircle]" denotes Hadamard (component-wise) multiplication, and the line search parameter [[tau].sub.k] in (2.3) is given by

[[tau].sub.k = min{[[tau].sub.uc], [[tau].sub.bd]}.

Here, [v.sub.k] = [u.sub.k] [dot encircle] [nabla][T.sub.wls]([u.sub.k]) and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This method has been shown to be effective on several astronomical imaging examples. For more details, see [1, 3]. In this paper, we will refer to this technique as weighted MRNSD, or WMRNSD.

2.2. The stopping rules. Next, we introduce three stopping rules for (2.3). Each can be motivated from the assumption that our statistical model has the form

(2.4) z - [gamma] = Au + n,

where n is distributed as a zero-mean Gaussian random vector with covariance matrix C. For (1.5), [gamma] = 0, C = [[sigma].sup.2]I, and model (2.4) is exact. For (1.7), on the other hand, we take C = diag(z), and the negative-log likelihood function (2.1) results. Note that alternative to the mathematical derivation of (2.1) above, one can use the statistical approximation (2.4) to derive the approximate maximum likelihood problem (2.1). More detailed derivations of these stopping rules can be found in [1].

2.2.1. The discrepancy principle. Assuming that (2.4) holds, [C.sup.-1/2] ([A.sub.ue] - (z - [gamma])), where [u.sub.e] is the true image, is approximately normally distributed with covariance I. Then by a standard result

[2T.sub.wls]([u.sub.e]) ~ [chi square](n),

where [chi square](n) denotes the chi-squared distribution with n degrees of freedom, and hence

E([T.sub.wls]([u.sub.e]))[approximately equal to])) 1/2 E([chi square](n)) = n/2.

Since in early iterations 2/n [T.sub.wls]([u.sub.k]) is typically much larger than 1, a stopping rule of the form

(2.5) 2/n [T.sub.wls]([u.sub.k]; z) [less than or equal to] 1 + [[member of].sub.n]

is therefore reasonable. We note that [[member of].sub.n] = 0 corresponds to Morozov's discrepancy principle [11, 23], and we recommend its use unless (2.5) isn't satisfied in a feasible number of iterations or if it uniformly yields over-regularized reconstructed images. Otherwise, [[member of].sub.n] can be taken to be, for example, [+ or -][square root of 2n/n], or [+ or -][square root of 2n/n], i.e., [+ or -] one or two times the standard deviation of [chi square](n).

We note that for model (1.5), this stopping rule is equivalent to that presented in [21].

2.2.2. Generalized cross validation. In [16, 18] it is shown that the generalized cross validation (GCV) method [5, 23] can be used to develop iteration dependent stopping rules for the steepest descent and RL iterations. Following their approach, in [1] a GCV method is developed for (2.3).

The GCV function for (2.1) at iteration k is defined by

(2.6) GCV(k) = n [[parallel][C.sup.-1/2]([A.sub.k] - (z - [gamma]))[parallel].sup.2]/trace[([I.sub.n] - [C.sup.- 1/2][AA.sub.k]).sup.2],

where [A.sub.k] is the iterative regularization matrix satisfying

[A.sub.k] (z - [gamma]) = [u.sub.k],

with [u.sub.k] the kth WMRNSD iterate. The idea is then to stop WMRNSD at the iteration k that results in an increase in the GCV function.

In order to evaluate GCV(k) in practice, trace([I.sub.n] - [C.sup.-1/2][AA.sub.k]) must be approximated. This can be accomplished using the Trace Lemma [23]: given B [member of] [R.sup.nxn] and v a discrete white noise vector, E([v.sup.T]Bv) = trace(B). Thus given a realization v from a white noise process,

trace([I.sub.n] - [C.sup.-1/2][AA.sub.k]) [approximately equal to] [v.sup.T]v - [v.sup.T][C.sjp.-1/2][AA.sub.k]v.

Thus, if we know [A.sub.k], we can efficiently estimate GCV(k) at each iteration. However, we don't know [A.sub.k], nor do we want to compute it. Instead, we follow [16] and define

[w.sub.k] = [A.sub.k]v,

which yields the following approximation of (2.6):

(2.7) [??](k) = n [[parallel][C.sup.-1/2][([Au.sub.k] - (z - [gamma]))[parallel].sup.2]/([v.sup.T]v - [v.sup.T] [C.sup.-1/2[Aw.sub.k]).

An iteration for [w.sub.k] can be derived, following [1], and is given by

(2.8) [w.sub.k+1] = [w.sub.k] - [[tau].sub.k][[w.sub.k] [dot encircle] [A.sup.T][C.sup.-1]([Au.sub.k] - (z - [gamma])) + [u.sub.k] [dot encircle] [A.sup.T][C.sup.-1] ([Aw.sub.k] - v)],

where [w.sub.0] = ([u.sub.0] [??] v)/(z - [gamma]).

Thus, we immediately have a stopping rule for WMRNSD iterations, namely, iterate (2.3) and (2.8) simultaneously and stop the iterations if [??](k) > [??] (k - 1), where [??](k) is defined in (2.7). As it is standard [23], we choose v so that its components are either -1 or 1 with equal probability.

[FIGURE 3.1 OMITTED]

2.2.3. The unbiased predictive risk estimator. Similar to GCV is the unbiased predictive risk estimator (UPRE) [23]. With UPRE, the goal is to choose the iteration k for which the predictive risk E (2/n[T.sub.wls]([u.sub.k]; [z.sub.e])), where [z.sub.e] = [Au.sub.e] + [gamma], is smallest. The UPRE function for (2.1) is

(2.9) UPRE(k) = 2/n[T.sub.wls]([u.sub.k]; z) + 2/n trace ([c.sup.-1/2][AA.sub.k])-1.

It is derived as in [23] for regular least-squares problems and is an unbiased estimator of the predictive risk, hence its name. Following the same approach as for GCV, we can approximate (2.9) by

(2.10) [??](k) = 2/n[T.sub.wls]([u.sub.k]; z) + 2/n [v.sup.T][C.sup.-1/2][Aw.sub.k] - 1,

where v and [w.sub.k] are as above.

Thus, we have our third stopping rule for WMRNSD iterations; namely, iterate (2.3) and (2.8) simultaneously, and stop the iterations if [??](k) > [??](k - 1), where UPRE(k) is defined in (2.10).

3. Numerical results. We now demonstrate the effectiveness of our methods on synthetically generated examples for each of the CT, PET, and SPECT imaging problems. For each of our tests we use the Shepp-Vardi phantom [20], given on the left in Figure 3.1, as the true image, generated by using MATLAB's phantom function. For data generation in CT, PET, and SPECT, we use the discretization of models (1.1), (1.3), and (1.2), respectively, described in Section 1.2. For PET, we ignore attenuation corresponding to [mu] = 0 in (1.3), which is common in synthetic numerical experiments done in the PET literature. For SPECT, we follow [4] and take [mu] = 1 where the absorption density of the body is positive and [mu] = 0 otherwise. The statistical model used for generating the CT data is (1.5), while for PET and SPECT it is (1.7). The noisy sinogram data for CT and SPECT are given in the middle and on the left in Figure 3.1. The sinogram for PET looks structurally very similar to the CT sinogram with noise more characteristic of the SPECT sinogram, thus we do not include it. The signal-to-noise level for these examples is approximately 30.

The inverse problem is to reconstruct the Shepp-Vardi phantom given the respective sinogram data z and the forward model matrix A. We do this by using the iterative method and stopping rule presented in Section 2 with an initial guess of [u.sub.0] = 1 and the three stopping rules, with [[epsilon].sub.n] =0 in (2.5). The computations were done in MATLAB.

[FIGURE 3.2 OMITTED]

[FIGURE 3.3 OMITTED]

For the CT example, the reconstruction given by the DP and UPRE stopping rules are given on the left and middle, respectively, in Figure 3.2. On the right is a plot of the relative error [parallel][u.sub.k] - [u.sub.true][parallel]/[parallel][u.sub.true][parallel] with the stopping iterations for DP, GCV, and UPRE labeled. In this example, the reconstructions for the three methods are visually very similar.

Before continuing, we note that GCV and UPRE depend, through iteration (2.8), upon the particular randomly drawn vector v. We have observed that different v yield different stopping iterations. Though the methods seem to be relatively stable with respect to the choice of v, additional stability can be obtained through averaging if multiple values of v are chosen, and multiple iterations of the form (2.8) are implemented.

Additionally, we have found that regularly restarting the iteration (2.8) within the WMRNSD iterations yields more effective stopping rules. In the results shown in Figure 3.2, as well as those that follow, we restart (2.8) every tenth WMRNSD iterations by setting [w.sub.k] = ([u.sub.k] [dot encircle] v)/(z - [gamma]).

For the PET example, the reconstruction given by the DP and GCV stopping rules are given on the left and middle in Figure 3.3. On the right is a plot of the relative error with the stopping iterations for DP, GCV, and UPRE labeled. In this case, the GCV and UPRE reconstructions appear slightly sharper than the DP reconstruction and have a slightly smaller relative error.

For the SPECT example, the reconstruction given by the DP and GCV stopping rules are given on the left and center in Figure 3.4. On the right is a plot of the relative error with the stopping iterations for DP, GCV, and UPRE labeled. In this case, once again, the GCV and UPRE reconstructions appear slightly sharper than the DP reconstruction, and have a slightly smaller relative error.

[FIGURE 3.4 OMITTED]

As it is the case for the discrepancy principle for regularization parameter selection in the case of Tikhonov regularization, the discrepancy principle consistently underestimates the stopping iteration. On the other hand, whereever we've chosen [[epsilon].sub.n] = 0 in (2.5), a negative value could reasonably be chosen as well, e.g., one standard deviation below the mean of the [chi square](n), or [[epsilon].sub.n] = [square root of 2n/n].

4. Conclusions. We have presented an iterative method and three statistically motivated stopping rules for solving inverse problems arising in the applications of CT, PET, and SPECT imaging. For PET and SPECT, a quadratic approximation of the negative-log Poisson likelihood function is needed in order to motivate the use of (2.1), whereas in the CT example Gaussian data, and hence a least-squares likelihood, are assumed. This allows for the implementation of the simple nonnegatively constrained WMRNSD algorithm of [1, 3] to be applied in all cases.

Stopping rules for WMRNSD follow from a Gaussian approximation of the data-noise model. For CT, such a model is assumed, whereas for PET and SPECT, an approximation (2.4) of a Poisson statistical model of the form (1.7) is used. The Gaussian/least-squares assumptions allow for the development of three stopping rules: the discrepancy principle, generalized cross validation, and unbiased predictive risk estimation. All three of these methods work in practice, as the numerical experiments show.

The codes used to perform the numerical experiments in this paper can be obtained from the author.

Acknowledgment. The author would like to thank Daniela Calvetti and Erkki Somersalo for their discussions of the PET problem and the National Science Foundation under grant DMS-0915107, for their support of this research.

REFERENCES

[1] J. M. BARDSLEY, Stopping rules for a nonnegatively constrained iterative method for ill-posed Poisson imaging problems, BIT, 48 (2008), pp. 651-664.

[2] J. M. BARDSLEY AND J. GOLDES, Regularization parameter selection methods for ill-posed Poisson maximum likelihood estimation, Inverse Problems, 25 (2009), 095005 (18 pages).

[3] J. M. BARDSLEY AND J. G. NAGY, Covariance-preconditioned iterative methods for nonnegatively constrained astronomical imaging, SIAM J. Matrix Anal. Appl., 27 (2006), pp. 1184-1198.

[4] P. GREEN, Bayesian reconstructions from emission tomography data using a modified EM algorithm, IEEE Trans. Med. Imag., 9 (1990), pp. 84-93.

[5] P. C. HANSEN, Rank-Deficient and Discrete Ill-Posed Problems, SIAM, Philadelphia, 1997.

[6] J. KAIPIO AND E. SOMERSALO, Statistical and Computational Inverse Problems, Springer, New York, 2005.

[7] L. KAUFMAN, Maximum likelihood, least squares, and penalized least squares for PET, IEEE Trans. Med. Imag., 12 (1993), pp. 200-214.

[8] L. KAUFMAN and A. NEUMAIER, PET regularization by envelope guided conjugate gradients, IEEE Trans. Med. Imag., 15 (1996), pp. 385-389.

[9] J. LLACER AND E. VEKLEROV, Feasible images and practical stopping rules for iterative algorithms in emission tomography, IEEE Trans. Med. Imag., 8 (1989), pp. 186-193.

[10] L. B. LUCY, An iterative technique for the rectification of observed distributions, Astron. J., 79 (1974), pp. 745-754.

[11] V. A. MOROZOV, On the solution of functional equations by the method of regularization, Soviet Math. Dokl., 7 (1966), pp. 414-417.

[12] J. NAGY AND Z. STRAKOS, Enforcing nonnegativity in image reconstruction algorithms, in Mathematical Modeling, Estimation, and Imaging, D. C. Wilson et al., eds., Proc. of SPIE, 4121, SPIE, 2000, pp. 182-190.

[13] F. NATTERER, The Mathematics of Computerized Tomography, SIAM, Philadelphia, 2001.

[14] F. NATTERER AND F. WUBBELING, Mathematical Methods in Image Reconstruction, SIAM, Philadelphia, 2001.

[15] J. M. OLLINGER AND J. A. FESSLER, Positron emission tomography, IEEE Signal Proc. Mag., 14 (1997), pp. 43-55.

[16] K. M. PERRY and S. J. REEVES, Generalized cross-validation as a stopping rule for the Richardson-Lucy algorithm, in The Restoration of HST Images and Spectra II., R. J. Hanisch and R. L. White, eds., Space Telescope Science Institute, Baltimore, 1994, pp. 97-103.

[17] S. PRASAD, Statistical-information-based performance criteria for Richardson-Lucy image deblurring, J. Opt. Soc. Amer. A, 19 (2002), pp. 1286-1296.

[18] S. J. REEVES, Generalized cross-validation as a stopping rule for the Richardson-Lucy algorithm, Int. J. Imag. Syst. Tech., 6 (1995), pp. 387-391.

[19] W. H. RICHARDSON, Bayesian-based iterative method of image restoration, J. Opt. Soc. Amer., 62 (1972), pp. 55-59.

[20] L. A. SHEPP AND Y. VARDI, Maximum likelihood reconstruction for emission tomography, IEEE Trans. Med. Imag., 1 (1982), pp. 113-122.

[21] H. J. TRUSSELL, Convergence criteria of iterative restoration methods, IEEE Trans. Acoust., Speech, Signal Process., 31 (1983), pp. 129-136.

[22] Y. VARDI, L. A. SHEPP, and L. KAUFMAN, A statistical model for positron emission tomography, J. Amer. Statist. Assoc., 89 (1985), pp. 8-37.

[23] C. R. VOGEL, Computational Methods for Inverse Problems, SIAM, Philadelphia, 2002.

JOHNATHAN M. BARDSLEY ([dagger])

* Received August 25, 2010. Accepted for publication October 18, 2010. Published online February 21, 2011. Recommended by R. Plemmons. This work was supported by the NSF under grant DMS-0915107, and by the Physics Department at the University of Otago, New Zealand, who housed the author during his 2010-11 sabbatical year.

([dagger]) Department of Mathematical Sciences, University of Montana, Missoula, MT, 59812 (johnathan.bardsley@umontana.edu).
COPYRIGHT 2011 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.