BSGD-TV: A parallel algorithm solving total variation constrained image reconstruction problems

1 Introduction

Our algorithm is inspired by applications in computed tomography (CT), where the efficient inversion of large sparse linear systems is required [1]: $y \approx A x_{t r u e}$ , where $x_{t r u e} \in R^{c}$ is the vectorised version of a 3D image that is to be reconstructed and $A \in R^{r \times c}$ is an X-ray projection model. $y \in R^{r}$ are the vectorised noisy projections. We are interested in minimizing $f (x) + g (x)$ , where $f (x)$ is quadratic and $g (x)$ is convex but non-smooth [2]. For example:

x^{⋆} = arg min x (y - A x)^{T} (y - A x)      f (x) + 2 λ TV (x)      g (x),

(1)

where $λ$ is a relaxation parameter and $TV (x)$ is the total variation (TV) of the image $x$ . For 2D images, it is defined as:

TV (x) = \sum s, t \sqrt{(x_{s, t} - x_{s - 1, t})^{2} + (x_{s, t} - s_{s, t - 1})^{2}},

(2)

where $x_{s, t}$ is the intensity of image pixel in row $s$ and column $t$ .

We recently introduced a parallel reconstruction algorithm called coordinate-reduced stochastic gradient descent (CSGD) to minimize quadratic objective function

f (x)

[3]. We here introduce a slight modification by simplifying the step length calculation and show that the modified version converges to the least squares solution of

f (x)

. We will call this modified algorithm block stochastic gradient descend (BSGD). We combine BSGD with an iterative shrinkage/thresholding (ISTA-type) step [4] to solve Eq.1. The new algorithm, called BSGD-TV, is compared with block ADMM-TV [5]

, an algorithm sharing the same parallel architecture and the same communication cost. Simulation results show that BSGD-TV is significantly faster as it requires significantly fewer matrix vector products compared to block ADMM-TV.

2 BSGD-TV Algorithm

2.1 Algorithm description

BSGD works on blocks of $x$ and $y$ . We assume that $A$ is divided into $M$ row blocks and $N$ column blocks. Let ${x_{J_{j}}}_{j = 1}^{N}$ and ${y_{I_{i}}}_{i = 1}^{M}$ be sub-vectors of $x$ and $y$ and let $A_{I_{i}}^{J_{j}}$ be the associated block of matrix $A$ so that $y_{I_{i}} \approx \sum_{j = 1}^{N} A_{I_{i}}^{J_{j}} x_{J_{j}}$ . Our algorithm splits the optimization into blocks, so that each parallel process only computes using a single block $x_{J_{j}}$ and $y_{I_{i}}$ for some $J_{j} \in {J_{j}}_{j = 1}^{N}$ and $I_{i} \in {I_{i}}_{i = 1}^{M}$

. Each process also requires an estimate of the current residual

r_{I_{i}}

and computes a vector

z_{I_{i}}^{j}

, both of which are of the same size as

y_{I_{i}}

. The main steps (ignoring initialisation) are described in Algo.1.

1: for

e p o c h = 1, 2, \dots

2: for

M \times N

pairs

{i, j}

drawn randomly without replacement in parallel do

{^g}_{J_{j}}^{i} = 2 (A_{I_{i}}^{J_{j}})^{T} r_{I_{i}}

z_{I_{i}}^{j} = A_{I_{i}}^{J_{j}} x_{J_{j}}

5: end for

r = y - \sum_{j} z^{j}

g = \sum_{i} {^g}^{i}

x = x + μ g

(constant

μ > 0

)

x = arg {min}_{t} ∥ t - x ∥^{2} + 2 μ λ % T V (t)

10: end for

Algorithm 1 BSGD-TV algorithm

To effectively solve line 9, we here adopt method proposed in [6].

2.2 BSGD Convergence

BSGD without the proximal operator ( $λ$ =0 in Eq.1), and with parallelization over all subsets can be shown to converge to the least squares solution. To see this, we write the update of $x$ as

$x^{k + 1}$	$= x^{k} + μ g^{k}$	(3)
	$= x^{k} + 2 μ A^{T} (y - N \sum j = 1 (z^{j})^{k - 1})$
	$= x^{k} + 2 μ A^{T} (y - A x^{k - 1})$

In this form, BSGD is similar to gradient descent but uses an old gradient. Assume that there is a fixed point $x^{⋆}$ defined by $x^{⋆} = x^{⋆} + 2 μ A^{T} (y - A x^{⋆})$ . Note that the fixed point condition implies that, if $A$ is full column rank, then $x^{⋆} = (A^{T} A)^{- 1} A^{T} y$ . Thus the fixed point is the least squares solution. Theorem 2.1 states the conditions on parameter $μ$ for convergence when all subsets ${I_{i}}_{i = 1}^{M}$ and ${J_{j}}_{j = 1}^{N}$

