# Accelerated Alternating Minimization for X-ray Tomographic Reconstruction

Peijian Ding This work was done while the author was an undergraduate student
at Emory University<EMAIL_ADDRESS>Advisor: Dr. James G. Nagy
Department of Mathematics, Emory University<EMAIL_ADDRESS>
(July, 2021)

###### Abstract

While Computerized Tomography (CT) images can help detect disease such as
Covid-19, regular CT machines are large and expensive. Cheaper and more
portable machines suffer from errors in geometry acquisition that downgrades
CT image quality. The errors in geometry can be represented with parameters in
the mathematical model for image reconstruction. To obtain a good image, we
formulate a nonlinear least squares problem that simultaneously reconstructs
the image and corrects for errors in the geometry parameters. We develop an
accelerated alternating minimization scheme to reconstruct the image and
geometry parameters.

###### keywords:

inverse problems, ill-posed problems, X-ray tomography, least squares problem,
medical imaging

##  1 Introduction

Tomography is a technique of displaying representations of a cross section
through an object through the use of some penetrating waves such as X-ray or
ultrasound. In simple words, it allows us to see the inside of an object
without breaking it. Thus, tomography is widely used in medical imaging,
seismic exploration, and material science. In medical imaging, a Computerized
Tomography (CT) Scan creates a cross-sectional image of human body by
combining X-ray images taken from different angles.

Figure 1: This is a simple illustration of our tomography problem.

Figure 2: $32\times 32$ true image of a Shepp–Logan phantom (left), image
computed by taking into account the correct source-to-object distance
(middle), image computed by using incorrect source-to-object distance (right)

During a CT scan, the patient lies on a bed that slowly moves through the
gantry while the X-ray tube rotates around the patient and shoots X-ray beams
through the human body, received by a detector. See Figure 1. Then, an image
of the cross section of the human body is reconstructed following a
mathematical procedure.

Although CT images can help doctors diagnose and monitor diseases such as
Covid-19, regular CT machines are heavy and expensive and not widely available
in less developed areas. The goal of this paper is to compensate for cheaper
and more portable machines by solving for geometry parameters such as the
source-to-object distance that may not be calibrated precisely during the
imaging process.

Source-to-object distance measures how far away the center of the object is
from the X-ray source. Since the source-to-object distance may vary from angle
to angle, the reconstructed image will be corrupted if incorrect values are
used, as illustrated in Figure 2. Similar conclusions can be drawn if
incorrect angles are used to reconstruct the image. Thus, our algorithm will
significantly improve image quality by taking into account the variation in
geometry parameters.

For the purpose of this paper, we primarily focus on 2D computed tomography to
simplify notations and to reduce computational costs required to perform a
substantial number of numerical experiments. But the methods discussed in this
paper can be extended to $3$D problems. An X-ray imaging problem can be
formulated as an inverse problem in linear algebra:

$Ax=b,$

where $A$ is a $m\times n$ forward operator, $b$ is projection data, also
referred to as a sinogram, and $x$ is the solution vector which represents the
image.

Since the forward operator $A$ is determined by geometry parameters $r$, we
form the following nonlinear least squares problem:

(1) $\underset{x,r}{\arg\min}||A(r)x-b||_{2}^{2}$

where $x$ represents the image vector, and $A$ is the forward operator that is
a function of $r$ and maps the true image $x$ to the sinogram $b$. In this
paper, we denote $r$ as a general geometry parameter vector that can represent
both source-to-object distance $d$ and errors in angles $\delta\theta$. Thus,
the conclusions we make about $r$ also apply to both $d$ and $\delta\theta$.
In this paper, we show that Equation 1 can be solved using an accelerated
Block Coordinate Descent method. In addition to incorporating an acceleration
scheme, we also exploit separability to further reduce the computational cost.

##  2 Alternating Minimization Scheme

In order to solve the nonlinear least squares problem, an intuitive approach
to consider is the Block Coordinate Descent (BCD) algorithm. In this section,
we discuss the linear and nonlinear least squares problem involved in BCD. We
also briefly discuss the variable projection approach and argue the advantage
of BCD over variable projection for our problem.

###  2.1 Block Coordinate Descent

Block Coordinate Descent is a simple approach to solve an optimization
problem. Its idea is based on the general Coordinate Descent (CD) algorithm.
Because of its lack of sophistication, most optimization researchers have not
focused on this approach until recently when CD approaches were found to be
computationally competitive to other reputable alternatives in various
applications such as machine learning [15]. Since we are able to exploit
separability in the geometry parameters, the BCD method, which is given in
Algorithm 1, is worth investigating for tomographic reconstruction.

1: Input: $r_{0}\in\mathbb{R}^{N_{A}}$, $x_{0}\in\mathbb{R}^{n}$

2: for $k=1,2,\dots$ until a stopping criterion holds do

3: $r_{k}=\underset{r}{\arg\min}||A(r)x_{k-1}-b||_{2}^{2}$

4: $x_{k}=\underset{x}{\arg\min}||A(r_{k})x-b||_{2}^{2}$

5: end for

Algorithm 1 BCD to Reconstruct Geometry and Image Parameters

$N_{A}$ denotes the number of angles, and $n$ is the length of the image
vector. Note that in practice an initial estimate $r_{0}$ is given. With this
information, we can easily obtain $x_{0}$ by solving the linear least squares
problem,

$x_{0}=\underset{x}{\arg\min}||A(r_{0})x-b||_{2}^{2}.$

We remark that for ill-posed problems, computing $x_{0}$ and Step $4$ of
Algorithm 1 typically requires incorporating regularization procedures to
avoid amplifying noise when solving the linear least squares problems. This is
discussed further in section 2.2.1. There is another property in our matrix
that makes this alternating minimization scheme favorable. The matrix $A$ and
vector $b$ can be partitioned into block of rows as

$A=\begin{bmatrix}A_{1}(r^{(1)})\\\ \vdots\\\ A_{i}(r^{(i)})\\\ \vdots\\\
A_{N_{A}}(r^{(N_{A})})\end{bmatrix}\quad b=\left[\begin{array}[]{cc}b_{1}\\\
\vdots\\\ b_{i}\\\ \vdots\\\ b_{N_{A}}\end{array}\right]$

