We present a new algorithm for penalized-likelihood emission image reconstruction. The algorithm monotonically increases the objective function, converges globally to the unique maximizer, and easily accommodates the nonnegativity constraint and nonquadratic but convex penalty functions. The algorithm is based on finding paraboloidal surrogate functions for the log-likelihood at each iteration: quadratic functions that are tangent to the log-likelihood at the current image estimate, and lie below the log-likelihood over the entire nonnegative orthant. These conditions ensure monotonicity. The paraboloidal surrogates are maximized easily using existing algorithms such as coordinate ascent. Simulation results show that the proposed algorithm converges faster than the SAGE algorithm, yet the new algorithm is somewhat easier to implement.Peer Reviewedhttp://deepblue.lib.umich.edu/bitstream/2027.42/85943/1/Fessler152.pd

Erdogan, Hakan

Fessler, Jeffrey A.

English

Deep Blue at the University of Michigan

A Paraboloidal Surrogates Algorithm for Convergent
Penalized-Likelihood Emission Image Reconstruction

We present a new algorithm for penalized-likelihood emission image reconstruction. The algorithm monotonically increases the objective function, converges globally to the unique maximizer, and easily accommodates the nonnegativity constraint and nonquadratic but convex penalty functions. The algorithm is based on finding paraboloidal surrogate functions for the log-likelihood at each iteration: quadratic functions that are tangent to the log-likelihood at the current image estimate, and lie below the log-likelihood over the entire nonnegative orthant. These conditions ensure monotonicity. The paraboloidal surrogates are maximized easily using existing algorithms such as coordinate ascent. Simulation results show that the proposed algorithm converges faster than the SAGE algorithm, yet the new algorithm is somewhat easier to implement. I. INTRODUCTION  Statistical methods for tomographic image reconstruction are increasingly being used clinically. For widespread acceptance of statistica..

Jeffrey A. Fessler

Hakan Erdogan

CiteSeerX

In Conference Record of 1998 IEEE Nuc. Sci. Symp. and Med. Im. Conf., Vol. 2, pp. 1132-5
A Paraboloidal Surrogates Algorithm for Convergent
Penalized-Likelihood Emission Image Reconstruction
Jeffrey A. Fessler and Hakan Erdo˘ gan
4240 EECS Bldg., University of Michigan, Ann Arbor, MI 48109-2122
Abstract
We present a new algorithm for penalized-likelihood
emission image reconstruction. The algorithm monotonically
increases the objective function, converges globally to the
unique maximizer, and easily accommodates the nonnegativity
constraint and nonquadratic but convex penalty functions.
The algorithm is based on ﬁnding paraboloidal surrogate
functions for the log-likelihood at each iteration: quadratic
functions that are tangent to the log-likelihood at the current
image estimate, and lie below the log-likelihoodover the entire
nonnegative orthant. These conditions ensure monotonicity.
The paraboloidal surrogates are maximized easily using
existing algorithms such as coordinate ascent. Simulation
results show that the proposed algorithm converges faster than
the SAGE algorithm, yet the new algorithm is somewhat easier
to implement.
I. INTRODUCTION
Statistical methods for tomographic image reconstruction
are increasingly being used clinically. For widespread
acceptance of statistical methods, fast converging algorithms
are essential. This paper summarizes a new algorithm for
reconstructing emission images from Poisson sinogram
measurements. The algorithm
 monotonicallyincreases the objective function,
 is globallyconvergent,
 naturally accommodates the nonnegativityconstraint,
 requires less CPU time per iteration than the space-
alternating generalize EM (SAGE) algorithm [1].
For reconstruction problems (such as 2D PET and SPECT)
where the system matrix G can be precomputed and stored,
we recommend this new algorithm over our previously
published algorithms for penalized-likelihood reconstruction.
For reconstruction problems where the system matrix is
represented in factored form [2], methods that can better
exploit this representation, such as the conjugate-gradient
algorithm, appear to remain preferable.
II. THE PROBLEM
We assume the usual Poisson statistical model for the
emission measurements Yi, i =1 ;:::;n d:
Yi  Poisson

 yi(
true)
	
(1)
 yi()=