are selected within one epoch.

Theorem 2.1.

If $μ \in (0, \frac{1}{2 u_{m a x}})$ , where $u_{m a x}$

is the maximum eigenvalue of

A^{T} A

and assume

A

is full column rank, then BSGD without the TV operator (

λ

=0 in Eq.1), and with parallelization over all subsets converges to the least squares solution

x^{⋆}

Proof of Theorem 2.1.

The iteration in Eq.3 can be written as

	$[\begin{matrix} x^{k} x^{k + 1} \end{matrix}] =$	$[\begin{matrix} 0 I - 2 μ A^{T} A I \end{matrix}] [\begin{matrix} x^{k - 1} x^{k} \end{matrix}] + [\begin{matrix} 0 2 μ A^{T} y \end{matrix}]$		(4)
		$= M [\begin{matrix} x^{k - 1} x^{k} \end{matrix}] + [\begin{matrix} 0 2 μ A^{T} y \end{matrix}] .$		(4)

Standard convergence results for iterative method of this type with fixed $M$ require the spectral radius of $M$ to be less than 1 [7]. Let $v$ be any (possibly complex valued) eigenvalue of $M$ , i.e. $v$ satisfies det $(M - v I) = 0$ . It is straightforward to obtain:

det ([\begin{matrix} - v I I - 2 μ A^{T} A I - v I \end{matrix}]) = det% (A^{T} A - \frac{v - v^{2}}{2 μ} I) = 0

(5)

By Eq.5, we see that eigenvalues $u$ of $A^{T} A$ correspond to

u = \frac{v - v^{2}}{2 μ},

(6)

Eigenvalues of $M$ are then given by

v_{1} = \frac{1 + \sqrt{1 - 8 μ u}}{2}, v_{2} = \frac{1 - \sqrt{1 - 8 μ u}}{} 2 .

(7)

As the spectral radius of $M$ corresponds to the largest magnitude of the eigenvalues of $M$ , we require $| v_{1} | < 1$ and $| v_{2} | < 1$ to ensure the convergence of the algorithm. $A^{T} A$ is a positive definite matrix and thus has only positive, real valued eigenvalues $u$ . Thus $v_{1}$ and $v_{2}$ are real valued if $0 < μ \leq \frac{1}{8 u}$ and complex valued if $μ$ is $\frac{1}{8 u} < μ$ . In the complex case, it is easy to see that $| v_{1} | < 1$ and $| v_{2} | < 1$ if $μ < \frac{1}{2 u}$ , implying that the acceptable range of $μ$ is $(0, \frac{1}{2 u_{m a x}})$ . ∎

Theorem 2.2 gives a general convergence condition when applying BSGD-TV to solve Eq.1.

Theorem 2.2.

If the constant step length $μ$ satisfies

f (x^{k + 1}) < f (x^{k}) + (x^{k + 1} - x^{k})^{T} \nabla f (x^{k - 1}) + \frac{1}{2 μ} ∥ x^{k + 1} - x^{k} ∥^{2},

(8)

where $f (x)$ is defined in Eq.1, then BSGD-TV converges to the optimal solution of Eq.1.

Proof of Theorem 2.2.

With parallelization over all subsets, BSGD-TV computes

		${^x}^{k + 1} = x^{k} - μ \nabla f (x^{k - 1})$		(9)
		$xk+1=argminx{2μg(x)+∥x−^xk+1∥}.$		(9)

We define a function $Q$ as

	$Q (x, x^{k}, x^{k - 1}) =$	$f (x^{k}) + (x - x^{k})^{T} \nabla f (x^{k - 1})$		(10)
		$+ \frac{1}{2 μ} ∥ x - x^{k} ∥^{2} + g (x),$		(10)

where $∥ \cdot ∥^{2}$ is the squared $ℓ_{2}$ norm. The fact that $arg {min}_{x} {Q (x, x^{k}, x^{k - 1})} \equiv x^{k + 1}$ means that:

Q (x^{k + 1}, x^{k}, x^{k - 1}) < Q (x^{k}, x^{k}, x^{k - 1}) \equiv f (x^{k}) + g (x^{k})

(11)

Finally, the definition of $Q$ and the requirement on the step length $μ$ in Eq.8, mean that $f (x^{k + 1}) + g (x^{k + 1}) < Q (x^{k + 1}, x^{k}, x^{k - 1}) < f (x^{k}) + g (x^{k})$ holds, so that BSGD-TV converges to the fixed point of Eq.1. ∎

3 Simulations

We show experimentally that the method also converges when only a fraction $α$ and $γ$ of subsets of ${x_{J_{j}}}_{j = 1}^{N}$ and ${y_{I_{i}}}_{i = 1}^{M}$ are randomly selected to calculate the $g$ at each iteration. The simulation geometry is shown in Fig.1.