where errors in geometry parameters associated with rows in each
$[A_{i}(r^{(i)}),b_{i}]$ block are assumed constant. For example, if we
collect projections at $0,1,2,\dots,359$ degrees around the center of the
object, and $N_{A}=10$, then an error is introduced into the geometry
parameters once every $36$ degrees. In reality, $N_{A}=360$, but if some of
the errors are small enough to ignore, using a smaller $N_{A}$ can reduce the
computational cost. The choice of $N_{A}$ is discussed further in Section 4.2.

The separability of the geometry parameters allows us to solve for a set of
much smaller systems at Step $3$ in Algorithm 1, namely,

$r_{k}^{(i)}=\underset{r}{\arg\min}||A_{i}(r)x_{k-1}-b_{i}||_{2}^{2},\quad
i=1,2,\dots,N_{A}$

and then

$r_{k}=\left[\begin{array}[]{cc}r_{k}^{(1)}\\\ \vdots\\\
r_{k}^{(N_{A})}\end{array}\right].$

Note that parameters in the vector $r_{k}$ are completely independent, so they
can be updated simultaneously on a parallel computing architecture, which
could dramatically lower the computing time. In fact, if we only consider
source-to-object distance as geometry parameters, the dimension of the
parameter $r^{(i)}$ is one. Comparing to this alternating minimization scheme,
Variable Projection [9][10] cannot utilize this matrix property because its
matrix multiplication process would take away the separability property. Thus,
it makes this alternating minimization scheme more worthwhile to investigate.

Next, we discuss linear least squares solvers and nonlinear least squares
solvers respectively.

###  2.2 Linear Least Squares Problem

In this subsection, we consider the linear least squares problem in Step $4$
of Algorithm 1, $\underset{x}{\min}||Ax-b||_{2}^{2}$, where we assume matrix
$A$ is fixed, e.g. $A=A(r_{k})$.

####  2.2.1 Regularization

Typically in inverse problems, the data we obtain is not the exact data. This
is also the case in our X-ray tomography problem, even in the case when
geometry parameters, and hence matrix $A$, are known exactly. Ideally, we want
to solve $Ax_{true}=b_{true}$ but the linear system that we actually solve is:

(2) $Ax=b=b_{true}+\eta$

where $\eta$ is the noise in our measurement data, $b_{true}$ is the noise-
free data, $A\in\mathbb{R}^{m\times n}$, and $b\in\mathbb{R}^{m}$. For
simplicity, we assume $A$ is full rank and $m=n$. Our discussion also holds
for more general cases when we compute a pseudo-inverse instead of an exact
inverse.

By using the SVD of $A=U\Sigma V^{T}$ where $\sigma_{i}$ is the $i^{th}$
singular value and $u_{i}$, $v_{i}$ are the $i^{th}$ column of the left and
right orthogonal matrix $U$ and $V$, we first recall that
$Av_{i}=\sigma_{i}u_{i}$, and we note that we can find scalars $\alpha_{i}$
and $\eta_{i}$ such that

$x_{true}=\sum_{i=1}^{n}\alpha_{i}v_{i}$

and  $\eta=\sum_{i=1}^{n}\eta_{i}u_{i}$

Using these relations, we obtain

(3)
$\begin{split}b&=\sum_{i=1}^{n}(\alpha_{i}\sigma_{i}+\eta_{i})u_{i}\end{split}$

Then, notice that

(4) $\begin{split}A^{-1}u_{i}&=V\Sigma^{-1}U^{T}u_{i}\\\
&=\frac{1}{\sigma_{i}}v_{i}\end{split}$

Thus,

(5) $\begin{split}x&=A^{-1}b\\\
&=\sum_{i=1}^{n}\left(\alpha_{i}+\frac{\eta_{i}}{\sigma_{i}}\right)v_{i}\end{split}$

From this result we observe that the noise in the data is magnified by the
small singular values. Since
$\sigma_{1}\geqslant\sigma_{2}\geqslant...\geqslant\sigma_{n}\geqslant 0$,
larger indices correspond to smaller singular values. The computed solution is
dominated by noise amplified by division of small singular values. Thus, we
need regularization schemes to filter out this noise. In particular, a
regularized solution can be written as:

(6)
$x_{reg}=\sum_{i=0}^{n}\phi_{i}\left(\alpha_{i}+\frac{\eta_{i}}{\sigma_{i}}\right)v_{i}$

where the scalar $0\leqslant\phi_{i}\leqslant 1$ is called a filter factor. As
$\sigma_{i}$ decreases, the filter factors should approach zero so that the
noise contributed by the small singular values are filtered out.

##### Truncated Singular Value Decomposition

Since the noise is magnified by small singular values, the most intuitive
approach is to cut off the small singular values by setting them to zero. This
is called Truncated SVD regularization. The TSVD solution to the inverse
problem is given by

(7) $x_{reg}=\sum_{i=1}^{k}\frac{b^{T}u_{i}}{\sigma_{i}}v_{i}$

where $k\leqslant n$. The critical part in the TSVD is identifying the
threshold $k$. One approach is choosing $k$ at a significant drop-off of
singular values, as illustrated by Figure 3.

Figure 3: The singular values plot of a $10\times 10$ diagonal matrix whose
singular values are respectively $99$, $98$ … $91$, and $1$

In this case, $k=9$ can be easily identified as the threshold for the TSVD
approach. However, typically in inverse problems, singular values decay
smoothly, as illustrated by Figure 4.

Figure 4: The singular values of the $256\times 256$ test problem _heat_
generated from Regularization Tools in MATLAB

In cases like Figure 4, a reliable cut-off threshold for singular values is
hard to find. Note that $k=250$ is not a good choice for the cut-off because
the singular values are very close to zero for smaller indices $k$ (e.g.
$\sigma_{200}\approx 10^{-8}$). Therefore, we need more advanced
regularization techniques to deal with smoothly decaying singular values.

##### Tikhonov Regularization

Classical Tikhonov Regularization can be used to solve ill-posed problems. The
regularized solution $x_{reg}$ is the unique solution to the following:

(8) $\underset{x}{\min}{||Ax-b||_{2}^{2}+\lambda^{2}||x||_{2}^{2}}$

where $\lambda$ is the regularization parameter that controls the smoothness
of the regularized solution. The above equation is equivalent to the
following:

(9) $\underset{x}{\min}\left\|\begin{bmatrix}A\\\ \lambda
I\end{bmatrix}x-\begin{bmatrix}b\\\ 0\end{bmatrix}\right\|_{2}^{2}$

Then the normal equations for this least squares problem can be written as:

(10) $(A^{T}A+\lambda^{2}I)x_{reg}=A^{T}b.$

From the normal equations, we obtain the following:

(11) $\begin{split}x_{reg}&=(A^{T}A+\lambda^{2}I)^{-1}A^{T}b\\\
&=\sum_{i=1}^{n}\phi_{i}\frac{b^{T}u_{i}}{\sigma_{i}}v_{i}\end{split}$

where the filter factor is
$\phi_{i}=\frac{\sigma_{i}^{2}}{\sigma_{i}^{2}+\lambda^{2}}$. For simplicity,
we assumed $A$ has full rank. The pseudo-inverse can be used for the rank
deficient case and Equation 11 would sum from $i=1$ to $rank(A)$ instead of
the full column rank $n$. The modified Tikhonov regularization is not too much
different from the classical Tikhonov except that it uses the 2-norm of $Lx$
instead of $x$ in Equation 8, where $L$ is a $p\times n$ matrix with
$p\leqslant n$.

$\underset{x}{\min}{||Ax-b||_{2}^{2}+\lambda^{2}||Lx||_{2}^{2}}$

In contrast to TSVD where singular values after $\sigma_{k}$ are cut off, the
Tikhonov regularization applies a smoother filter to all the singular values.
Given a good regularization parameter, the large singular values would not be
affected too much whereas the small values would be gradually filtered more as
they approach zero. The quality of the filtering depends on how we choose the
regularization parameter. If $\lambda=0$, then all $\phi_{i}=1$ and we are
directly calculating the inverse (or pseudo-inverse) solution. If we select a
very large $\lambda\gg\sigma_{1}$, all $\phi_{i}$ would approach zero and we
would have over-smoothed the solution.

####  2.2.2 Parameter Choice Methods

The choice of regularization parameter is critical to the quality of the
regularized solution. Parameter choice methods can usually be divided into two
classes depending on their assumption about error norm
$||\eta||_{2}^{2}=||b-b_{true}||_{2}^{2}$ where $b$ is the measured data and
$b_{true}$ is the noise free data [6]. The first class contains methods based
on a good estimate of $||\eta||_{2}^{2}$. The second class includes methods
that are not based on a good estimate of $||\eta||_{2}^{2}$, but seek to
extract this information from the given right hand side $b$. In this section,
we introduce the Generalized Cross-Validation method (GCV) which is a popular
method in the second class. The underlying idea of GCV is that a good
regularization parameter should predict missing values. For example, if some
data point $b_{i}$ is missing in the right hand side, the regularized solution
should predict this missing value well. GCV can be written as:

(12)
$G(\lambda)=\frac{n||Ax_{reg}-b||_{2}^{2}}{trace(I_{m}-AA^{\dagger}_{F})^{2}}$

where $x_{reg}$ is the regularized solution,
$A^{\dagger}_{F}=\sum_{i=1}^{n}\phi_{i}\frac{v_{i}u_{i}^{T}}{\sigma_{i}}=V\Sigma^{\dagger}_{F}U^{T}$,
and
$\Sigma^{\dagger}_{F}=diag(\frac{\phi_{1}}{\sigma_{1}},\frac{\phi_{2}}{\sigma_{2}},...,\frac{\phi_{n}}{\sigma_{n}})$.
The goal is to find $\lambda$ such that $G(\lambda)$ is minimized. Now we
simplify the numerator and the denominator of $G(\lambda)$ respectively:

$\begin{split}||Ax_{reg}-b||_{2}^{2}&=||\Sigma
V^{T}x_{reg}-U^{T}b||_{2}^{2}\\\
&=||(\Sigma\Sigma^{\dagger}_{F}-I)U^{T}b||_{2}^{2}\end{split}$

$\begin{split}trace(I_{m}-AA^{\dagger}_{F})&=trace((I_{m}-U\Sigma\Sigma_{F}^{\dagger}U^{T}))\\\
&=trace(I_{m}-\Sigma\Sigma_{F}^{\dagger})\end{split}$

Thus, $G(\lambda)$ becomes:

(13)
$G(\lambda)=\frac{n||(\Sigma\Sigma^{\dagger}_{F}-I)U^{T}b||_{2}^{2}}{trace(I_{m}-\Sigma\Sigma_{F}^{\dagger})}$

Since $\phi_{i}=\frac{\sigma_{i}^{2}}{\sigma_{i}^{2}+\lambda^{2}}$ in the
Tikhonov case,

$\Sigma_{F}^{\dagger}=diag\left(\frac{\phi_{1}}{\sigma_{1}},\frac{\phi_{2}}{\sigma_{2}},...,\frac{\phi_{n}}{\sigma_{n}}\right)=diag\left(\frac{\sigma_{1}}{\sigma_{1}^{2}+\lambda^{2}},\frac{\sigma_{2}}{\sigma_{2}^{2}+\lambda^{2}},...,\frac{\sigma_{n}}{\sigma_{n}^{2}+\lambda^{2}}\right)$

Then we can write $G(\lambda)$ in the case of Tikhonov regularization as:

(14)
$G(\lambda)=\frac{n\sum_{i=1}^{n}(\frac{\hat{b_{i}}}{\sigma_{i}^{2}+\lambda^{2}})^{2}}{(\sum_{i=1}^{n}\frac{1}{\sigma_{i}^{2}+\lambda^{2}})^{2}}$

where $\hat{b}=U^{T}b$. We can solve for $\lambda$ by using the function
fminbnd in MATLAB which is based on golden section search and parabolic
interpolation. In practical applications, several studies have found that
occasionally GCV would drastically under-smooth the solution by choosing the
regularization parameter too small [3] [4] [13]. In [1], GCV has been found to
over-smooth the solution in the Lanczos-hybrid methods, discussed in the next
subsection. To alleviate this difficulty, it was proposed to use the weighted
GCV method:

