Superiorization of Preconditioned Conjugate Gradient Algorithms for Tomographic Image Reconstruction
Elias S. Helou Gabor T. Herman Chuan Lin Marcelo V. W. Zibetti


Properties of Superiorized Preconditioned Conjugate Gradient (SupPCG)
algorithms in image reconstruction from projections are examined.
Least squares (LS) is usually chosen for measuring data-inconsistency
in these inverse problems. Preconditioned Conjugate Gradient algorithms
are fast methods for finding an LS solution. However, for ill-posed
problems, such as image reconstruction, an LS solution may not provide
good image quality. This can be taken care of by superiorization.
A superiorized algorithm leads to images with the value of a secondary
criterion (a merit function such as the total variation) improved
as compared to images with similar data-inconsistency obtained by
the algorithm without superiorization. Numerical experimentation shows
that SupPCG can lead to high-quality reconstructions within a remarkably
short time. A theoretical analysis is also provided.


A Common Source for Scalars:

Axiflavon-Higgs Unification
    Florian Goertz
    Received 22 May 2018 / Accepted 17 July 2018
=========================================================





§ INTRODUCTION


Superiorization <cit.> is an algorithmic
framework that consists in perturbing iterative methods. It can be
applied to feasibility problems (problems of finding a point belonging
to a given set) or to optimization problems (problems of finding the
best point among a set of candidates, according to some criterion).
Its aim is to improve the iterates, according to a secondary criterion
(which is some sort of merit function), produced by an iterative method
for solving the original problem without the secondary criterion.