np X
j=1
aijj + ri; (2)
where j is the unknown activity in the jth pixel,
j =1 ;:::;n p,t h eri’s denote background events such as
random coincidences and scatter, and A = faijg represents
the system model, including ray-dependent factors such as
attenuation and detector efﬁciency. We assume the ri’s and
aij’s are known nonnegative constants. We assume ri > 0;
which is reasonable and realistic since any real PET scan
will have nonzero randoms. (It also turns out to simplify our
algorithm.) We wish to estimate  =[ 1;:::; np] from a
measurement realization fYi = yig
nd
i=1. The log-likelihoodis:
L() 
nd X
i=1
yi log  yi() −  yi(); (3)
neglectingconstantsindependentof. Forpenalized-likelihood
image reconstruction, one seeks the image that maximizes an
objective function as follows:
^  =a r gm a x
 0
(); ()=L() − R(); (4)
where R() is a roughness penalty included for regularization.
The resolution/noisetradeoff is controlledby the regularization
parameter  [3]. The problem is to ﬁnd the maximizer ^ 
using an algorithm that converges rapidly and that requires
as little CPU time per iteration as possible. Usually these
are conﬂicting requirements, since often one can modify an
algorithm to converge in fewer iterations, but at the expense of
more work per iteration.In Conference Record of 1998 IEEE Nuc. Sci. Symp. and Med. Im. Conf., Vol. 2, pp. 1132-5
0 2 4 6 8 10
−20
−15
−10
−5
0
5
Parabolic Surrogate Function
h
i
(
l
)
a
n
d
q
i
(
l
;
l
n
i
)
Log-likelihoodhi(l)
Surrogateqi(l;ln
i )
ln
i l !
Figure 1: Illustration of 1D parabolic surrogate function. Note that qi(l;l
n
i )  hi(l) for l  0.
III. THE NEW ALGORITHM
We had previously recommended the SAGE algorithm [1]
as a fast globally-convergent algorithm for this problem.
The construction of the SAGE algorithm requires certain
minimizations that are somewhat unusual in the tomographic
literature, and must be implemented carefully to achieve
reasonable CPU time per iteration.
The new algorithm we propose is based on the following
simple idea, illustrated in one dimension in Fig. 1. Since
the log-likelihood L() is difﬁcult to maximize directly, we
endeavor to ﬁnd a surrogate function Q(;
n) that is easier to
maximize, and maximize that functionat the nth iteration,i.e.:

n+1 =a r gm a x
 0
