Added LSQR algorithm

[MaikeMM]

Description

Implementation of the LSQR algorithm (Paige and Saunders, 1982, ACM TOMS). LSQR is an algorithm that solves a least squares problem $\min_x|Ax-b|$. It is mathematically equivalent to CGLS, which is already implemented in CIL, but more stable in finite-precision arithmetic. This implementation of LSQR works in exactly the same way as how CGLS is implemented in CIL, i.e. wherever CGLS is called LSQR can be called without changing any syntax.

NEW: In a recent edit, the optional regularisation parameter alpha can now be passed to LSQR to become the Tikhonov regularised version.

Example Usage

Below an example heavily inspired by the CIL Demo on Tikhonov regularisation.

# Load components needed
from cil.framework import ImageGeometry, AcquisitionGeometry, ImageData
from cil.optimisation.algorithms import LSQR
from cil.plugins.astra import ProjectionOperator
from cil.utilities.display import show2D

# Load third-party inputs (using phantominator because the TomoPhantom plug-in does not work on mac for me)
import numpy as np    
import matplotlib.pyplot as plt
from phantominator import shepp_logan

# Set up toy problem
# Geometry
n = 256
ig = ImageGeometry(voxel_num_x=n, 
                   voxel_num_y=n, 
                   voxel_size_x=2/n, 
                   voxel_size_y=2/n)

num_angles = 180
ag = AcquisitionGeometry.create_Parallel2D()  \
                   .set_angles(np.linspace(0, 180, num_angles, endpoint=False))  \
                   .set_panel(n, 2/n)

# Shepp-Logan phantom
phantom2DArray = shepp_logan(n)
phantom2D = ImageData(array = phantom2DArray, deep_copy=False, geometry=ig)

# Load operator
device = "cpu"
A = ProjectionOperator(ig, ag, device)

# Create noise-less sinogram
sinogram = A.direct(phantom2D)

# Add noise to sinogram
background_counts = 10000 
counts = background_counts * np.exp(-sinogram.as_array())
noisy_counts = np.random.poisson(counts)
sino_out = -np.log(noisy_counts/background_counts)
sinogram_noisy = ag.allocate()
sinogram_noisy.fill(sino_out)

# Solve problem with LSQR
maxit = 25
initial = ig.allocate(0)
lsqr_simple = LSQR(initial=initial, operator=A, data=sinogram_noisy)
lsqr_simple.run(maxit, verbose=True)

# Plot the solution
plots = [lsqr_simple.solution, lsqr_simple.solution - phantom2D]
titles = ["LSQR solution","Difference from ground truth" ]
show2D(plots, titles, fix_range=[(-0.2,1.2),(-0.2,0.2)])

# Include Tikhonov regularisation without building a block operator
alpha = 0.1
lsqr_reg = LSQR(initial=initial, operator=A, data=sinogram_noisy, alpha = alpha)
lsqr_reg.run(maxit, verbose=True)

# Plot the solution
plots = [lsqr_reg.solution, lsqr_reg.solution - phantom2D]
titles = ["Regularised LSQR solution","Difference from ground truth" ]
show2D(plots, titles, fix_range=[(-0.2,1.2),(-0.2,0.2)])

:heart: Thanks for your contribution!

• [ ] I have performed a self-review of my code
• [ ] I have added docstrings in line with the guidance in the developer guide
• [ ] I have updated the relevant documentation
• [ ] I have implemented unit tests that cover any new or modified functionality
• [ ] CHANGELOG.md has been updated with any functionality change
• [ ] Request review from all relevant developers
• [ ] Change pull request label to 'Waiting for review'

Contribution Notes

Please read and adhere to the developer guide and local patterns and conventions.

• [ ] The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in CIL (the Work) under the terms and conditions of the Apache-2.0 License
• [ ] I confirm that the contribution does not violate any intellectual property rights of third parties

[jakobsj]

THanks so much @MaikeMM for adding this LSQR implementation! As discussed, it would be great and a big help for reviewing this if you could add a brief code snippet demonstrating how a user might run the algorithm. This could be some lines of code here or as a python/notebook attachment here. Also, we'd be very interested in including your notebook in the showcase :) https://github.com/TomographicImaging/CIL-User-Showcase it would make an excellent addition demonstrating how to implement a new CIL algorithm, assess it against an existing algorithm (CGLS) as well as compare the block-free version of Tikhonov

[MaikeMM]

Hi @jakobsj, so sorry for the delay but I've updated this comment now. I hope this kind of thing is what is required but please do let me know if you'd like anything more!

[MargaretDuff]

Closes #1992

[jakobsj]

Thank you very much @MaikeMM for adding an example. We will take a look and get back to you.

Any chance of contributing the nice notebook you demonstrated to the showcase? :)

[MaikeMM]

Hi Margaret, please see the updates to include Tikhonov regularisation!

[MargaretDuff]

Thanks Maike! Will take a look

[MargaretDuff]

I had a play with this to try and optimise memory usage. Compared to CGLS, I still need to get rid of one stored image copy. Will discuss with Gemma and Mariam next week