(15)
$G(w,\lambda)=\frac{n||Ax_{reg}-b||_{2}^{2}}{trace(I_{m}-wAA^{\dagger}_{F})^{2}}$

where $w$ is the weight parameter that determines the function $G(w,\lambda)$
along with $\lambda$. When $w=1$, we have the non-weighted version of GCV like
in Equation 13. When $w>1$, the solution is smoother. When $w<1$, the solution
is less smooth. So far in weighted-GCV literature, only experimental
approaches have been used to determine the value for $w$.

####  2.2.3 Hybrid LSQR

We have so far discussed regularization and parameter choice methods using the
SVD. For large-scale problems, such as in image reconstruction, directly
applying SVD based methods is not computationally feasible. In this subsection
we describe a hybrid LSQR scheme that combines an efficient iterative method
with SVD based approaches that enforce regularization on small projected sub-
problems.

The standard LSQR algorithm projects the linear least squares problem onto a
sequence of Krylov subspace of small and increasing dimensions. The best
approximation $x_{k}$ to the least squares problem in the Krylov subspace
$K_{k}$ is given by $x_{k}=V_{k}y_{k}$, where $V_{k}$ is a $n\times k$ matrix
with orthonormal columns at the $k^{th}$ step of the Golub-Kahan
bidiagonalization process [11]. $y_{k}$ can be solved by

$y_{k}=\arg\min_{y_{k}}||B_{k}y_{k}-\beta_{1}e_{1}||_{2}$

where $B_{k}$ is the bidiagonal matrix at the $k^{th}$ step of the Golub-Kahan
bidiagonalization process, $\beta$ is a scalar, and $e_{1}$ is the standard
basis vector. When being applied to ill-posed problems, LSQR exhibits a semi-
convergence behavior which means that early iterations construct information
related to the solution while later iterations construct information related
to noise [1].

This can be compensated by applying a direct regularization method such as
Tikhonov or TSVD, which can be solved cheaply on a small scale problem of the
reduced linear least squares in the Krylov subspace. So, we can write the
Hybrid LSQR using Tikhonov regularization as:

$\underset{y_{k}}{\min}||B_{k}y_{k}-\beta_{1}e_{1}||_{2}^{2}+\lambda_{k}^{2}||y_{k}||_{2}^{2}$

where $\lambda_{k}$ is a regularization parameter chosen at the $k^{th}$
iteration using the weighted GCV method discussed in Section 2.2.2. A method
like GCV can be used to choose a stopping iteration so that $k$ will not be
too large; details can be found in [1].

Comparing to LSQR, this hybrid method can effectively stabilize the iterations
[1]. Although at each iteration a new regularization parameter must be chosen,
it is not computationally expensive for the projected problem.

To summarize the method, the hybrid LSQR method projects the large scale
linear least squares problem onto a low-dimensional Krylov subspace where we
can inexpensively apply a direct regularization method like the adaptive
weighted-GCV.

###  2.3 Nonlinear Least Squares

In our alternating minimization scheme, we iteratively solve the image and the
geometry parameters. While we have discussed methods to solve the linear least
squares problem in the previous section, we need other tools to solve the
nonlinear least squares problem:

(16) $\underset{r}{\min}||A(r)x-b||_{2}$

where $x$ is approximated by the linear least squares solution we obtained by
using hybrid LSQR in Section 2.3.

We utilize the implicit filtering method which solves the bound-constraint
optimization problem for which the derivative information is not available
[8]. Since we do not have the derivative information of our objective function
and a reasonable bound can be established for the geometry parameters in our
tomographic reconstruction problem, implicit filtering serves as a good tool
to solve our problem. Implicit filtering builds the local model of the
objective function using a quasi-Newton method.

In our numerical experiments, we compare implicit filtering to the MATLAB
function fminbnd.

##  3 Acceleration Algorithms

In the previous section, we have introduced the alternating minimization
scheme and the methods we use to solve least squares problems. In this
section, we introduce methods that will accelerate the convergence of our
minimization scheme.

###  3.1 Accelerated Block Coordinate Descent

Since we can divide variables in our least squares problem into two blocks –
geometry parameters $r$ and image $x$, it makes sense for us to directly
investigate methods that accelerate the BCD algorithm. We implemented
Accelerated Block Coordinate Descent (ABCD). This method can be applied to a
four-block problem by dividing it into two larger blocks, but in our problem
we do not have to do so. Theoretically, the proposed acceleration method in
[2] has a complexity of $O(\frac{1}{k^{2}})$. In our implementation, we
simplify the algorithm as Algorithm 2.

Algorithm 2 Accelerated Block Coordinate Descent

1: Inputs: $t_{0}=1$, $r_{0}\in\mathbb{R}^{N_{A}}$ and
$x_{0}\in\mathbb{R}^{n}$

2: for $k=1,2,\dots$ until a stopping criterion holds do

3: $\tilde{r}_{k}=\underset{r}{\arg\min}||A(r)x_{k-1}-b||_{2}^{2}$
$\tilde{x}_{k}=\underset{x}{\arg\min}||A(\tilde{r}_{k})x-b||_{2}^{2}+\lambda^{2}||x||_{2}^{2}$

4: $\tilde{w}_{k}=(\tilde{x}_{k},\tilde{r}_{k})$
$t_{k}=\frac{1}{2}(1+\sqrt{1+4t^{2}_{k-1}})$
$w_{k}=\tilde{w}_{k-1}+\frac{t_{k-1}}{t_{k+1}}(\tilde{w}_{k}-\tilde{w}_{k-1})$,
where $w_{k}=(x_{k},r_{k})$

5: end for

$N_{A}$ denotes the number of angles, and $n$ is the length of the image
vector. We remark that we can exploit separability when solving for
$\tilde{r}_{k}$, as discussed in Section 2.1. As mentioned in Algorithm 1,
since $r_{0}$ is given,

$x_{0}=\underset{x}{\arg\min}||A(r_{0})x-b||_{2}^{2}+\lambda^{2}||x||_{2}^{2}.$