Many problems in science or engineering can be mathematically modeled
as feasibility <cit.> or optimization<cit.>
problems. In this paper we use tomographic imaging as our example
application of the superiorized algorithms that we introduce, but
the basic idea is useful in numerous other contexts.[The reader can find a continuously updated bibliography about superiorization
at http://math.haifa.ac.il/YAIR/bib-superiorization-censor.html.]

We investigate a superiorized Preconditioned Conjugate Gradient (PCG)
method, applying it to tomographic imaging and comparing it to other
approaches to this application. In Section <ref> we
connect tomographic image reconstruction to optimization, as in <cit.>.
Section <ref> reviews superiorization, the Conjugate Gradient (CG)
and the PCG algorithms. Superiorization of these algorithms results
in the algorithms Superiorized Conjugate Gradient (SupCG) and Superiorized
Preconditioned Conjugate Gradient (SupPCG), and superiorization of
one of the Algebraic Reconstruction Techniques (ART) <cit.>
leads to SupART. Section <ref> reports on numerical
experiments comparing the performance in image reconstruction of these
algorithms, as well as a version of ART that uses blob basis functions <cit.>.
Section <ref> presents a theoretical analysis
of the SupPCG algorithm. Finally, Section <ref>
provides the concluding remarks.




§ OPTIMIZATION IN TOMOGRAPHIC RECONSTRUCTION




Assume that the image x∈ℝ^n generates the data b∈ℝ^m
through

    Rx=b_exact+e,

where R∈ℝ^m× n is the projection (also called
Radon) matrix, b=b_exact+e, and e∈ℝ^m
is the unknown error vector introduced into the measurements. Therefore,
we need to solve Rx≈ b, where the meaning of the approximation
must be well defined mathematically. There are several ways of doing
this, in the present paper we follow the least squares approach:

    x∈_y∈ℝ^nRy-b^2.

Because the least squares approach is well accepted and extremely
general, it is not surprising that several techniques have been developed
for its solution. The minimizers in (<ref>) are
given by the solutions of the normal equations

    R^TRx=R^Tb,

which is a linear system of equations. Due to the huge size and sparsity
of the projection matrix (as can be the case in image reconstruction),
iterative methods that successively approximate a solution may be
the most practical.

Among the best-known iterative techniques for the general problem
of least squares, one can point at CG and at PCG, to be applied to
the system (<ref>); see <cit.>.
For tomographic image reconstruction, algebraic reconstruction techniques (ART) <cit.>
are widely applied.

To put this into a general context, let f be a function defined
over images such that f(x) is some measure of the inconsistency
of x with the given data. The primary aim of a proposed reconstruction
method should be to produce an x for which the value of f is
relatively small, therefore f is referred to as the primary criterion.
Here we concentrate on the squared error, f(x)=Rx-b^2, as
the primary criterion. We also discuss a primary criterion based on
a Bayesian approach to image reconstruction <cit.>, namely

    r^2‖ Rx-b‖ ^2+‖ x-μ_X‖ ^2,

where the number r is the so-called signal-to-noise ratio and μ_X
is a uniformly gray image with the gray value being the average value
of the image as can be accurately estimated based on all the measured
data <cit.>. Both these need to be specified
for a particular reconstruction task.




§ SUPERIORIZATION OF PCG AND ART







 §.§ Superiorization


As stated in <cit.>:

The superiorization methodology is used for improving the efficacy
of iterative algorithms whose convergence is resilient to certain
kinds of perturbations. Such perturbations are designed to “force”
the perturbed algorithm to produce more useful results for the intended
application than the ones that are produced by the original iterative
algorithm. The perturbed algorithm is called the “superiorized version”
of the original unperturbed algorithm. If the original algorithm is
computationally efficient and useful in terms of the application at
hand and if the perturbations are simple and not expensive to calculate,
then the advantage of this method is that, for essentially the computational
cost of the original algorithm, we are able to get something more
desirable by steering its iterates according to the designed perturbations.

Several iterative algorithms use an updating of the current approximation
to the problem at hand that is of the form

    (x_k+1,u_k+1^1,…,u_k+1^ν)=𝑈(x_k,u_k^1,…,u_k^ν),

where the u_k^i for 1≤ i≤ν are some auxiliary vectors
and x_k is the image at iteration k. A superiorized version
of such a method can have the form

    (x_k+1,u_k+1^1,…,u_k+1^ν)=𝑈(x_k+1/2,u_k^1,…,u_k^ν),

where x_k+1/2=x_k+s_k, with {s_k} referred
to as a superiorization sequence. For our experiments in Section <ref>
we used the method from <cit.>, as specified in Algorithm <ref>,
for generating the superiorization sequence. The actual values produced
by that specification depend on the chosen secondary criterion and
the procedure to compute “nonascending vectors” for this criterion;
here we report on results for total variation (TV) as the secondary
criterion (see <cit.>). 

Total variation is defined for a vector y∈ℝ^n
that represents a discrete image [y] with l rows and c columns
of pixels with lc=n. The value [y]_i,j of the pixel in the
ith row and jth column of this image is defined to be the (c(i-1)+j)th
component of vector y. The definition of total variation is

    TV(y):=∑_i=1^l-1∑_j=1^c-1√(([y]_(i,j)-[y]_(i+1,j))^2+([y]_(i,j)-[y]_(i,j+1))^2).



Algorithm <ref> uses nonascending
vectors. Following <cit.>, given a point y∈ℝ^n,
we say that t∈ℝ^n is a nonascending vector for TV
at y if t≤1 and there is Δ>0 such that

    δ∈[0,Δ]⇒ TV(y+δ t)≤ TV(y).

We denote the set of nonascending vectors for TV at y as ∂̃TV(y).
Note that ∂̃TV(y) is never empty because the zero
vector is in ∂̃TV(y). In order to get the t∈∂̃TV(y)
that we actually use, first define t∈ℝ^n
componentwise. For any component of t, its value is
the negative of the partial derivative of TV at y with respect
to that component, if this partial derivative is well defined, and
is 0 otherwise. We finally define the nonascending vector t∈∂̃TV(y)
that we use as

    t=t/t,    if ‖t‖≠0,
    the zero vector,    otherwise.

That the t specified this way is indeed a nonascending vector for
TV at y is an immediate consequence of <cit.>.

The specification of Algorithm <ref>
involves further parameters: K∈ℕ (the number of nonascending
steps), a∈(0,1) (the step-length diminishing factor for each
nonascending trial) and γ∈(0,∞) (the starting step-length);
for complete specification, these parameters need to be selected by
the user. Next we provide details of the purpose and operation of
Algorithm 1.




Algorithm <ref> specifies the s_k to be
used to define x_k+1/2=x_k+s_k, which is then
further used to provide us with the next image x_k+1 of the superiorized
iterative algorithm, see (<ref>). As the iterations
proceed, changes are made not only to the image vector but also to
an integer variable ℓ (see Step 7 of Algorithm 1), these changes
are needed to make the superiorized algorithm behave as desired (see
<cit.> and Section <ref>
below). In the superiorized algorithm, ℓ is initialized to be
0 (that is why γ is referred to as the starting step-length).
In each step (<ref>) of the superiorized iterative
algorithm, Algorithm 1 is called with the current values of the image
vector x_k, the integer ℓ and the user-specified parameters
a, γ and K (these do not change during an execution
of the superiorized algorithm). Algorithm <ref>
returns s_k(=y-x_k) and the new value of ℓ,
to be used in the next iterative step of the algorithm. During the
execution of the superiorized algorithm, the step lengths γ̃
form a summable sequence as needed in the theoretical results on the
behavior of superiorized algorithms (see, for example, <cit.>).

Superiorization has, by design, an influence over algorithmic convergence,
but ideally it maintains the required properties of the limit point.
The influence of the superiorization should be that the sequence of
iterates {x_k} (or {x_k+1/2}) for the superiorized
method has reduced secondary criterion (in our case, TV) when compared
to the equivalent sequence produced by the unsuperiorized method.
For our method this is true because the superiorization step, although
only guaranteed not to increase the value of the undifferentiable
TV criterion, is frequently capable of actually reducing the value
of the TV. That is, one usually finds that TV(x_k+1/2)<TV(x_k).

We provide a theoretical discussion of the convergence of the SupPCG
algorithm in Section <ref>. A convergence
proof of the SupCG algorithm appears in <cit.> and that result
is extended naturally to the superiorized version of the PCG algorithm,
which is a method that we now proceed to discuss.




 §.§ Preconditioned Conjugate Gradient


The CG algorithm is designed to solve a linear system of equations
Ax=y, where A is symmetric positive semi-definite. If A=R^TR
and y=R^Tb, then CG can solve (<ref>) via (<ref>).
The CG method may converge slowly depending on the distribution of
the eigenvalues of the system matrix A. The idea of PCG is to replace
the system Ax=y by a preconditioned system MAx=My, where M
is a symmetric positive-definite matrix, so that the two systems have
the same solutions but MA has better spectral properties for the
application of the CG algorithm. We return to this topic later in
Section <ref>.




Given U_PCG defined by Algorithm <ref>,
the update for PCG algorithm is

    (x_k+1,p_k,h_k)=U_PCG(x_k,p_k-1,h_k-1),

where x_0∈ℝ^n is an arbitrary vector, but x_1,
p_0 and h_0 must follow precise rules, as described for
SupPCG in the next subsection. The preconditioning matrix M that
we use here is of the form F^-1DF, where F is the Discrete
Fourier Transform (DFT) and D is a diagonal matrix that represents
a generalized Hamming window <cit.> (using parameters
μ and ρ) combined with a ramp filter for frequency ω
with -π≤ω≤π:

    h(ω)=(|ω|+μ)·(ρ+(1-ρ)cosω).






 §.§ Superiorized Preconditioned Conjugate Gradient


In this subsection we describe, as Algorithm <ref>, the
SupPCG method, which is the main subject of the present paper. In
its specification, we use the notations f(x)=Rx-b^2, A=R^TR
and y=R^Tb; see (<ref>) and (<ref>).
The inputs of Algorithm <ref> are the initial image x_0,
the a, γ and K that have the same roles as in Algorithm <ref>
and a user-specified ε>0 that determines the termination
of Algorithm <ref> (see Step 10). The algorithm returns
the number of iterations k at the time of its termination and an
output image x_k-1/2.




For the termination of Algorithm <ref> to be guaranteed
by the results we present in Section <ref>,
we will need appropriate bounds on the elements of the superiorization
sequence { s_k}; we now provide such bounds.

Upon the execution of Step 11 in Algorithm <ref>,
we have that s_k≤ Kγ a^(k-1)K. 

To verify this assertion, we first note that the
vectors t∈∂̃TV(y) in Algorithm 1 satisfy t≤1,
by the definition of nonascending vectors as given above. Next we
prove inductively that ℓ_k≥(k-1)K. This inequality holds
for k=1 as an equality. Now assume, for the induction, that it
holds for iteration k. The value of ℓ_k+1 is determined
by Step 11 of Algorithm <ref>, which consists in a call
S_TV(x_k,ℓ_k,a,γ,K) to Algorithm <ref>.
In each of the K loops of Algorithm <ref>
(see its Step 2), ℓ is increased at least once (in Step 7).
This and the induction hypothesis implies that for the ℓ_k+1
returned by S_TV(x_k,ℓ_k,a,γ,K),
we have ℓ_k+1≥ℓ_k+K≥(k-1)K+K=kK=((k+1)-1)K;
this completes the inductive proof.

We introduce some extra notations. Consider the
call S_TV(x_k,ℓ_k,a,γ,K) to Algorithm <ref>
in Algorithm <ref>. In the execution of that call, there
is an inner loop (indexed by i) consisting of Steps 3-9 of Algorithm <ref>.
We now use t_k,i to denote the value of t obtained by Step
3 in that inner loop and ℓ_k,i as ℓ-1 at the time when
the condition TV(ỹ)≤ TV(y) in Step 8 of the inner loop
is satisfied. Looking at S_TV(x_k,ℓ_k,a,γ,K),
we see that for the s_k returned by Step 11 of Algorithm <ref>
we have that s_k=γ∑_i=1^Ka^ℓ_k,it_k,i. Since
‖ t_k,i‖≤1, ℓ_k,i≥ℓ_k≥(k-1)K
and a∈(0,1), we conclude that ‖ s_k‖≤γ∑_i=1^K‖ a^ℓ_k,it_k,i‖≤γ∑_i=1^Ka^ℓ_k,i‖ t_k,i‖≤γ∑_i=1^Ka^ℓ_k≤ Kγ a^(k-1)K.




 §.§ Superiorized Algebraic Reconstruction Technique


The Algebraic Reconstruction Techniques consist of sequentially projecting
the current image toward the hyperplanes defined by the equations
of a linear system and we consider a full cycle through the equations
to be a single iteration. Several variants are possible, we use the
version originally proposed in <cit.> with pixel basis functions
(we refer to this as ART) and also with blob basis functions <cit.>,
see alternatively <cit.> (we refer to this
as BlobART). We also report on superiorization using the
TV secondary criterion for the pixel basis, we refer to the so
superiorized ART as SupART.

To be more specific, the ART method proposed in <cit.> (see
alternatively <cit.>) aims at finding
the minimizer of the primary criterion given in (<ref>).
Therefore the signal-to-noise ratio r is a parameter that needs
to be specified for both ART and for BlobART. There is an additional
parameter λ for these algorithms, it is the so-called relaxation
parameter that determines the size of the step in projecting the current
image toward the next hyperplane. Mathematical convergence theory
allows the relaxations to vary with iterations <cit.>,
and this may be useful in practice, but here we choose a single λ
for a reconstruction run.




§ NUMERICAL EXPERIMENTATION







 §.§ Tomographic data acquisition simulation


We used data obtained from simulation of a phantom that includes features
such as seen in medical images <cit.>.
The phantom has low-contrast features of interest inside a high attenuation
skull-like structure. If shown in its full grayscale range, the inner
features are nearly invisible, therefore the phantom and reconstructions
are displayed using a gray scale in which linear attenuation values
smaller than 0.204 cm^-1 are displayed as black and values
higher than 0.21675 cm^-1 appear as white. Phantom and
tomographic data were obtained using SNARK14.[SNARK14 may be downloaded free of charge from http://turing.iimas.unam.mx/SNARK14M/.]
(SNARK14 is the most recent version of a series of software packages
for the reconstruction of 2D images from 1D projections. From the
point of view of the current paper, its main advantage over the previous
version SNARK09 <cit.> is the new ability to superiorize automatically
any iterative algorithm.) The phantom can be seen at the bottom-right
of Fig. <ref>. It consists of an array
of 243×243 pixels, each of which corresponds to an area of
0.0752×0.0752 cm^2. Data were collected for 360
equally spaced angles in [0,π) and each of these angles was sampled
for 345 parallel rays spaced at 0.0752 cm between them.
Random photon emission at the x-ray source and scattering of the x-ray
photons was simulated.







 §.§ Algorithmic parameter selection


The compared algorithms, with the exception of CG, have parameters
to be tuned. To select these parameters, we evaluated the quality
of the reconstructed image x using the figure of merit called selective
error (SE) and then chose, for each algorithm separately, the parameters
that provided the best reconstructed image within the first 15
iterations. SE is defined by

    SE(x)=C√(∑_(i,j)∈ S([x]_i,j-[x^†]_i,j)^2),

where C is a constant dependent on the phantom x^† and
S is the set of pixels contained in a centered ellipse with a 5 cm
horizontal and a 7 cm vertical semi-axis. This ellipse is just
inside the “skull” in the phantom.

All the superiorized methods share three parameters (K, a and
γ, see Algorithm 1 and 3 ) that were drawn from the following
sets: K∈{10,20,40}, a∈{1-10^-5,1-10^-4,…,1-10^-1},
γ∈{10^-2,5·10^-2}. PCG and SupPCG have two preconditioning
parameters: μ∈{10^-5,10^-4,10^-3} and ρ∈{0.4,0.6,0.8};
see (<ref>). ART, BlobART and SupART use the
signal-to-noise ratio r∈{5,10}, which corresponds to the parameter
with the same symbol in <cit.> (where the version
of ART we use here is described) and a relaxation parameter λ∈{10^-2,5·10^-2,10^-1,5·10^-1,1};
see also <cit.>. The combination of parameters,
within the sets specified above, that provided the best results are
listed in Table <ref>.
















 §.§ Numerical results


In Fig. <ref> we present the evolution of the
SE of (<ref>) as iterations proceed, while Figure <ref>
exhibits a time-wise comparison among the methods. The computer used
to run the methods had an Intel i7 7700HQ processor running at up
to 3.4GHZ with 32GB of RAM available for the computations. We notice
that a superiorized method obtains better values for the SE figure
of merit than the corresponding unsuperiorized version for almost
every iteration. This does not come as a surprise, for two reasons.
First, the superiorized version is supposed to incorporate better
prior knowledge about the sought-after image to be reconstructed through
the secondary criterion. Second, in the experiments we present, the
algorithmic parameter settings were tailored to obtain an improved
value of the figure of merit SE.




§ THEORETICAL ANALYSIS


In the present section we prove that the SupPCG method (Algorithm <ref>)
terminates under a suitable hypothesis on the superiorization sequence
{ s_k} that is defined by calling S_TV(x_k,a,γ,K).
Our proof is based on facts established in <cit.>, which contains
a proof that SupCG terminates. The needed relation between the SupPCG
and the SupCG is established via the so-called Transformed PCG algorithm
(TPCG) and its superiorized version (SupTPCG). Next we discuss these
algorithms.

Since the preconditioning matrix M is symmetric positive-definite,
there exists an invertible matrix N such that M=N^TN and,
for such a matrix, the eigenvalues of MA are the same as the eigenvalues
of NAN^T (if v is an eigenvector of MA, then N^-Tv
is an eigenvector of NAN^T with the same associated eigenvalue).
Applying one step of the CG algorithm (which is exactly Algorithm <ref>
with M the identity matrix) to the linear system NAN^Tx̂=Ny
provides one step of TPCG, as given in Algorithm <ref>.
This step is also called from SupTPCG, see Step 13 of Algorithm <ref>.

There is a strict equivalence between Algorithm <ref>
and Algorithm <ref>. To see this consider Table <ref>.
We claim, and this can be verified by following step-by-step SupPCG(x_0,a,γ,K,ε)
and SupTPCG(N^-Tx_0,a,γ,K,ε), that
the behavior of the two algorithms are identical via the stated equivalences.
In particular, since f(x_k-1/2)=f(N^Tx̂_k-1/2),
the two algorithms terminate for the same value of k and the two
images returned by the respective algorithms can be obtained from
each other using x_k+1/2=N^Tx̂_k+1/2.










Next we show the equivalence of Algorithm <ref> and Algorithm
7 of <cit.> under some assumptions. To make the current paper
self-contained, we reproduce here (as Algorithm <ref>)
Algorithm 7 of <cit.>. We keep the boldface notation of <cit.>
in order to avoid confusion with the symbols used in the algorithms
introduced in the present paper. Consistently with the notation of
that previous paper, the function f' in Algorithm 7 is defined
by

    f'(𝐱)     =1/2𝐀𝐱-𝐲^2,

where 𝐀=RN^T and 𝐲=b.

As already defined in Subsection <ref>,
A=R^TR and y=R^Tb. By setting 𝐀
and 𝐲 in Algorithm 8 of <cit.> (restated
here as Algorithm <ref>) to be as we have just defined
them, we see that Algorithm <ref> is equivalent to
Algorithm 8 of <cit.> in the sense that if we have in the input
of the two algorithms the equivalences 𝐱=x̂, 𝐩=p̂,
and 𝐡=ĥ, then 𝐠=ĝ (to see this,
notice that in Algorithm <ref> we have Ny=NR^Tb=𝐀^T𝐲,
and NAN^T=NR^TRN^T=𝐀^T𝐀), the β
in the two algorithms are equal (since ẑ=ĝ) and, therefore,
we also have that the output values satisfy 𝐩=p̂,
𝐡=ĥ (because NAN^T=𝐀^T𝐀)
and 𝐱=x̂ (since, as is easily verified, the value
of α is the same for the two algorithms). Thus, for identical
inputs, Algorithms 4 and 6 will return identical outputs.







To show the equivalence of Algorithm <ref> and Algorithm
<ref>, we associate 𝐱_k,
𝐱_k+1/2, 𝐠_k,
𝐡_k, 𝐩_k in Algorithm <ref>
with x̂_k, x̂_k+1/2, ĝ_k,
ĥ_k, p̂_k in Algorithm <ref> by
𝐱_k=x̂_k, 𝐱_k+1/2=x̂_k+1/2,
𝐠_k=ĝ_k, 𝐡_k=ĥ_k,
𝐩_k=p̂_k. We now compare the step-by-step
executions of SupTPCG(x̂_0,a,γ,K,ε)
and SupCG(𝐱_0,a,γ,K,ε'), with
𝐱_0=x̂_0 and ε'=ε/2.
Clearly, after Step 1 in the two algorithms 𝐱_1/2=x̂_1/2.
That 𝐠_0=ĝ_0 follows from the facts that Ny=𝐀^T𝐲,
and NAN^T=𝐀^T𝐀. Continuing in this fashion
it is trivial to check that just before entering the 𝐰𝐡𝐢𝐥𝐞
statement for the first time in the two algorithms the values of the
associated vectors match as stated above.

We now consider the condition that appears in the 𝐰𝐡𝐢𝐥𝐞
statement in the algorithms. We claim that if just before entering
the 𝐰𝐡𝐢𝐥𝐞 in the two algorithms the values of the associated
vectors match as stated above, then the condition is either satisfied
in both algorithms or is not satisfied in both algorithm. Indeed,
due to the definition of f as the squared error and (<ref>),

    f(N^Tx̂_k-1/2)=‖ RN^Tx̂_k-1/2-b‖ ^2=‖𝐀𝐱_k-1/2-𝐲‖ ^2=2f'(𝐱_k-1/2).

Since ε'=ε/2, our claim on the satisfaction
of the conditions is valid.

We specify the perturbation operator in Step 9 of Algorithm <ref>
by

    perturbed(𝐱_k)=𝐱_k+N^-Ts_k,

where s_k is the same as given by Algorithm <ref>.
Using the above associations and facts, we see that the behavior of
Algorithm <ref> is equivalent to that of Algorithm <ref>,
and, therefore, to that of Algorithm 7 of <cit.>.

Theorem A.1 of <cit.> says that SupCG(𝐱_0,a,γ,K,ε')
terminates provided that ε'>ε_0', where

    ε_0'=min_𝐱1/2𝐀𝐱-𝐲^2.

Recalling the definitions of 𝐀 and 𝐲 after
(<ref>) and the fact that N is invertible, we get that
ε'_0=1/2ε_0, for

    ε_0=min_x̂‖ RN^Tx̂-b‖ ^2=min_x‖ Rx-b‖ ^2.

This together with the comparison of the the step-by-step executions
of SupTPCG(x̂_0,a,γ,K,ε) and SupCG(𝐱_0,a,γ,K,ε'),
with 𝐱_0=x̂_0 and ε'=ε/2,
gives us our promised termination result for SupPCG:

Given any positive number ε such that ε>ε_0,
with

    ε_0=min_x‖ Rx-b‖ ^2,

Algorithm <ref> terminates within a finite number of iterations.




§ CONCLUDING REMARKS




We have discussed the superiorized version of the Preconditioned Conjugate
Gradient method (SupPCG) for tomographic image reconstruction. Experimental
work has been presented indicating that, when compared to several
other methods, SupPCG produces images of good quality within a short
time. Furthermore, we have proved that the algorithm produces, in
a finite number of steps, an image whose data-inconsistency is no
worse than what is specified by the user, provided that such an image
exists.




§ ACKNOWLEDGMENTS


The research of E. S. Helou was funded by FAPESP grants 2013/07375-0
and 2016/24286-9. The research of M. V. W. Zibetti was partially supported
by CNPq grant 475553/2013-6. The research of C. Lin was partially
supported by China Scholarship Council. We are grateful
to Christoph Scharf for essential help with programming the PCG algorithm
in the SNARK14 environment.

10
chj17 Y. Censor, G. T. Herman, and M. Jiang. Superiorization: theory and applications. Inverse Problems, 33, 040301, 2017.
com96 P. L. Combettes. The convex feasibility problem in image recovery. Advances in Imaging and Electron Physics, 95:155–270, 1996.
gah14 E. Garduño and G. T. Herman. Superiorization of the ML-EM algorithm. IEEE Transactions on Nuclear Science, 61:162–172, 2014.
gh17 E. Garduño and G. T. Herman. Computerized tomography with total variation and with shearlets. Inverse Problems, 33, 044011, 2017.
hzm16 E. S. Helou, M. V. W. Zibetti, and E. X. Miqueles. Superiorization of incremental optimization algorithms for   statistical tomographic image reconstruction. Inverse Problems, 33, 044010, 2017.
HERMAN:2009a G. T. Herman. Fundamentals of Computerized Tomography: Image Reconstruction   from Projections. Springer, London, UK, 2nd. edition, 2009.
Herman2012 G. T. Herman, E. Garduño, R. Davidi, and Y. Censor. Superiorization: An optimization heuristic for medical physics. Medical Physics, 39:5532–5546, 2012.
hhl79 G. T. Herman, H. Hurwitz, A. Lent, and H.-P. Lung. On the Bayesian approach to image reconstruction. Information and Control, 42:60–71, 1979.
KLUK13 J. Klukowska, R. Davidi, and G. T. Herman. SNARK09 - A software package for the reconstruction of 2D   images from 1D projections. Computer Methods and Programs in Biomedicine, 110:424–440,   2013.
mhc98 R. Marabini, G. T. Herman, and J. M. Carazo. 3D reconstruction in electron microscopy using ART with smooth   spherically symmetric volume elements (blobs). Ultramicroscopy, 72:53–65, 1998.
Vogel-2002 C. R. Vogel. Computational Methods for Inverse Problems. SIAM, Philadelphia, PA, 2002.
zlh17 M. V. W. Zibetti, C. Lin, and G. T. Herman. Total variation superiorized conjugate gradient method for image   reconstruction. Inverse Problems, 34, 034001, 2018.