Figure 1: Using a fan-beam x-rays geometry to scan a 2D Shepp-Logan phantom. Projections are taken at $10^{\circ}$ intervals. The image $x$ is partitioned into 4 subsets ${J_{j}}_{j = 1}^{4}$ and the total projections are also partitioned into 4 subsets.

We add Gaussian noise to the projections so that the SNR of $y$ is 17.7 dB. We define the relative error as $\frac{∥ x_{d i f} ∥}{∥ x_{t r u e} ∥}$ , where $∥ x_{d i f} ∥$ is the $ℓ_{2}$ norm of the difference between reconstructed image vector and the original vector $x_{t r u e}$ . Convergence is shown in Fig.2a. We plot relative error against epochs, where an epoch is a normalised iteration count that corrects for the fact that the stochastic version of our algorithm only updates a subset of elements at each iteration.

BSGD-TV and ADMM-TV are faster than ISTA in terms of epochs. However, ADMM-TV is significantly slower than BSGD-TV, because ADMM-TV requires matrix inversions at each iteration, while BSGD does not. Even when implementing ADMM-TV using as few conjugate gradient iterations per step as possible, as shown in Fig.2b, BSGD-TV is more computationally efficient in terms of the number of required matrix vector multiplications [3]. Compared to ISTA and GD, our block method allows these computations to be fully parallelised which would enable to reconstruct large scale CT reconstructions while the computation node have limited storage capacity.

Figure 2: The step length $μ$ for ISTA, gradient descent (GD, which solves $f (x)$ in Eq.1 ) and BSGD is $6 e - 4$ and $λ$ in Eq.1 is $0.1$ . For BSGD and ADMM, $A$ is divided into $4 \times 4$ sub-matrices while ISTA processes $A$ as a whole. (a): Relative error vs. epoch. The high relative error of GD suggests the necessity of incorporating the TV norm. (b): BSGD-TV uses significantly fewer matrix-vector multiplications compared to ADMM-TV.

4 Conclusion

BSGD-TV is a parallel algorithm for large scale TV constrained CT reconstruction. It is similar to the popular ISTA algorithm but is specially designed for optimisation in distributed networks. The advantage is that individual compute nodes only operate on subsets of $y$ and $x$ , which means they can operate with less internal memory. The method converges significantly faster than block-ADMM methods.

References

[1] X. Guo, “Convergence studies on block iterative algorithms for image reconstruction”, Applied Mathematics and Computation, 273: 525–534, 2016.
[2] E. Sidky and X.Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization”, Physics in Medicine $&$ Biology, 53(17):4777–4807, 2008.
[3] Y. Gao and T.Blumensath, “A Joint Row and Column Action Method for Cone-Beam Computed Tomography”, IEEE Transactions on Computational Imaging, 2018.
[4] A. Beck and M.Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems”, SIAM journal on imaging sciences, 2(1):183–202, 2009.
[5] N. Parikh and S. Boyd, “Block splitting for distributed optimization”, Mathematical Programming Computation, 6(1):77–102, 2014.
[6] A. Beck and M.Teboulle, “Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems”, IEEE Transactions on Image Processing, 18(11):2419–2434, 2009.
[7] Y. Saad, “Iterative methods for sparse linear systems”, siam,82, 2003.

BSGD-TV: A parallel algorithm solving total variation constrained image reconstruction problems

Authors

Enhanced total variation minimization for stable image reconstruction

Higher-order total variation approaches and generalisations

Block stochastic gradient descent for large-scale tomographic reconstruction in a parallel network

Compressed MRI Reconstruction Exploiting a Rotation-Invariant Total Variation Discretization

Hessian Schatten-Norm Regularization for Linear Inverse Problems

Image reconstruction from few views by L0-norm optimization

On and beyond Total Variation regularisation in imaging: the role of space variance

1 Introduction

2 BSGD-TV Algorithm

2.1 Algorithm description

2.2 BSGD Convergence

Theorem 2.1.

Proof of Theorem 2.1.

Theorem 2.2.

Proof of Theorem 2.2.

3 Simulations

4 Conclusion

References

BSGD-TV: A parallel algorithm solving total variation constrained image reconstruction problems

Authors

Related Research

Enhanced total variation minimization for stable image reconstruction

Higher-order total variation approaches and generalisations

Block stochastic gradient descent for large-scale tomographic reconstruction in a parallel network

Compressed MRI Reconstruction Exploiting a Rotation-Invariant Total Variation Discretization

Hessian Schatten-Norm Regularization for Linear Inverse Problems

Image reconstruction from few views by L0-norm optimization

On and beyond Total Variation regularisation in imaging: the role of space variance

This week in AI

1 Introduction

2 BSGD-TV Algorithm

2.1 Algorithm description

2.2 BSGD Convergence

Theorem 2.1.

Proof of Theorem 2.1.

Theorem 2.2.

Proof of Theorem 2.2.

3 Simulations

4 Conclusion

References