The Tikhonov regularized least squares problems for $x_{0}$ and
$\tilde{x}_{k}$ are solved using the hybrid scheme IRhybrid_lsqr provided in
IRTools[5] in MATLAB. The nonlinear least squares problem for $\tilde{r_{k}}$
is solved using a MATLAB package called imfil [8]. We have found that
accelerating the solution vector $x$ alone has yielded stabler results with
slightly better accuracy than performing acceleration on both $x$ and $r$. We
show the numerical experiments in section 4.

###  3.2 Anderson Acceleration

Anderson Acceleration, also called Anderson mixing, is a method used to
accelerate the convergence of fixed point iteration. Note that we can write
out alternating minimization scheme as a fixed point iteration,

Algorithm 3 Fixed point iteration of the image vector

1: Input: $x_{k}\in\mathbb{R}^{n}$, Output: $x_{k+1}=g(x_{k})$

2: $r_{k+1}=\underset{r}{\arg\min}||A(r)x_{k}-b||_{2}^{2}$

3:
$x_{k+1}=\underset{x}{\arg\min}||A(r_{k+1})x-b||_{2}^{2}+\lambda^{2}||x||_{2}^{2}$

where $g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is the fixed point iteration
of image vector $x$, as shown in Algorithm 3.

For this fixed point iteration, the general form of Anderson Acceleration is
formed as the following:

Algorithm 4 Anderson Acceleration

1: Inputs: $x_{0}$ and $m\geqslant 1$

2: Set $x_{1}=g(x_{0})$, using Algorithm 3

3: for $k=1,2,\dots$ until a stopping criterion holds do

4: $m_{k}=min(m,k$)

5: Set $F_{k}=(f_{k-m_{k}},...,f_{k})$, where $f_{i}=g_{i}(x_{i})-x_{i}$ and
$g_{i}(x_{i})$ comes from Algorithm 3

6: Determine $\alpha^{(k)}=(\alpha_{0}^{(k)},...,\alpha_{m_{k}}^{(k)})^{T}$
that solves $\underset{\alpha}{\min}||F_{k}\alpha||_{2}$ s.t
$\sum_{i=0}^{m_{k}}\alpha_{i}=1$

7: Set $x_{k+1}=\sum_{i=0}^{m_{k}}\alpha_{i}^{(k)}g(x_{k-m_{k}+i})$, where
$g(x_{k-m_{k}+i})$ is from Algorithm 3

8: end for

We can cast the linear constrained optimization problem in Step $7$ of
Algorithm 4 into an unconstrained form which is straightforward to solve and
convenient for efficient implementation [14].

We define $\nabla f_{i}=f_{i+1}-f_{i}$ for each $i$ and set $\nabla F=(\nabla
f_{k-m_{k}},...,\nabla f_{k-1})$. Then the least squares problem is equivalent
to

$\underset{\gamma=(\gamma_{0},...\gamma_{m_{k}-1})^{T}}{min}||f_{k}-F_{k}\gamma||_{2},$

where $\alpha_{0}=\gamma_{0}$ and $\alpha_{i}=\gamma_{i}-\gamma_{i-1}$, for
$1\leqslant i\leqslant m_{k}-1$ and $\alpha_{m_{k}}=1-\gamma_{m_{k}-1}$.

This unconstrained least squares problem leads to a modified version of
Anderson Acceleration in Algorithm 5,

Algorithm 5 Modified Anderson Acceleration

1: Given $x_{0}$ and $m\geqslant 1$

2: Set $x_{1}=g(x_{0})$, using Algorithm 3

3: for k=1,2,… do

4: $m_{k}=min(m,k$)

5: Determine $\gamma^{(k)}=(\gamma_{0}^{(k)},...,\gamma_{m_{k}-1}^{(k)})^{T}$
that solves
$\underset{\gamma=(\gamma_{0},...\gamma_{m_{k}-1})^{T}}{min}||f_{k}-F_{k}\gamma||_{2}$

6: Set $x_{k+1}=g(x_{k})-G_{k}\gamma^{(k)}$, where $g(x_{k})$ comes from
Algorithm 3

7: end for

where

$x_{k+1}=g(x_{k})-\sum_{i=1}^{m_{k}-1}\gamma_{i}^{(k)}[g(x_{k-m_{k}+i+1})-g(x_{k-m_{k}+i})]=g(x_{k})-G_{k}\gamma^{(k)}$

with $G=(\nabla g_{k-m_{k}},...,\nabla g_{k-1})$, $\nabla
g_{i}=g(x_{i+1})-g(x_{i})$.

Homer Walker proposed implementation that efficiently updates the QR factors
in the decomposition $F_{k}=Q_{k}R_{k}$ [14]. The basic logic is the
following: every $F_{k}$ is obtained from $F_{k-1}$ with a column added on the
right. If the resulting matrix has more columns than $m$, then delete one from
the left. The column addition can be achieved by a modified Gram–Schmidt
process. The deletion process is a little more complicated. We delete the
first column on the left when $m_{k-1}=m$. If $F_{k-1}=QR$, then
$F_{k-1}(:,2:m)=QR(:,2:m)$, where $R(:,2:m)$ is upper-Hessenberg. Then, we can
determine $m$ Givens rotations to cancel out the entries in the sub-diagonal.

##  4 Numerical Experiments

In this section, we make a few comparisons of different methods to solve the
X-ray tomography problem. Firstly, we compare the speed of BCD exploiting the
separability of geometry and the speed of BCD without using such property.
Secondly, we compare results produced from different number of angles.
Thirdly, we make comparisons of all the acceleration schemes. Fourthly,
different regularization parameters in the linear least squares solvers are
compared. Fifthly, we compare implicit filtering with the MATLAB function
fminbnd as the nonlinear least squares solver. Lastly, we show our algorithm
works for both geometry parameters $d$ and $\delta\theta$. For experiments
1-5, we solve problems with only unknown source-to-object distance $d$, where
geometry parameters $r=d$. For the last experiment, $r=(d,\delta\theta)$.
Comparisons in all experiments are made about geometry errors and
reconstruction errors. Geometry errors are defined as the relative errors of
geometry parameters, $\frac{||r-r_{true}||_{2}}{||r_{true}||_{2}}$. For
experiment 1-5, where $r=d$, geometry errors are
$\frac{||d-d_{true}||_{2}}{||d_{true}||_{2}}$. For experiment 6, where
$r=(d,\delta\theta)$ geometry errors are represented by both
$\frac{||d-d_{true}||_{2}}{||d_{true}||_{2}}$ and
$\frac{||\delta\theta-\delta\theta_{true}||_{2}}{||\delta\theta_{true}||_{2}}$.
The reconstruction errors are defined as the relative errors of the image,
$\frac{||x-x_{true}||_{2}}{||x_{true}||_{2}}$.

