A Primal-Dual Method for Large-Scale Image Reconstruction in Emission Tomography.
In emission tomography, images can be reconstructed from a set of measured projections using a maximum likelihood (ML) criterion. In this paper, we present a primal-dual algorithm for large-scale three-dimensional image reconstruction. The primal-dual method is specialized to the ML reconstruction problem. The reconstruction problem is extremely large; in several of our data sets the Hessian of the objective function is the product of a 1.4 million by 63 million matrix and its scaled transpose. As such, we consider only approaches that are suitable for large-scale parallel computation. We apply a stabilization technique to the system of equations for computing the primal direction and demonstrate the need for stabilization when approximately solving the system using an early-terminated conjugate gradient iteration.We demonstrate that the primal-dual method for this problem converges faster than the logarithmic barrier method and considerably faster than the expectation maximization algorithm. The use of extrapolation in conjunction with the primal-dual method further reduces the overall computation required to achieve convergence.
Introduction
. In this paper we consider the image reconstruction problem
in emission tomography. This problem is encountered in the eld of nuclear medicine,
which is concerned with the study of organ function through radioactively labeled
\tracer" compounds. The quantity of interest in this problem is the spatial concentration
of radioactive emissions within the object under study. The quality of the
reconstructed image can depend upon a number of factors including the number of
emission events (i.e., counts) collected by the scanner and the method used to reconstruct
the image. In studies that are characterized by poor counting statistics (that
is, few counts), statistical reconstruction methods that model the Poisson nature
of the emission process have been shown to improve image quality over traditional,
non-statistical reconstruction methods [35, 57]. The low-count problem has generated
considerable interest in the medical imaging community because low radiotracer doses
and short scanning durations are highly desirable.
The estimation of emission density in an organ is an inherently three-dimensional
(3-D) process. Volume, or 3-D acquisition improves the counting statistics compared
with 2-D acquisition (in which axially oblique coincidences are either physically or
electronically blocked from detection) but increases the problem size considerably.
Since the 3-D problem may involve image and measurement vectors with millions
of elements, the amount of computation required to perform 3-D statistical reconstructions
can be quite substantial. In our computational studies for example, the
larger reconstructions consist of 1.4 million image variables which are reconstructed
from a measurement vector with 63 million elements. As such, it is important to use
reconstruction methods that converge rapidly. The statistical image reconstruction
Ariela Sofer is partly supported by National Science Foundation grants DMI-9414355 and DMI
9800544.
y Center for Information Technology, National Institutes of Health, Bethesda, Maryland, 20892-
5624 (johnson@mail.nih.gov).
z Department of Systems Engineering and Operations Research, George Mason University, Fairfax,
Virginia 22030-4444 (asofer@gmu.edu).
C.A. JOHNSON AND A. SOFER
problem can be posed as a constrained nonlinear optimization problem. In this paper
we present a primal-dual method for performing statistical 3-D reconstructions in
emission tomography that has been specialized to the intricacies of the application.
We demonstrate the rapid convergence of our primal-dual method in computational
studies on low-count, 3-D positron emission tomography (PET) data.
This paper is organized as follows. In Section 2 we present the statistical model
and develop the objective function. Section 3 reviews the EM method for ML recon-
struction. In Section 4 we develop a primal-dual method for ML reconstruction and
discuss initialization, stabilization, and extrapolation enhancements. Computational
tests comparing the primal-dual results to a logarithmic barrier approach and the EM
method on small animal data are presented in Section 5. Some concluding remarks
are made in Section 6.
2. Statistical model and objective function. We begin our discussion by
forming a nite parameter space for the image estimates, as is customary [20]. Consider
the situation depicted in Figure 2.1 where a grid of boxes or voxels has been
imposed over the emitting object (for simplicity, the Figure is depicted in 2-D; the
concept is readily extended to 3-D). Given a set of measurements along lines of coin-
cidence, we seek to estimate x expected number of counts
emitted from voxel i. Let X i be the number of radioactive events emitted from voxel
are assumed to be independent Poisson-distributed random variables with mean
system matrix C 2 < nN is used to model a number of physical eects
including spatially dependent resolution and attenuation. The elements C i;j of the
system matrix represent the probability that an event emitted from voxel i will be
detected by detector pair (coincidence line). The number of events emitted from voxel
and detected at coincidence line j is therefore  are also independent
Poisson variables. The measurements y j are thus realizations of sums of
independent Poisson variables y
The above is a considerably simplied model of the actual measurement process; for
further discussion on its validity to the present situation, see [24].
Given our simplied Poisson model, the likelihood may be written as
Y
Y
e
The ML objective function is formed by taking the log likelihood
log
Ignoring the constant term, we dene our objective function fML (x) as
is a vector of 1's, so that q is the sum of the
columns of C (which need not necessarily be 1). Dening
Fig. 2.1. Relationship between estimate x i and measurement y j . Shown here is the case of
PET, where emission-count measurements are taken along coincidence lines from pairs of detectors.
A nite parameter space is formed by imposing a grid of voxels over the emitting region. The
estimate of the expected emission intensity within voxel i is x i .
to be a forward transformation, we can write the gradient and Hessian of the objective
function, respectively, as
The Hessian is negative semidenite (since
so the objective function (2.1) is concave. Thus, any local maximum
will also be a global maximum.
Equation (2.2) sheds some insight into the computational costs associated with
maximizing the objective function. Given a current solution estimate x k ; computing
the gradient requires rst computing a forward transformation ^
then computing a backward transformation
k y from the forward transformation.
The costs of performing the forward transformation and backward transformation are
similar and together dominate the computation associated with iterative reconstruction
methods, especially in large scale. We shall revisit this computational structure,
which is common to all iterative reconstruction methods.
Since the underlying activity distribution is non-negative, the ML reconstruction
problem is a constrained optimization problem with lower-bound constraints:
maximize fML (x)
subject to x  0:
The ML objective function has a nite maximum and compact level sets on x  0
[36].
2.1. Maximum a posteriori reconstruction. Without regularity conditions
on x, estimating the spatial emission distribution is a statistically ill-posed problem
[7, 33]. The fully converged ML reconstruction, being dominated by noise and edge
artifact, is not generally of biomedical interest [55]. Regularization can be included
4 C.A. JOHNSON AND A. SOFER
in the objective function by introducing a Bayesian formulation [20, 37]. Given prior
probabilities P fxg and P fyg for the image and measurements, respectively, we dene
the posterior probability
The estimate of x is then obtained by maximizing the posterior probability P fxjyg.
A common choice for the image prior is the Gibbs distribution P
although other priors (e.g., Gaussian, Gamma) have been investigated [37, 39]. The
popularity of Gibbs priors stems in part from their ability to capture the local correlation
property of images [19]. The energy function R is dened as a sum of "potential"
functions designed to discourage non-smoothness in a neighborhood
denotes the neighborhood of voxel i. In order to maintain concavity and
twice continuous dierentiability in the objective function, the potential function V i;l
is chosen to be convex with continuous rst and second derivatives. In our studies we
have used the potential function V i;l
z
log
z
and - is a shaping constant that we typically set to 1 [38].
For maximum a posteriori (MAP) reconstructions, the objective function is the
log-posterior likelihood log P fxjyg. Ignoring a constant, our objective function become

The MAP reconstruction problem can also be posed as a constrained optimization
problem
maximize fMAP (x)
subject to x  0:
We note for future reference the following:
Although the function R (with the potential function (2.5)) is concave it is not strictly
concave. only for vectors v that are a scalar multiple of
the unit vector e N , and since e T
negative denite and that fMAP is strictly concave [38]. In addition, fMAP has a
nite maximum and bounded level sets on x  0 [37].
2.2. The optimization problem. For convenience of notation, let us pose the
reconstruction problem as a constrained minimization problem:
subject to x  0;
(2.
0. The case
corresponds to the unregularized
ML; in general we shall be more interested in the fully converged MAP where
The Karush-Kuhn-Tucker (KKT) rst-order necessary conditions for optimality
of (2.8) at a point x are existence of Lagrange multipliers  so that
x is the Lagrangian function. Due to the strict convexity of
f , the second-order su-ciency conditions are satised, and x is the unique minimizer
of f .
3. The EM algorithm. The expectation maximization (EM) method, as presented
by Dempster, Laird, and Rubin [8] for ML estimation, is an iterative algorithm
for computing ML estimates when the measurements are viewed as incomplete data.
Shepp and Vardi [53] and Lange and Carson [36] applied the EM method to emission
and transmission tomography problems, respectively. The EM algorithm has been
proven to converge to an optimal solution of (2.4) [36, 56].
The EM algorithm for emission tomography can be derived [56, 27] from the
optimality conditions for the reconstruction problem. For the unregularized problem
can be written as
diag (x), and Premultiplication by X , and utilizing the complementary
slackness condition yields
or since
Applying a xed-point algorithm x
to the above equation yields the
EM update equation
where x k is the current image estimate,
diag
. Given a positive initial solution x the algorithm maintains non-negativity
at every iteration and converges to a xed point x which is
an optimal solution of (2.4). The asymptotic rate of convergence is governed by the
spectral radius of rM which is typically very close to unity. In one example
using reasonable assumptions about the scanner geometry, the lower bound of the
spectral radius was calculated to be .99938 [17]. Indeed, EM has been observed to
converge very slowly, especially close to the optimal solution. The slow convergence of
the EM algorithm has limited its clinical applicability. The cost of one EM iteration
is equivalent to the cost of one gradient calculation.
In MAP-EM, the presence of the regularizing term in (2.6) precludes a closed-form
update equation such as (3.1) for ML-EM. We mention two algorithms that are commonly
used for MAP-EM reconstructions: the \one step late" (OSL) algorithm and
6 C.A. JOHNSON AND A. SOFER
DePierro's algorithm. Green's OSL algorithm approximates R (x) with the constant
, thereby permitting a closed-form approximated update [16, 17]
. OSL converges to the MAP solution provided that
, where
is an upper threshold for the prior strength. DePierro's algorithm is
a \true" MAP-EM implementation that substitutes the convex function R (x) with a
separable, convex, and twice continuously dierentiable function R x; x k
R (x), so
that separable maximizations can be performed on the variables [9, 10]. Regularization
improves the convergence rate of EM, with larger prior strengths resulting in lower
spectral radii. However, for reasonable prior strengths (mild to moderate smoothing),
the convergence rates of OSL and DePierro's algorithm are still quite close to unity.
The EM update formula on the right-hand side of (3.1) follows Kaufman [27],
who was the rst to pose the EM algorithm as an optimization algorithm (namely, a
scaled steepest-ascent method). This representation allows for the inclusion of a line
search [27, 28] to accelerate the method's performance. Likewise, the MAP update
(3.2) can be enhanced by a line search.
Several other approaches for solving the maximum likelihood estimation problem
have been proposed. These include preconditioned conjugate gradient techniques
[27, 28, 34, 42] or truncated-Newton methods [27, 28]. The nonnegativity constraints
are maintained either be limiting the step length or by using a bending line search.
The paper [44] explores active set methods, while [43] enforces nonnegativity via
a quadratic penalty in the objective. In other work [29, 30, 31] a penalized least-squares
objective is used instead of the maximum likelihood. These problems are
solved by a preconditioned conjugate gradient and use specialized techniques to drive
the complementary slackness to zero.
There is considerable debate within the PET community regarding the appropriate
model for reconstruction. It has long been observed that the unregularized
maximum likelihood estimator gives grainy images. However if the EM algorithm is
stopped early, the resulting solution often produces images of acceptable quality. For
this reason some researchers argue that early termination is a form of smoothing, and
that no regularization is needed. Proponents of MAP argue that the approach allows
the user to control the amount of regularization through the parameter, and that the
regularized objective function is better conditioned. In either case, it has been observed
that EM-type algorithms may lead to nonuniform convergence. In particular,
the algorithms may converge slowly in \cold spots" (regions of low activity within
regions of activity) and in areas of isolated activity within cold spots. The use of an
interior-point algorithm oers the hope of more uniform convergence.
4. A primal-dual approach. The drawbacks of the EM algorithm motivate
our investigation into interior-point approaches for the ML and MAP reconstruction
problems. As is clear from (2.1), the objective function can be undened outside the
feasible region x  0. Thus the ML and MAP reconstruction problems would appear
to be \natural" candidates for interior-point algorithms. The reconstruction problem
is especially suited to an interior-point approach, because its output is a gray-scale
image. Whether a particular value is exactly \zero" or just very close to zero is
immaterial. Slight inaccuracies below the gray scale threshold are inconsequential;
obtaining an image rapidly is a neccessity.
Primal-dual methods have enjoyed considerable success in linear programming
[18, 32, 40], and have recently been proposed for nonlinear programming [5, 13, 41].
Although they are closely related to the logarithmic barrier method, primal-dual
methods may pose some advantages. In the logarithmic barrier method, the Lagrange
multiplier estimates may be inaccurate when the primal variables are not close to the
barrier trajectory [11]. Primal-dual methods oer the potential of improved \center-
ing" over barrier methods. Given the size of the current problem, the developments
presented here must be suitable for large-scale parallel computation.
In a manner similar to classical barrier methods, primal-dual methods attempt
to follow the \barrier trajectory," a smooth trajectory characterized by a barrier parameter
[12]. The points along the trajectory satisfy a perturbed
version of the KKT conditions:
Dening ng and our method
maintains (4.3) while attempting to solve (4.1), (4.2), that is
Xen en
0:
Given the point x
and the barrier parameter , the search direction
prescribed by Newton's method satises the \unsymmetric" primal-dual
equations [41]:
I
Elimination of the (1,2) block of the matrix in (4.5) yields the reduced system
where the \condensed" primal-dual matrix is given by
We have implemented an algorithm in which the primal and dual variables are
permitted to take separate steplengths:
The primal steplength  x is chosen to ensure su-cient decrease in the merit function
log
Observe that F (x; ) is simply the logarithmic barrier function and that
8 C.A. JOHNSON AND A. SOFER
is identical to the right-hand side of (4.6) for . The unconstrained
minimizer x () of F (x; ) satises the perturbed KKT conditions (4.1)-(4.3) with
corresponding multiplier  i Furthermore, the solution
of the condensed primal-dual Newton equation (4.6) is guaranteed to be a descent
direction of the merit function for  > 0, since
and M is positive denite. We shall discuss in further detail the computation of the
primal search direction and step length.
The formula for the dual step length follows a suggestion by Conn, Gould, and
Toint (CGT) [5]. If  lies component-wise in the interval
(where  is a constant parameter that we have set to 100) then
otherwise nd 0 <   < 1 such that
subject to  k+1 being in the interval (4.9). These conditions on the dual step might
appear at rst glance to be overly restrictive but are actually designed to give maxi-
mum
exibility in the choice of  k+1 . CGT use these bounds on  and nonsingularity
of M to prove that, for any xed parameter value
, the minimization of F (x;
must be successful, that is, eventually a solution is found that satises the perturbed
KKT conditions (4.1)-(4.3).
In general it is neither necessary nor desirable to reach full subproblem conver-
gence. Rather, we have implemented a \short-step" algorithm in which only one
primal-dual step is usually needed before adjusting . Setting the barrier parameter
is an important consideration in primal-dual algorithms, and has a strong in
uence
on the convergence rate. A reduction in  k is performed whenever the \-criticality"
conditions [5, 54] are satised:

k+1
are constant parameters. If the above conditions are satised, the
barrier parameter is reduced according to
where  is a constant parameter such that
A consequence of (4.14) is that  k cannot increase. Furthermore, since the minimization
of F (x; ) must be successful, a -critical solution (a weaker requirement) must
eventually be found. Thus it is impossible for  k to be non-decreasing. Using this
argument, CGT prove that the algorithm must converge to a KKT solution [5].
In practice we nd that both the primal and dual direction vectors are well scaled,
and that  x and   are both typically close to 1. By far the most costly operations
are computing the primal direction p x and updating the gradient rF (x), as we shall
explore. In contrast, the costs of the line search for the primal steplength, the computation
of the dual search direction (4.7), and the dual line search (4.9) are relatively
insignicant. From empirical evidence in our computational studies, we have
found that a \short-step" algorithm with gradual reduction in  achieves the fastest
convergence to the KKT conditions. Specically, we dene #
These parameter values enable the -critical conditions to be met after
only one primal-dual step for most subproblems.
4.1. Computing the primal direction. For large problems, factoring the condensed
primal-dual matrix M or even forming the Hessian r 2 f (x) would be prohibitive
due to the size of the matrix (376,000376,000 for even the smaller reconstructions
being considered in this paper) and the enormous amount of computation
that would be required. Thus we must consider methods for approximating the Newton
direction in (4.6). The approach we have successfully applied to this problem is
motivated by the truncated-Newton [6] method of unconstrained optimization. The
search direction is an approximate or truncated solution to the Newton equations
[47, 49]
An early-terminated conjugate gradient (CG) iteration [21] is used to obtain an approximate
solution to (4.15).
An equivalent statement of (4.15) is we seek to nd the direction p x that approximately
minimizes the quadratic Q (p x
x reasonable and
eective truncation point for (4.15), based on the monotonicity of Q (p x ), is proposed
in [48]; the CG is terminated at subiteration l if
x
x
The CG termination rule (4.16) has been an important component of the reconstruction
software in that it consistently yields a well-scaled primal direction vector as long
as    s , where  s is a threshold value below which stabilization is required (we
shall discuss the  <  s case in Section 4.2).
The CG method does not require storage of the Hessian or condensed primal-dual
matrix, but rather only application of matrix-vector products. From (2.3) we
can write the rst term of the matrix-vector product
for an arbitrary vector v 2 < n . Computationally, (4.17) consists of a forward transformation
(C T v) followed by a diagonal scaling (^y is already available from the computation
of rf (x)), followed by a backward transformation (premultiplication by C).
To be explicit, recalling (4.8), we have
C.A. JOHNSON AND A. SOFER
where r 2 R (x) v can be computed exactly without incurring signicant computational
expense. The forward-and-back-transformation operation in (4.17) dominates the
computational cost of a CG iteration. This operation is computationally similar to
computing the gradient, or one EM iteration.
Some authors advocate solving simultaneously for p x and p  , using the full unsymmetric
primal-dual equations (4.5), or an equivalent symmetrized system [13, 14, 52,
61]. The unsymmetric primal-dual matrix in particular remains nonsingular, and its
condition number remains bounded as  ! 0 [12, 41], when the standard conditions
of a constraint qualication, strict complementarity, and the second-order su-cient
conditions are satised at the solution. In our application, due to the size of our
problem, we must use an iterative method. We believe that solving a symmetric system
via a symmetric solver such as the CG would be more e-cient than solving the
full unsymmetric system via an unsymmetric iterative solver such as GMRES (even
though our symmetric system is ill-conditioned), since the the amount of work and
storage required per iteration in GMRES increases linearly with the iteration count.
An advantage of using the condensed system (4.6){(4.7) is that although the primal
search direction is computed inexactly, the equation for maintaining complementarity
(4.7) is maintained. In practice we nd that the resulting primal and dual direction
vectors are both well scaled, and that  x and   are typically close to 1.
4.1.1. Preconditioning. The use of a preconditioner with the CG is essential
for a competitive algorithm. Since every CG subiteration is as costly as a gradient
evaluation or EM iteration, it is highly desirable to obtain a quality direction vector
in as few CG iterations per subproblem as possible. We have investigated a number of
preconditioners, including FFT-based preconditioners that model the approximately
Toeplitz-block-Toeplitz nature of CC T with a circulant-block-circulant approximation
[2, 3], high-pass lter approximations to the FFT-based preconditioner [4], the EM
preconditioner XQ 1 [34], the exact diagonal of M , and diagonal Hessian approximations
[46].
Of the above preconditioners, by far the best-performing was the exact diagonal
of M , which can be computed at reasonable cost:
Note that the rst right-hand side term in (4.18) is similar in form to a backward
transformation, although a bit more expensive due to the squaring operations. We
have found that the preconditioned CG method using an exact diagonal preconditioner
in the form of (4.18) almost always requires using fewer than 10 iterations to
achieve (4.16), regardless of the size of the problem. In many cases, only 3 or 4 CG
iterations are required. Moreover, the directions produced using an exact diagonal
preconditioner are well scaled (usually resulting in primal step sizes of near 1), and
lead to rapid descent.
In contrast, the other preconditioners did not perform well. Already in the initial
subproblems they tended to yield a poorly-scaled search direction, which in turn,
resulted in small steplengths. Subsequent calls to the CG suered further from this
problem, and the algorithm made little progress. This behavior was particularly
surprising for the block-circulant FFT-based preconditioners. These preconditioners
perform very well in other reconstruction methods, especially in least-squares methods
where the block-circulant approximation is well matched to the Hessian structure. We
were motivated to try them for our problem because the ML Hessian is almost block
circulant. But because of the strong diagonal component in M and its spatially-
variant dependence on y, ^
y, x, and , shift-invariant Toeplitz models of M yield a
poor approximation in our method.
4.1.2. Line search. For ML and MAP reconstructions, knowledge of the structure
of the objective function can lead to a substantial reduction in the cost of implementing
a line search over a more naive approach. Specically, after the search
direction p x has been found, and once a forward transformation ^
been computed, it is possible to compute the objective function and rst and second
directional derivative values at the trial points x k
at nearly negligible cost.
To see this, note that ^
and therefore [27, 28]
Similar expressions exist for the directional rst and second derivatives [24].
After the initial forward transformation to compute ^
no further forward- or
back-transformation operations are required during the line search at any of the trial
points. The forward transformation ^
w can be re-used, so that only one backward
transformation is subsequently required to update the gradient. The above observations
and the well behaved convex nature of the objective function have permitted
us to implement a highly accurate but low-cost Newton line search. Due to the low
cost of each step we have chosen a relatively strict tolerance of 0:05 on the Wolfe
condition for termination of the line search: We nd this line search technique to be
highly eective and, in no small part, responsible for the positive results we report.
4.2. Stabilization. A well known property of the Hessian of the primal barrier
function is its increasingly ill-conditioned nature as  ! 0 [45]. Analogous results
hold for the condensed primal-dual matrix: as the solution is approached the matrix
becomes increasingly ill-conditioned. (For a detailed analysis see the paper by Wright
[60]).
In [50], Nash and Sofer developed an approximation to the Newton direction
for the logarithmic barrier, that avoids the structural ill-conditioning of the barrier
Hessian and is suitable for large-scale problems. The direction is the sum of two
vectors, one in the null space of the Jacobian of the active constraints, and the other
orthogonal to it. The associated decoupling is based on a prediction of the binding
set at the solution.
We have recently adapted this approximation to the condensed Newton equations
arising in primal-dual methods. Although our derivation is valid for general nonlinear
constraints, we present it here for the special case of bound constraints in the context
of (4.6).
We will assume in the following that strict complementarity holds at the solution,
that is,
0g to be the index set of binding
constraints at the solution, and ^
n to be the number of binding constraints at the
solution. We will assume that 0 <
always the case in reconstructions of
practical interest. Dene I
0g the set of nonbinding constraints. Let x I
be the subvector of variables that are positive at the optimal solution, and x J the
subvector of variables that are zero at the optimal solution. Assume also that the
12 C.A. JOHNSON AND A. SOFER
variables are ordered so that the positive variables are rst, i.e.,
x I
x J
The Hessian of the objective function will then be similarly partitioned,
as will the condensed primal-dual matrix
I;J MJ ;J
I  I H I;J
where X I , XJ ,  I and J are the diagonal matrices of the associated components
of x and .
We will assume that the sequence of iterates (x; ) generated by the primal-dual
satises the following properties, when  is su-cently small:
Here we dene there exist constants 0 <  l <  u so that  l   kk   u
for all su-ciently small  > 0. We say that a vector or matrix is () if its norm
is (). We also dene  = O() if there exists some positive constant  u so that
kk   u  for all su-ciently small  > 0.
We will also assume that near the solution the Hessian is reasonably well condi-
tioned, so that Now the diagonal terms of MJ ;J are O(1=), and become
unbounded as  ! 0. In contrast, the diagonal terms of M I;I dier from those of the
reduced Hessian H I;I by O(), and the condition of M I;I thus re
ects that of the constrained
problem. The condensed primal-dual matrix M can then be shown to have
\large" eigenvalues of magnitude (1=), and
\small" eigenvalues that dier
from those of H I;I by O(), and have magnitude (1). The condensed primal-dual
matrix thus suers from the same structured ill-conditioning as the barrier Hessian.
For small values of  we propose approximating the primal Newton direction p x ,
by a direction ~
whose null- and range-space components are computed as follows:
~
The system for computing the component ~
I
x involves the well conditioned matrix
I;I , and can be solved exactly or inexactly via the conjugate gradient method. The
computation of ~
x is straightforward. Thus, the ill-conditioning of the condensed
primal-dual is avoided. We will show now that under the assumptions above, ~
so that the accuracy of the approximation increases as the solution is
approached and the potential harm from ill-conditioning increases.
Using the well known formula for the inverse of a partioned matrix (see e.g.
[51, 61]) it follows that
~
I
I;I rF I +M I;J (XJ  1
I;I rF I (XJ  1
where
I;I
Now by denition
so that G
Note further, that
I
whereas
It follows that
~
I
and
~
so that ~
In [50], Nash and Sofer prove (for the case of the Newton direction arising from
the logarithmic barrier objective function) that, for su-ciently small ; the vector
computed using an approximation similar to (4.19) and (4.20) yields a descent direction
with respect to the logarithmic barrier objective function. The proof is readily
extended to the present primal-dual case; thus p x is a descent direction for the merit
function F (x; ). We have found that, for the present problem, the above approximation
to the Newton direction is useful for values of  of order 10 4 or less.
Recently Wright [60] showed that the errors generated by backward-stable numerical
methods (various Cholesky factorizations and Gaussian elimination with partial
pivoting) for solving (4.6) are not magnied by the structured ill-conditioning. These
methods are inappropriate for our large problems which involve potentially millions of
variables. Instead we nd an approximate solution using a CG iteration. When working
in inexact arithmetic with large numbers of variables, the convergence rate of the
CG method depends on the condition of M [15]. Thus the structural ill-conditioning
in M can lead the CG iteration to spend an unnecessary amount of work in computing
Further, as we have observed, the criterion for terminating the CG may be overly
optimistic in an ill-conditioned system, so that the resulting direction is poorly scaled
as
The potential eect of ill-conditioning is illustrated through an example in Table
4.1. This example was encountered during development and motivated the incorporation
of stabilization into the algorithm. Starting at the subproblem
the primal steplength, dual steplength, and ncg (the number of CG it-
erations), are listed for both the non-stabilized and stabilized cases. This test was
terminated at  T x=n  7:5  10 5 . Note that in the non-stabilized case, the number
of CG iterations from the rst subproblem in the test to termination is signicantly
lower in the stabilized test than the non-stabilized test. Note also that in many of the
non-stabilized subproblems, either the primal or dual steplength is small, indicating
a poorly scaled direction or loss of accuracy.
14 C.A. JOHNSON AND A. SOFER

Table
An example of the eect of stabilization. The number of CG iterations, ncg, is counted from
the beginning of the subproblem. The termination condition in this example is
non-stabilized stabilized
7.08E-5 0.392 1.000 46 3.25E-5
5.83E-5 1.000 0.166 51
4.77E-5 1.000 62There has been much recent interest in stabilization methods that do not require
a prediction of the active set [13, 14, 59]. These approaches are based on factorization
methods which are unsuitable for a problem as large as the present one. The
argument against stabilization methods that require a prediction set is that the active
set is unknown in interior-point methods. We argue that, close to the solution
in the emission tomography reconstruction problem, an accurate prediction of the
active set can be made. In our problem, the constraints have a simple interpretation.
The positive variables correspond to those voxels containing at least some radioactive
tracer, while the zero-valued variables correspond to those voxels that lack any tracer
activity. Close to the solution, when  becomes su-ciently small that stabilization
is appropriate, the set of binding constraints is obvious and can be conservatively
identied with a -dependent threshold.
4.3. Extrapolation. Fiacco and McCormick showed that the solutions x () at
the perturbed KKT solutions form a unique dierentiable trajectory in  [12]. The
perturbed KKT conditions (4.1){(4.3) dene a \central path" as  ! 0. Thus, a
successful algorithm may be able to move both \along" and \toward" the path. As
discussed in [12], from the subproblem solutions fx ( l the trajectory
can be approximated as a polynomial
l=k r
c l  l ;
where r is the degree of the approximating polynomial and c k r are r
vectors of coe-cients. Using the approximation in (4.21), we nd a direction x such
that
l=k r
c l  l x
and set
to be a prediction to the next subproblem's primal solution. Here x k is the computed
(approximate) subproblem solution for Primal feasibility is maintained by
the steplength
is the maximum steplength that does not
violate non-negativity in x. Then, in the manner of (4.7), we compute a dual direction
vector according to
The dual vector is then moved according to
4x;
which requires another dual line search to minimize (4.10). The resulting point
serves as a starting point for the 1)st subproblem, a prediction
to the solution at  k+1 . The extrapolated primal-dual method can be viewed as
a predictor-corrector algorithm, with the extrapolation (4.22 and 4.24) serving as the
\predictor" step, and the subproblem minimization serving as the centering or \cor-
rector" step [23]. The degree r of the approximating polynomial is 1 when predicting
the 3rd subproblem, 2 for the 4th, and 3 for the 5th and beyond.
We have experimented with line searches in conjunction with (4.22), but often
1, and hence the line search just yields
. For this reason, we have found that
(4.24) yields a more eective dual direction than does the equivalent of (4.7) in the
context of extrapolation. Although the extrapolated search direction x can often
be poorly scaled (i.e.,
1), we have observed that the directions produced are
always descent directions to the merit function and lead to a signicant decrease in
the objective function f . A number of reconstructions were performed in which
was computed by extrapolating the dual solution vector (rather than computing it via
(4.24)); the discouraging nature of the results led us to abandon direct extrapolation
of the dual vector in favor of (4.24) which is highly eective in comparison.
Following extrapolation, a gradient evaluation is required to update the vector
the primal-dual algorithm requires between 12 and 25
subproblems to perform a 3-D MAP reconstruction, extrapolation adds that many
gradient evaluation operations to the computational cost. So extrapolation is only
economical if it reduces the computational burden by at least as much as it adds.
Our experience has been that for some data sets, the cost of extrapolation is well
worthwhile but for other data sets the benets were only marginal. Extrapolation
thus appears to serve as somewhat of a safeguard against di-cult problems. In an
extrapolated primal-dual reconstruction, the convergence measure does not
decrease as monotonically as in a primal dual reconstruction without extrapolation.
Certain extrapolated steps seem to cause the algorithm to \get ahead of itself," but
this eect is transient. On the studies we've performed, the algorithm does ultimately
converge to an accurate solution with extrapolation.
4.4. Initialization. The choice of the initial barrier parameter may have a substantial
eect on the algorithm. If the parameter is too small, the rst subproblem
may have extreme di-culty due to ill conditioning; if the parameter is too large,
then many (unnecessary) subproblems will be required to solve the problem. Proper
initialization of the barrier parameter  involves nding the most suitable point on
the barrier trajectory based on the initial solution x and the measurement data y.
Recalling the perturbed necessary conditions in (4.1), if the initial solution were to
be on the central path, it would satisfy
C.A. JOHNSON AND A. SOFER
Pre-multiplying by ^
T we arrive at
r
This suggests the following rule for initialization, which we nd quite eective:
r
Another, similar, initialization rule is motivated by the goal of nding an initial
value  0 so that
While (4.26) cannot be solved exactly, we can try to nd a  0 that results in a point ^
that is close to the barrier trajectory according to, say, the 2-norm. This motivation
leads to an alternative initialization rule [51]


During the course of development, both initialization rules were tried on certain data
sets. Although both initialization rules performed well, reconstructions initialized
with (4.25) usually reached the optimal solution in slightly less overall work than
those initialized with (4.27).
The initial estimate for ^
used most frequently was in each case a
positive uniform eld. A discussion on the rationale of using a uniform eld for
and on criteria for choosing the constant value of the primal initial solution may be
found in [24]. Alternative choices for the initial dual vector may be preferable, and
an investigation into this question may be worthwhile.
4.5. Termination. Given that subproblem termination is based on the -criti-
cality conditions (4.11) and (4.12), the closeness of each subproblem solution can be
measured by . If subproblems are solved exactly, jf [12]. The
-criticality conditions, however, are designed for a \short-step" algorithm in which
one truncated-Newton step should satisfy each subproblem for su-ciently small . To
ensure the accuracy of the nal solution, nal termination is based on the following
two requirements:
We have found that reasonably accurate solutions are ensured when "
The traditional view in tomographic reconstruction is that a highly accurate solution
is unnecessary. This view stems in part from the ill-posedness of the problem
and the computational cost of taking a reconstruction to full convergence. From empirical
evidence in our studies, the ability to perform certain imaging tasks such as
\cold spot detectability" improves with accuracy of the solution. Although the termination
criteria we propose above may not appear particularly strict, they are from
a tomographic reconstruction perspective.

Table
Properties aecting computation, memory, and storage costs for two dierent-sized reconstruction
problems. Gradient evaluation costs are based on a 2.5M-count study on 10 120-MHz IBM
RISC/6000 SP processors.
size class n N elements density storage cost of
in C in C cost of C gradient
thick-slice 376,882 5:36
thin-slice
5. Computational studies. To test our algorithm we have performed a number
of reconstructions on data acquired from a small animal scanner, and on data
generated by Monte Carlo simulations on the same animal scanner.
5.1. Size of the problem. Our studies involved two dierent-sized problems.
Raw coincidence data from the scanner can be binned into either \thick-slice" or
\thin-slice" measurement spaces, or both. \Thick-slice" reconstructions, in which
minutes for a gradient evaluation
using processors (120 MHz) on a 2.5M-count study. For
a \thin-slice" reconstruction with on the same
data and processors, a gradient evaluation requires 6.75 minutes. These properties
are summarized in Table 5.1. The cost of storing the full n  N system matrix is
prohibitive, even for thick-slice reconstructions. Extensive exploitation of the sparsity
and symmetries inherent in the system matrix makes its storage and retrieval possible
[24, 25].
The dominant computational operations of the reconstruction problems are the
forward- and back-transformation operations that underlie EM iterations, gradient
evaluations, Hessian-vector products, and diagonal Hessian calculations. These operations
have been implemented in parallel via a data decomposition strategy that partitions
the \measurement-space" vectors y and ^
y across the processors. The \image-
space" vectors such as x and  are replicated over all processors. Our data decomposition
is justiable under the observation that N >> n. On a data set with 2.5M
counts, at most 47% of the elements of y will be nonzero in the thick-slice case; at
most 4% in the thin-slice case. (The thin-slice conguration has over 10 times as
many lines of response as the thick-slice.) The dominant computational operations
have been implemented in such a way to exploit sparsity in y and further conserve
computation [24].
5.2. Cost metrics. We have devised metrics to measure the cost of an interior
point reconstruction. Dene the number of subproblems to be npr, the number of
truncated-Newton iterations nit, the number of conjugate gradient subiterations ncg.
The cost of one CG iteration (dominated by the Hessian-vector product) is equivalent
to the cost of one gradient calculation or EM iteration. One truncated-Newton iteration
requires, in addition to the ncg operations, one diagonal Hessian evaluation plus
one forward transformation and one backward transformation. The exact cost of these
operations varies depending on the size of the problem and number of counts, but we
shall approximate the cost of one truncated-Newton iteration to be the equivalent of
two gradient calculations beyond the cost of the conjugate gradients.
Using this approximation, the total cost of unextrapolated interior-point reconstructions
can be measured in units of equivalent number of gradient calculations (or
C.A. JOHNSON AND A. SOFER

Table
Summary of thick-slice primal-dual results and comparison with MAP-EM and LSEM. Extrapolation
was not used, and in all cases 2.
study f  npr nit ncg ngr MAP-EM LSEM
A 2,465,770 19 19 110 148 1000 344
G 3,660,344 24 24 127 175 > 1000 724
average ngr 183

Table
Summary of thick-slice extrapolated primal-dual results and comparison with MAP-EM and
LSEM; in all cases 2.
study f  npr nit ncg ngr MAP-EM LSEM
A 2,465,772 17 17 94 145 960 332
F 3,296,029
G 3,660,384 20 20 100 160 >1000 430
average ngr 156
EM iterations):
Extrapolation requires an additional gradient calculation following the extrapolation
in order to update the gradient vector. With extrapolation we modify the formula to
5.3. Computational results. We have performed a number of 3-D reconstructions
on data acquired from a small animal scanner and data generated by a Monte
Carlo simulation of the same small animal scanner. Reconstructions of seven datasets
were taken to full convergence, as dened by the termination criteria (4.28) and (4.29)
with . The various datasets used in our computational
studies represent a fairly diverse sample of the types of scans that might be
encountered in practice. The number of counts in the datasets used in these studies
ranged from 850K to 5.1M. The number of binding constraints at the optimal solution
ranged from approximately 20% to 80%.
Our main results are summarized in Tables 5.2 and 5.3 for the non-extrapolated
and extrapolated primal-dual cases, respectively. Studies A through D are reconstructions
of data acquired from a small animal PET scanner, while studies E through G
are reconstructions of Monte Carlo simulated data. These reconstructions were performed
in \thick-slice" mode (376,832 variables) with the regularization parameter
set at
In these tables, the column "MAP-EM" indicates the number
of DePierro MAP-EM iterations that were required to achieve the value of f  in the
same row. The column "LSEM" indicates the number of iterations required for an
EM algorithm, where the search direction on the last term of (3.2) is enhanced by a

Table
Summary of thick-slice logarithmic barrier results and comparison with MAP-EM and LSEM.
Extrapolation was used on all data sets, and in all cases
study f  npr nit ncg ngr MAP-EM LSEM
A 2,465,832 5 28 159 218 880 194
average ngr 265
line search. (To avoid excessive computation, the function values were only calculated
every iterations, and the nal count was rounded down, to favor this
method.) Since the cost of one gradient evaluation is equivalent to the cost of one EM
iteration, the numbers in the columns ngr and MAP-EM and LSEM can be compared
directly. We nd that the primal-dual method consistently reaches convergence much
more rapidly than either MAP-EM or LSEM.
Another interesting observation can be made in the comparison between Tables
5.2 and 5.3. Consider the number of EM iterations required to reach f  for study C.
In

Table

5.2, the LSEM algorithm reached iterations. In Table
5.3 on the same data set, the LSEM algorithm reached
tions. Thus, the algorithm took 64 iterations to reduce the function value by only 10
units near the solution. MAP-EM did even worse, requiring 180 iterations to reduce
the function value by 10. This is in fact a typical example of the slow limit behavior
of the EM algorithm. In all studies, the EM method did not achieve the same convergence
results obtained by the primal-dual method at termination. The Lagrangian
gradient norm and complementary slackness values of the terminated MAP-EM and
LSEM iterates were consistently much higher than those of the terminated primal-dual
solution.
We have also performed these reconstructions using a stabilized logarithmic barrier
algorithm based on the method presented in [50] and specialized to the present
reconstruction problem. Many of the computational features of our logarithmic barrier
implementation are identical to our primal-dual implementation, e.g., truncated
Newton, line search, computation of the gradient, Hessian-vector product, etc. For
a more detailed discussion, see [24]. The logarithmic barrier results are summarized
and compared against MAP-EM in Table 5.4. Termination of the logarithmic barrier
was dened by (4.29) and
These termination criteria for the logarithmic barrier correspond to roughly the same
accuracy as (4.28) and (4.29) do for the primal-dual method. Being a \long-step"
method, the logarithmic barrier gives the user less control over the exact stopping
point than does the \short-step" primal-dual. All of the logarithmic barrier reconstructions
in Table 5.4 used extrapolation. In all logarithmic barrier reconstructions,
was reduced by a factor of 10 between subproblems.
The eect of extrapolation is illustrated in Figures 5.1 and 5.2. In Figure 5.1,
the equivalent number of gradient evaluations (ngr) to reach termination is plot-
C.A. JOHNSON AND A. SOFER
100 150 200 250 300 350
PD
barrier
f-f*
ngr
Fig. 5.1. \Distance" from optimal solution at termination, as measured by dierence in objective
function f f  (where f  is here dened to be the lowest objective function obtained per study),
versus work required to reach termination, as measured by ngr, the equivalent number of gradient
evaluations. The studies included are those listed in Table 5.2. PD stands for non-extrapolated
primal-dual, PDX for extrapolated primal-dual.
ngr
PD
barrier
Fig. 5.2. Average value of  at subproblem termination versus average ngr (equivalent number
of gradient evaluations) for the seven studies listed in Table 5.2. PD stands for non-extrapolated
prmal-dual, PDX for extrapolated primal-dual.
ted against objective function \distance" f f  , the dierence between the function
value of the terminated solution and the lowest function value obtained for that recon-
struction. In all seven test cases (those listed in Tables 5.2{5.4), the unextrapolated
primal-dual method achieved the lowest objective function value. Thus, f f  is
zero for all unextrapolated primal-dual (PD) results but greater than zero for the
extrapolated primal-dual (PDX) and barrier results. The PDX results are clustered
in a region of lower ngr than the PD results. This indicates that extrapolation lowers
the computational expense to the solution at a slight deterioration in the nal
objective. Compared with the barrier method, either extrapolated or unextrapolated
primal-dual produces equivalent or better accuracy with less computation required.
In

Figure

5.2, the average number of equivalent gradient evaluations at subproblem
termination is plotted against the average value of  for each subproblem. Both
averages (ngr and ) were taken from the same seven test cases of Tables 5.2{5.4.
Compared with either unextrapolated primal-dual (PD) or extrapolated primal-dual
Study C
MAP-EM
LSEM
PD
|f-f*|
ngr
Fig. 5.3. Improvement in objective function as a function of gradient evaluations, Study F.
PDX denotes extrapolated primal-dual, PD denotes unextrapolated primal-dual, both using
Study F
MAP-EM
LSEM
PD
|f-f*|
ngr
Fig. 5.4. Improvement in objective function as a function of gradient evaluations, Study C.
PDX denotes extrapolated primal-dual, PD denotes unextrapolated primal-dual, both using
(PDX), the logarithmic barrier is clearly on a slower trajectory. The PD and PDX
trajectories are quite similar until approximately 0:01, at which point the PD
curve \swings out", while the PDX curve continues to descent log-linearly. This result
conrms that the prediction (extrapolation) step becomes more accurate near the so-
lution, resulting in more rapid convergence. However, a comparison of the objective
functions indicates that the value of PDX  is perhaps one step \ahead of itself,"
compared with the unextrapolated case.
The progress of the reconstruction on a study of a rat skull, Study C, is compared
for the various algorithms in Figure 5.3. The measure used is kf f  k (plotted on
a logarithmic scale). In the initial iterations DePierro Map EM and LSEM progress
rapidly and are ahead of the primal dual method. However the interior-point methods
rapidly reach the DePierro and LSEM objective values, and hence on, surpass them.
In the primal-dual methods depicted the value of  is 2. The methods achieve faster
initial progress using however the overall computational eort for full convergence
with this parameter setting is greater. The progress of the reconstruction in
another example, Study F, is compared in Figure 5.4.
We have also reconstructed a number of very large-scale \thin-slice" reconstructions
involving variables. Table 5.5 summarizes a number of properties of
22 C.A. JOHNSON AND A. SOFER

Table
Summary of thin-slice extrapolated results, including convergence measures and computational
costs to optimal solution.
study
f kr'k
F 1E-5 7,721,001 1.29E-9 2.87E-11 14 14 71 113
average ngr 115
these extrapolated primal-dual reconstructions at the converged solution. A smaller
group of datasets (the more visually \interesting" studies) were selected for the thin-
slice work, and certain reconstructions were repeated with dierent values of the
prior strength
. Thin-slice reconstructions seem to require a lower prior strength
than the corresponding thick-slice reconstructions. The most visually pleasing results
were from reconstructions using
which is 1/30 the prior strength that
was generally found to be most satisfactory in thick-slice reconstructions. The total
amount of work (as measured in ngr) required to reach termination in Table 5.5 is also
quite pleasing. The number of variables in a thin-slice reconstruction is approximately
3.7 times the number in thick-slice. The number of nonzero-valued measurements in
thin-slice mode is only marginally greater than in thick-slice mode, however, since
the number of counts is the same in both cases. These thin-slice reconstructions may
thus be better conditioned than their thick-slice counterparts.
In closing, we should comment that the tolerance we have used in our tests is
stricter than that usually necessary. Indeed, less accurate solutions may still give
acceptable images. When the EM method is applied to the (unregularized) ML ob-
jective, it is usually terminated after 50 or 100 iterations, and the images produced
are often good. Thus EM-ML remains a practical method that can sometimes reach
a solution of desirable image quality faster than an interior-point method. The difculty
with EM-ML is that its convergence is object-dependent [1]. Convergence in
areas of high activity amidst low activity or vice versa is notoriously slow, and a xed
termination rule based on (say) 50 or 100 iterations cannot guarantee acceptable image
quality. This has been observed in a number of reconstructions, including some
of high biomedical interest. In contrast to ML-EM, the primal-dual algorithm has
object-independent convergence characteristics. Furthermore it is
exible, and can be
adapted to solve a problem e-ciently both to the strict tolerance in the studies above
by setting a modest rate of decrease for the barrier parameter, say 2, and to a
looser tolerance by setting a more aggressive reduction rate such as
6. Conclusion. From the results of the previous Section, it is clear that the
primal-dual method can converge signicantly faster than the EM algorithm for regularized
ML reconstructions in emission tomography. The results also indicate that
the primal-dual method converges faster than the logarithmic barrier method. The
use of extrapolation in conjunction with the primal-dual method further reduces the
amount of computation required to achieve convergence.
Given that the negative regularized ML objective function that we minimize is
convex, approximately solving the reduced unsymmetric primal-dual Newton equa-
tions is appropriate. Symmetrizing the unsymmetric system, while potentially useful
for nonconvex problems, would in this case require solving for 2n variables without
avoiding the potential for ill-conditioning. Our stabilization technique avoids the
structural ill-conditioning of the condensed primal-dual matrix, and therefore solving
the reduced system poses no asymptotic di-culty as the barrier parameter approaches
zero. The computational e-ciency and relative simplicity of formation of the reduced
system of equations pose such a strong advantage that our choice of primal-dual
method almost seems obvious for this problem.
Since Newton's method converges quadratically near the solution, for a well-conditioned
system in the limit as  ! 0, one truncated-Newton step per subproblem
should yield an increasingly accurate and well scaled direction to the subproblem
solution for  k . As  is decreased, the subproblem solutions should become \close" to
each other for a convex problem [14]. Yet, the example in Table 4.1 illustrates that the
direction produced by the early-terminated CG can in fact become less accurate for
smaller  due to the structured ill-conditioning in M . In practice, we do not require
the accuracy of the test example in Table 4.1. Our termination conditions are dened
to be near the point on the trajectory where the stabilization approximation becomes
accurate enough to guarantee descent. These termination criteria are quite accurate
by the standards of the tomography community. Thus, although most reconstruction
problems are unlikely to be severely aected by ill-conditioning, the potential for slow
convergence near the solution due to ill-conditioning does exist. Our experience has
been that stabilization has been an eective safeguard against poor performance for
small values of the barrier parameter.
7.

Acknowledgments

. The study utilized the high-performance computational
capabilities of the IBM RISC/6000 SP system at the Division of Computer Research
and Technology, National Institutes of Health, Bethesda, MD. We are grateful to
Jurgen Seidel of the Department of Nuclear Medicine, National Institutes of Health,
for kindly providing us with the small animal data and Monte Carlo simulation data.
Our thanks go to two anonymous referees and the associate editor for their careful
reading and helpful comments.



--R

Noise properties of the EM algorithm: I.
Conjugate gradient methods for Toeplitz systems
A general class of preconditioners for statistical iterative reconstruction of emission computed tomography
Preconditioning methods for improved convergence rates in iterative reconstructions
A primal-dual algorithm for minimizing a nonconvex function subject to bound and linear equality constraints

Image reconstruction and restoration: overview of common estimation structures and problems
Maximum likelihood from incomplete data via the EM algorithm

On the convergence of an EM-type algorithm for penalized likelihood estimation in emission tomography
Numerical stability and e-ciency of penalty algorithms
Sequential Unconstrained Minimization Techniques

Stability of symmetric ill-conditioned systems arising in interior methods for constrained optimization
Matrix Computations
On use of the EM algorithm for penalized likelihood estimation


A generalized EM algorithm for 3-D Bayesian reconstruction from Poisson data using Gibbs priors
Image Reconstruction from Projections: the Fundamentals of Computerized Tomography
Methods of conjugate gradients for solving linear systems
Accelerated image reconstruction using ordered subsets of projection data
A practical interior-point method for convex programming
Nonlinear optimization for Volume
A system for the 3D reconstruction of retracted-septa PET data using the EM algorithm
Evaluation of 3D reconstruction algorithms for a small animal PET camera
Implementing and accelerating the EM algorithm for positron emission tomog- raphy


PET regularization for envelope guided conjugate gradients
Constrained reconstruction by the conjugate gradient method
A primal-dual interior point algorithm for linear programming
Probability measure estimation using
The importance of preconditioners in fast Poisson-based iterative reconstruction algorithms for SPECT
Practical tradeo
EM Reconstruction Algorithms for Emission and Transmission Tomography
A theoretical study of some maximum likelihood algorithms for emission and transmission tomography
Convergence of EM image reconstruction algorithms with Gibbs smoothing
A maximum a posteriori probability expectation maximization algorithm for image reconstruction in emission tomography

The superlinear convergence of a nonlinear primal-dual algorithm
Statistical Modeling and Fast Bayesian Reconstruction in Positron Tomog- raphy
Fast gradient-based methods for Bayesian reconstruction of transmission and emission PET images
Bayesian reconstruction of PET images: methodology and performance analysis
Analytic expressions for the eigenvalues and eigenvectors of the Hessian matrices of barrier and penalty functions
Preconditioning of truncated-Newton methods
Block truncated-Newton methods for parallel optimization




Barrier methods for large-scale quadratic programming
Maximum Likelihood Reconstruction for Emission Tomography
An infeasible interior-point method for linear complementarity problems

A statistical model for positron emission tomography
Noise properties of
Interior methods for constrained optimization


Stability of linear equation solvers in interior point methods.
--TR
applications of nonlinear programming;parallel computation;tomography;primal-dual methods;estimation;large-scale problems