Q(;
n) − R() (5)
If we choose the sequence of surrogate functions Q(;
n)
properly, then the sequence of iterates f
ng will converge to
the maximizer ^ . The SAGE algorithm is indirectly based
on this idea; the expected conditional log-likelihood of the
“complete” data space given the observed data is a type of
surrogate function that indeed satisﬁes the conditions sufﬁcient
to ensure convergence [1]. However, the statistical construction
of the surrogate functions for the SAGE and other EM
algorithms can seem somewhat mysterious.
In this paper, we propose a new approach to constructing
surrogatefunctions that uses only basic calculus principles. We
can rewrite the log-likelihoodL() in (3) as follows:
L()=
nd X
i=1
hi([A]i)
where
[A]i =
np X
j=1
aijj
is the projection of the object along the ith ray, and
hi(l)=yi log(l + ri) − (l + ri)
is the marginal log-likelihood of the ith measurement. The hi
functions are concave, and strictlyconcave if yi > 0.
Our new proposed strategy for choosing the surrogate
function Q is to ﬁrst ﬁnd one-dimensional parabolic surrogate
functions qi(l;ln
i ), as illustratedin Fig. 1, where
ln
i
, [A
n]i:
We then combine these 1D surrogate functions to form an
overall surrogate function as follows:
Q(;
n)
,
nd X
i=1
qi([A]i;ln
i ): (6)
Provided that the 1D surrogate functions satisfy the following
three properties:
 qi(ln
i ;ln
i )=hi(ln
i );In Conference Record of 1998 IEEE Nuc. Sci. Symp. and Med. Im. Conf., Vol. 2, pp. 1132-5
 _ qi(ln
i ;ln
i )=_ hi(ln
i )
 qi(l;ln
i )  hi(l); 8l  0;
then it can be shown that the recursive algorithm given by (5)
will monotonically increase L(), and in fact can be shown to
converge globally by a proof similar to that in [1].
Using the fact that the hi’s are concave and that the ﬁrst
derivatives of the hi’s are convex, one can show [4] that the
above three conditions will be satisﬁed if we choose parabolic
surrogate functions as follows:
qi(l;ln
i )=hi(ln
i )+_ hi(ln
i )(l − ln
i ) −
1
2
ni(ln
i )(l − ln
i )2; (7)
where
ni(l)
,
( 2
l2
h
hi(l) − h(0) − l_ hi(l)
i
;l > 0
−¨ hi(l);l =0 :
(8)
The basic idea is illustrated in Fig. 1. The parabolic
surrogate function qi(l;ln
i ) has the same value as the marginal
log-likelihood hi(l) at the current projection value l = ln
i ,
and has the same slope at that point. This is evident from
(7). In addition, the parabolic function lies below hi(l) for
all nonnegative l. This is the key to having a monotonic
algorithm1. The proof that this qi choice satisﬁes the third
of the above three conditions is somewhat detailed, and
is described in [4]. When the parabolic qi functions are
“assembled” as in (6) to form an overall surrogate function,
the ﬁnal form of Q is a paraboloid, i.e. a quadratic form. One
can maximize this quadratic form by any number of different
methods. (Exact maximization is not necessary, due to the
iteration(5).) We have chosen to use successive over-relaxation
or coordinate ascent [5] for this maximization, since it easily
accommodates both the nonnegativity constraint as well as
nonquadraticpenalty functions.
In words, the algorithm alternates between two steps. The
ﬁrst step is to ﬁnd the coefﬁcients of the surrogate parabolas in
(7) using (8). This involves a single-pass over the sinogram
with trivial computation. Then the paraboloid (6) is partially
maximized by a single cycle of coordinate ascent. This
step requires roughly the equivalent of one forward and one
backprojection, similar to most iterative algorithms.
IV. RESULTS
We used the same brain emission simulation reported in [1]
to evaluate the proposed algorithm. In [1], the SAGE algorithm
was compared to many alternatives, includingGreen’s one-step
late (OSL) method [6], Kaufman’s boundedline search [7], and
the generalized EM (GEM) algorithm of Hebert and Leahy [8].
We found that the SAGE algorithm converged faster than all of
theabovemethods,so herewefocusoncomparingtheproposed
algorithm to just the SAGE algorithm.
1One could also ﬁnd approximating parabolas using Newton’s
method, but these parabolas can cross hi(l), and the result is a
nonmonotonic algorithm that is not guaranteed to converge. Our
construction using (8) avoids this problem.
Figure 2 plots the increase in the log-likelihood
(
n) − (
0) versus CPU time on a DEC AlphaStation
600 5/333 workstation. The proposed algorithm (paraboloidal
surrogates coordinate ascent (PSCA) algorithm) converges a
little bit faster than SAGE, in part because it requires about
10% less CPU time per iteration, and in part because it
increases  more each iteration. The differences are modest
since both algorithms converge quite quickly, so there is
limited room remaining for improvement. The new algorithm
is simpler to implement than SAGE as well.
Theresultingimages areindistinguishablefromthoseshown
in[1]sincewe used thesame objectivefunctionand bothPSCA
and SAGE algorithms are globally convergent.
V. DISCUSSION
The proposed algorithm is based on the “optimization
transfer” principle. Since the original objective function  is
cannot be maximized directly, we instead maximize a sequence
of surrogate functions (;
n). The key is to choose surrogate
functions that are easier to maximize than , but have low
curvature (high curvature surrogate functions lead to slow
convergence rate [4, 9]). Many algorithms in the literature
are based (implicitly or explicitly) on optimization transfer
ideas, include EM, SAGE, grouped coordinate ascent [10],
the convex algorithm [11], and ISRA [12–14]. In almost all
cases, thesurrogatefunctions are separable, which makes them
trivial to maximize, but also means that they have very high
curvature and hence poor convergence rate. The paraboloidal
surrogate functions that we have proposed here and in [4] are
the ﬁrst that we are aware of that are nonseparable. Generally
nonseparable functions are harder to maximize. Fortunately, a
notable exception is quadratic surrogate functions, which is the
choice we have made. This choice may be suboptimal; it is
possible that there exist other nonseparable surrogate functions
that are easily maximized but have even lower curvature and
hence yield faster convergence rates.In Conference Record of 1998 IEEE Nuc. Sci. Symp. and Med. Im. Conf., Vol. 2, pp. 1132-5
0 1 2 3 4 5
200
300
400
500
600
700
800
900
CPU Seconds, DEC AlphaStation 600 5−333
F
(
l
n
)
−
F
(
l
0
)
Penalized−Likelihood
FBP Initialization
Quadratic Penalty
35% Randoms
PSCA
SAGE
Figure 2: Comparison of convergencerate of proposed paraboloidalsurrogate coordinate ascent(PSCA) algorithm versus SAGE algorithm.
VI. REFERENCES
[1] J. A. Fessler and A. O. Hero, “Penalized maximum-
likelihood image reconstruction using space-alternating
generalized EM algorithms,” IEEE Tr. Im. Proc.,v o l .4 ,
no. 10, pp. 1417–29, October 1995.
[2] J. Qi, R. M. Leahy, S. R. Cherry, A. Chatziioannou,
and T. H. Farquhar, “High resolution 3D Bayesian image
reconstructionusingthemicroPET small-animalscanner,”
Phys. Med. Biol., vol. 43, no. 4, pp. 1001–14, April 1998.
[3] J. A. Fessler and W. L. Rogers, “Spatial resolution
properties of penalized-likelihood image reconstruction
methods: Space-invariant tomographs,” IEEE Tr. Im.
Proc., vol. 5, no. 9, pp. 1346–58, September 1996.
[4] H. Erdo˘ gan and J. A. Fessler, “Monotonic algorithms for
transmission tomography,”I E E ET r .M e d .I m . , vol. 18, no.
9, pp. 801–14, September 1999.
[5] J. A. Fessler, “Penalized weighted least-squares image
reconstruction for positron emission tomography,” IEEE
T r .M e d .I m . , vol. 13, no. 2, pp. 290–300, June 1994.
[6] P. J. Green, “Bayesian reconstructions from emission
tomography data using a modiﬁed EM algorithm,” IEEE
T r .M e d .I m . , vol. 9, no. 1, pp. 84–93, March 1990.
[7] L. Kaufman, “Implementing and accelerating the EM
algorithm for positron emission tomography,” IEEE Tr.
Med. Im., vol. 6, no. 1, pp. 37–51, March 1987.
[8] T. Hebert and R. Leahy, “A generalized EM algorithm
for 3-D Bayesian reconstruction from Poisson data using
Gibbs priors,” IEEE Tr. Med. Im., vol. 8, no. 2, pp. 194–
202, June 1989.
[9] H. Erdo˘ gan and J. A. Fessler, “Accelerated monotonic
algorithms for transmission tomography,” in Proc. IEEE
Intl. Conf. on Image Processing, volume 2, pp. 680–4,
1998.
[10] J. A. Fessler, E. P. Ficaro, N. H. Clinthorne,and K. Lange,
“Grouped-coordinate ascent algorithms for penalized-
likelihood transmission image reconstruction,” IEEE Tr.
Med. Im., vol. 16, no. 2, pp. 166–75, April 1997.
[11] K. Lange and J. A. Fessler, “Globally convergent
algorithms for maximum a posteriori transmission
tomography,” IEEE Tr. Im. Proc., vol. 4, no. 10, pp. 1430–
8, October 1995.
[12] M. E. Daube-Witherspoon and G. Muehllehner, “An
iterativeimagespace reconstructionalgorithmsuitablefor
volume ECT,” IEEE Tr. Med. Im., vol. 5, no. 2, pp. 61–66,
June 1986.
[13] J. M. Ollinger, “Iterative reconstruction-reprojection and
the expectation-maximization algorithm,” IEEE Tr. Med.
Im., vol. 9, no. 1, pp. 94–98, March 1990.
[14] A. R. De Pierro, “On the relation between the ISRA
and the EM algorithmfor positronemission tomography,”
I E E ET r .M e d .I m . , vol. 12, no. 2, pp. 328–333,June 1993.
See http://www.eecs.umich.edu/fessler.
This work was supported in part by NIH grants CA-60711
and CA-54362, and by the Whitaker Foundation.

A Paraboloidal Surrogates Algorithm for Convergent Penalized-Likelihood Emission Image Reconstruction

We present a new algorithm for penalized-likelihood emission image reconstruction. The algorithm monotoni-cally increases the objective function, converges globally to the unique maximizer, and easily accommodates the nonnegativity constraint and nonquadratic but convex penalty functions. The algorithm is based on finding paraboloidal surrogate functions for the log-likelihood at each iteration: quadratic functions that are tangent to the log-likelihood at the current image estimate, and lie below the log-likelihood over the entire nonnegative orthant. These conditions ensure monotonicity. The paraboloidal surrogates are maximized easily using existing algorithms such as coordinate ascent. Simulation results show that the proposed algorithm converges faster than the SAGE algorithm, yet the new algorithm is somewhat easier to implement. I

A paraboloidal surrogates algorithm for convergent penalized-likelihood emission image reconstruction

http://people.sabanciuniv.edu/~haerdogan/pubs/fessler98aps.pdf

A Paraboloidal Surrogates Algorithm for Convergent Penalized-Likelihood Emission Image Reconstruction

Abstract

Similar works

Full text

Available Versions

Deep Blue at the University of Michigan

CiteSeerX

CiteSeerX