We use fan-beam projection for all our tomography problems for the sake of
consistency. Note that we can also easily adapt our code to solve parallel
beam projection problems by using the IRtools [5] and AIR Tools [7] MATLAB
packages.

In practice, good initial guesses of geometry parameters $r$ are available and
prior knowledge can help us set the bounds for them. For the source-to-object
distance $d$, we generate a test problem where true $d$ values, $d_{true}$,
are random numbers (chosen from a uniform distribution) between $1.5$ and
$2.5$. We use a constant initial guess of $d=2$ for all angles and set the
bounds for $d$ from $1.5$ to $2.5$. For the errors in angles $\delta\theta$,
we assume for each set of angles the error bound is from $-0.5$ to $0.5$. The
test image we use is Shepp-Logan Phantom [12], with image size $n=32\times
32$. The noise level is $\frac{||\eta||_{2}}{||b||_{2}}=0.01$, where $\eta$ is
a vector with random entries chosen from a normal distribution. Budget is a
hyper-parameter in imfil that stands for the maximum number of function
evaluations in the nonlinear least squares solver. $N_{A}$ stands for number
of angles. Moreover, $0^{th}$ iteration is included in the relative errors
figures below. This represents the relative error of the initial guesses with
regard to the true solution.

###  4.1 BCD Exploiting Separability vs BCD

Let BCD that exploits separability be called BCDS. In this section, we compare
the running time of BCD and that of BCDS. In this subsection’s numerical
experiments, $budget=1000$ (used in imfil to put a limit on the number of
function evaluations), and $N_{A}=10$. As we see in Figure 5, BCDS
dramatically speeds up the convergence because the separability allows us to
solve a much smaller problem independently for one $r(i)$ at a time. The
average running time of image reconstruction using BCD is $12.8s$, around the
same as BCDS’s $12.6s$. The time of geometry reconstruction using BCD is
$1857.5s$, more then ten times longer than BCDS’s $157.1s$. This is also the
reason BCDS is discussed first in this section. In the remaining experiments,
we always use the BCDS to reduce the running time. Moreover, the geometry
errors and reconstruction errors of BCDS are both better than that of BCD. In
Figure 6, the phantom reconstructed by BCDS is much closer to the true image
than the one from BCD.

Figure 5: Comparison between BCD and BCDS.

Figure 6: Comparison of $32\times 32$ Shepp-Logan phantom: true image (left),
BCDS image (middle), BCD image (right).

A typical CT image using fan beam x-ray sources collects data at angles of one
degree increment, from 0 to 359 degrees. In a perfect machine, the geometry
parameters are precisely known each time the source is rotated to a new angle.
In our experimental scenario, we assume the geometry parameters are only known
approximately. To experiment with various scenarios, we assume that errors are
introduced into the geometry parameters once every $\frac{360}{N_{A}}$
degrees. That is, with $N_{A}=10$, errors occur once every 36 degrees.

Partially, $N_{A}$ depends on the precision of machine calibration. For a good
machine, $N_{A}$ may be small. Also, $N_{A}$ may depend on how precise we
measure the data. For example, two geometry parameters that are different in
terms of double precision may be rounded to the same number in single
precision. In this scenario, the difference could be small enough that we can
treat the two geometry parameters as being equal.

If the number of angles is larger than the true number of angles, we may end
up solving a larger problem than we need. For example, if for every $36$
degrees only one geometry error is introduced, then $N_{A_{true}}=10$. If we
assume $N_{A}=20$, the average of the two geometry errors per $36$ degrees
would approximate the one true geometry error introduced in that set of
angles. If the number of angles is smaller than the true number of angles, we
end up solving a low dimension approximation. We do not seek to solve for the
image exactly but aim to compute good approximations that yield much better
results than not considering geometry parameters at all. In practice, we can
consider the number of angles as a hyper-parameter that practitioners can set
based on their expertise and knowledge of the machine calibration. In our
problems, we assume to know the number of angles.

###  4.2 Number of Angles

We compare the case where the number of angles is $5$, $10$, and $20$
respectively. We want to explore the differences in relative error of $r$,
relative error of $x$, the convergence, and the image quality.

In this comparison, $budget=10$. The reason that we can use a such small
budget is we are essentially solving one scalar of $r$ at a time. The
recommended formula [8] to set up a good guess for budget is $10*N^{2}$, where
$N$ is the length of the solution vector. Since we solve for each component of
the solution vector separately, the budget is just $10$ based on the
recommended formula. Adjusting budget according to each particular problem as
a hyper-parameter could further improve the quality of the result. Since we
try to make apple-to-apple comparisons, we keep budget the same for all three
cases.

h

Figure 7: Comparison of number of angles: geometry errors (left),
reconstruction errors (right)

Both geometry and image parameters converge for all three different number of
angles, as shown in Figure 7. Although the errors when $N_{A}=10$ are not as
low as when $N_{A}=5$, the reconstruction errors when $N_{A}=10$ still have a
$40\%$ reduction comparing to its initial error level. When $N_{A}=20$, the
geometry error has a huge spike at the first iteration. Although our algorithm
has gradually decreased the geometry error after the first iteration, the
resulting geometry error is almost the same as the geometry error of the
initial guess. The reconstruction error drops from more than $100\%$ at the
initial guess to below $85\%$, but the error climbed back up in the next few
iterations. In the end, there was a $10\%$ drop in reconstruction errors
comparing to the initial guess. There is one possible explanation for this
unusual phenomenon when $N_{A}=20$. Let $r_{k}$ and $x_{k}$ represent the
value of $r$ and $x$ at $k^{th}$ iteration in BCDS. First, given $x_{0}$,
$r_{1}$ is calculated from the nonlinear least squares problem in Algorithm 1.
$r_{1}$ is particularly bad because how bad $x_{0}$ is. Despite this $x_{0}$,
BCDS is still improving the geometry parameter after the first iteration.

###  4.3 Acceleration

In this section, we compare BCDS, ABCDS, and BCDS with Anderson acceleration,
where ABCDS represents the accelerated block coordinate descent that exploits
separability of geometry parameters. The number of columns of the matrix
$F_{k}$, also referred to as memory size in this paper, in Algorithm 5 is
$m=5$. Since the linear least squares problem at Step $5$ of Algorithm 5 is
relatively small, we directly use MATLAB backslash operator to solve it.

There are two implementations of ABCDS that we compare. Then, we compare
Anderson acceleration, ABCDS, and the BCDS method.

####  4.3.1 ABCDS

Since we can choose to apply the acceleration scheme on the image vector $x$
only or apply it on both geometry parameter $r$ and image vector $x$, we
compare the two versions of ABCDS. We denote the former version of ABCDS as
ABCDS-1. We denote the latter version that applies the acceleration on both
image and geometry parameters as ABCDS-b. The experiment setup for this
subsection is $budget=100$ and $N_{A}=10$.

Figure 8: Comparison between ABCDS-1 and ABCDS-b: geometry errors (left),
reconstruction errors (right)

From Figure 8, we can tell that although both approaches converge and produce
similar results at the end of iterations, the geometry errors of ABCDS-1 are
slightly lower and more stable than that of ABCDS-b. Thus, applying the
acceleration scheme on the image vector $x$ alone is enough to produce stable
and convergent results. For the sake of simplicity, when we mention ABCDS in
the rest of the paper, it always refers to ABCDS-1 .

####  4.3.2 BCDS, ABCDS, and Anderson

Figure 9: Comparison of BCDS, ABCDS, and Anderson: geometry errors (left),
reconstruction errors (right)

In Figure 9, geometry errors and reconstruction errors converge for all three
methods. The geometry errors and the reconstruction errors of BCDS converge
around $15^{th}$ iteration. The both geometry and reconstruction errors of
ABCDS converge around $10^{th}$ iteration. Anderson acceleration converges
slowest because there are some spikes around $15^{th}$ iteration for both
geometry and reconstruction errors. There are safeguards that can ensure
global convergence of the type I Anderson acceleration [16]. However, in our
problem, Anderson acceleration does converge but does not accelerate the
convergence. This could be an interesting challenge to be explored in the
future. For our problem, ABCDS is a better acceleration scheme its convergence
and acceleration effects.

###  4.4 Regularization

In this section, we compare BCDS with different regularization parameters: no
regularization, GCV, and weighted-GCV, as shown in Figure 10. The reason we
use BCDS without the acceleration methods is that we want to see the direct
impact of regularization on the alternating minimization scheme. In this
experiment, $N_{A}=10$ and $budget=100$.

Figure 10: Comparison of regularization: geometry errors (left),
reconstruction errors (right)

When there is no regularization, IRhybrid_lsqr is essentially an LSQR
algorithm, which exhibits semi-convergence behavior. When GCV is used instead
of the weighted-GCV, the geometry and image parameter errors also exhibit
semi-convergence behavior. The weighted-GCV helps stabilize LSQR’s convergence
and thus produces the best result.

###  4.5 Imfil Budget

We explore the effect of evaluation budgets in the nonlinear least squares
solver on the geometry and reconstruction errors. We set $N_{A}=10$, and use
$budget=10,100,1000,10000$. Since the budget size may greatly affect the
nonlinear least squares solutions, we explore its effects on BCDS without any
acceleration.

Figure 11: Comparison of budgets: geometry errors (left), reconstruction
errors (right)

As we can see in Figure 11, $100$, $1000$, and $10000$ are equally good. Both
geometry and reconstruction errors are very small when $100$, $1000$, $10000$
are used. Thus, $100$ is the best budget out of the four because any more
evaluation beyond $100$ does not make the solution better. When budget is
$1000$ and $10000$ respectively, we wasted many evaluations without making any
progress. When $10$ is used, even though we get convergent results earlier,
the small budget causes the algorithm to terminate before it finds a better
solution.

###  4.6 Nonlinear Least Squares Solver

In this section, we compare BCDS with different nonlinear least squares
solvers – fminbnd and imfil. Set $N_{A}=10$ and $budget=100$.

Figure 12: Comparison of nonlinear least squares solvers: geometry errors
(left), reconstruction errors (right)

In Figure 12, fminbnd exhibits semi-convergence behavior whereas imfil
converges. Under our problem settings, fminbnd reconstructs geometry
parameters around $10$s faster than imfil on average. However, fminbnd tends
to exhibit semi-convergence behavior whereas iterations in imfil are more
stable. Since in practical cases it is impossible to know the exact iteration,
at which the geometry errors are lowest in the semi-convergence case, imfil is
more reliable in solving the nonlinear least squares problem.

###  4.7 Geometry Parameters: $d$ and $\delta\theta$

We show that our alternating scheme as well as the accelerated version can
also be applied to the case when geometry parameters $r=(d,\delta\theta)$.The
nonlinear least squares problem becomes:

$\underset{x,\delta\theta,d}{\arg\min}||A(\delta\theta,d)x-b||.$

In all previous experiments we assume there was no error in angles. If the
angles contain errors as well, the image quality could also be affected, as
illustrated in Figure 13. The experiment setting is $N_{A}=10$ and
$budget=100$.

Figure 13: True image (top left), image obtained by considering both
$\delta\theta$ and $d$ (top right), image obtained by using correct $d$ but
incorrect $\delta\theta$ (bottom left), image obtained by using incorrect
geometry parameters $d$ and $\delta\theta$ (bottom right).

Figure 14: Comparison of BCDS and ABCDS when correcting for errors in source-
to-object distance and angles: geometry errors (left), reconstruction errors
(middle), errors in angles (right)

As shown in Figure 14, both ABCDS and BCDS converge for all three parameters
$d$, $\delta\theta$, and $x$. ABCDS has accelerated the convergence for all
three geometry parameters, especially for the reconstruction errors. The only
difference between this problem, $r=(d,\delta\theta)$, and the previous
problem, $r=d$, is we solve for two geometry parameters for each nonlinear
least squares problem instead of one.

##  5 Conclusion

In this paper, we propose an accelerated alternating minimization scheme to
solve X-ray tomography problems. The linear least squares problem is solved by
a weighted hybrid LSQR algorithm with Tikhonov regularization. The nonlinear
least squares problem is solved by implicit filtering. We also investigated
ABCDS and Anderson mixing to accelerate convergence. We have the following
findings:

  * 1.

BCDS runs faster than the normal BCD because the separability of parameters
allow us to solve each entry of $r$ independently. Also, since solving for
each parameter separately makes the nonlinear least squares problem much
easier than solving all parameters at once, BCDS is faster at finding a good
solution and converge faster than BCD.

  * 2.

ABCDS is much less computationally costly than BCDS with Anderson acceleration
because Anderson acceleration requires solving a linear least squares problem
in each iteration when choosing the weight. Although both methods converge,
ABCDS have better acceleration effects than Anderson acceleration.

  * 3.

The weighted hybrid LSQR algorithm with Tikhonov regularization stabilizes the
convergence more than applying an unweighted GCV to the LSQR algorithm. When
no regularization method is used, the LSQR algorithm exhibits semi-convergence
behavior. It reconstructs the solution at earlier iterations but noise at
later iterations.

  * 4.

Choosing an appropriate budget for implicit filtering is important. When the
budget is chosen too small, better solutions are not explored. When budget is
chosen too big, many function evaluations are wasted without making progress.
A suggested number in the imfil documentation is $10*N^{2}$, where $N$ is the
length of $r$. We found this formula does not always give the appropriate
budget. Since we solve for each geometry parameter using separability, the
dimension of each small problem is one. But, we have found $100$ to be the
best budget for $32\times 32$ test problem with $N_{A}=10$.

  * 5.

imfil is a better nonlinear least squares solver for our algorithm than
MATLAB’s fminbnd because the latter tends to exhibit semi-convergence
behaviors.

The ABCDS method has shown its success in our tomographic reconstruction
problems. We believe this algorithm can be used to effectively solve X-ray
tomography problems that have variations in the geometry parameter. A future
direction towards more improvement would be adapting, applying, and advancing
this algorithm on X-ray images produced in clinical trials.

## Acknowledgments

This work was done as part of an undergraduate honors thesis project at Emory
University, and was partially supported by the U.S. National Science
Foundation under Grant DMS-1819042. The author would like to thank Chang Meng
and James Nagy for their guidance and mentorship.

## References

  * [1] J. Chung, J. G. Nagy, and D. P. O’Leary, A weighted GCV method for Lanczos hybrid regularization, Electronic Transactions on Numerical Analysis, 28 (2008).
  * [2] Y. Cui, D. Sun, and K. Toh, Computing the best approximation over the intersection of a polyhedral set and the doubly nonnegative cone, SIAM Journal on Optimization, 29 (2019), pp. 2785–2813, https://doi.org/10.1137/18M1175562.
  * [3] D. Cummins, T. Filloon, and D. Nychka, Confidence intervals for nonparametric curve estimates: Toward more uniform pointwise coverage, Journal of the American Statistical Association, 96 (2001), pp. 233–246, https://doi.org/10.2307/2670362.
  * [4] J. H. Friedman and B. W. Silverman, Flexible parsimonious smoothing and additive modeling, Technometrics, 31 (1989).
  * [5] S. Gazzola, P. C. Hansen, and J. G. Nagy, IR tools: A MATLAB package of iterative regularization methods and large-scale test problems, 2018, https://arxiv.org/abs/1712.05602.
  * [6] P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion, Society for Industrial and Applied Mathematics, Philadelphia, 1998.
  * [7] P. C. Hansen and M. Saxild-Hansen, AIR tools — a MATLAB package of algebraic iterative reconstruction methods, Journal of Computational and Applied Mathematics, 236 (2012), pp. 2167–2178, https://doi.org/10.1016/j.cam.2011.09.039.
  * [8] C. T. Kelley, Implicit filtering, Society for Industrial and Applied Mathematics, (2011).
  * [9] D. P. O’Leary and B. W. Rust, Variable projection for nonlinear least squares problems, Computational Optimization and Applications, 54 (2013), pp. 579–593.
  * [10] V. Pathuri-Bhuvana, S. Schuster, and A. Och, Joint calibration and tomography based on separable least squares approach with constraints on linear and non-linear parameters, in 2020 28th European Signal Processing Conference (EUSIPCO), IEEE, 2021, pp. 1931–1935.
  * [11] Å. Björck, Numerical Methods for Least Squares Problems, Society for Industrial and Applied Mathematics, 1996, https://doi.org/10.1137/1.9781611971484.
  * [12] L. Shepp and B. F. Logan, The Fourier reconstruction of a head section, IEEE Transactions on Nuclear Science., (1974).
  * [13] R. Vio, P. Ma, W. Zhong, J. G. Nagy, L. Tenorio, and W. Wamsteker, Estimation of regularization parameters in multiple-image deblurring, A&A, 423 (2004), pp. 1179–1186, https://doi.org/10.1051/0004-6361:20047113.
  * [14] H. F. Walker, Anderson acceleration: Algorithms and implementations, WPI Math. Sciences Dept. Report MS-6-15-50, (2011).
  * [15] S. J. Wright, Coordinate descent algorithms, Mathematical Programming, 151 (2015), pp. 3–34, https://doi.org/10.1007/s10107-015-0892-3.
  * [16] J. Zhang, B. O’Donoghue, and S. Boyd, Globally convergent type-i anderson acceleration for nonsmooth fixed-point iterations, SIAM Journal on Optimization, 30 (2020), pp. 3170–3197, https://doi.org/10.1137/18M1232772.