Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

SIRT objective #1790

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

SIRT objective #1790

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
fa559d8
SIRT objective fix and unit tests
MargaretDuff Apr 23, 2024
a6f37d3
Tolerance of the test
MargaretDuff Apr 23, 2024
6b2e404
Tolerance of the test - attempt 2
MargaretDuff Apr 23, 2024
4e751ae
Merge branch 'master' into SIRT_objective
MargaretDuff May 3, 2024
6426c8b
Merge branch 'master' into SIRT_objective
MargaretDuff Jul 5, 2024
1216d54
SIRT objective
MargaretDuff Jul 5, 2024
e3daef1
Jakob's suggestion
MargaretDuff Jul 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 10 additions & 3 deletions Wrappers/Python/cil/optimisation/algorithms/SIRT.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,8 @@ class SIRT(Algorithm):
:math:`\mathrm{prox}_{C}` is the projection over a set :math:`C`,
and :math:`\omega` is the relaxation parameter.

Note that SIRT is equivalent to preconditioned gradient descent on a weighted least squares objective (weighted by :math:`M`) preconditioned by the sensitivity matrix (:math:`D`).
Thus the objective calculation for SIRT is the weighted least squares objective :math:`\frac{1}{2}\|A x - b\|^{2}_{M}` where :math:`\|y\|_M^2:=y^TMy`.
Parameters
----------

Expand Down Expand Up @@ -106,6 +108,7 @@ def set_up(self, initial, operator, data, lower=None, upper=None, constraint=Non
self.data = data

self.r = data.copy()


self.constraint = constraint
if constraint is None:
Expand Down Expand Up @@ -201,9 +204,13 @@ def update(self):
self.x=self.constraint.proximal(self.x, tau=1)

def update_objective(self):
r"""Returns the objective
r"""Returns the weighted least squares objective

.. math:: \frac{1}{2}\|A x - b\|^{2}
.. math:: \frac{1}{2}\|A x - b\|^{2}_{M}

"""
self.loss.append(0.5*self.r.squared_norm())
y = self.operator.direct(self.x)
y.subtract(self.data, out = y)
wy=self.M.multiply(y)
self.objective.append(0.5* y.dot(wy))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we want the 0.5? Could be omitted for (slight) speed-up and objective function rewritten to omit this. Should think about consistency with defaults in related operations including LeastSquares and CGLS though.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There isn't a 0.5 in CGLS (it uses \min || A x - b ||^2_2) and the default for least squares is 1.0. So in some ways, SIRT would be odd to include the 0.5. Mathematically, it is optimising the same thing. Shall we remove it?

Comment on lines +212 to +215
Copy link
Member

@gfardell gfardell Jul 10, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is adding an additional 2 copies of the data in to memory so we need to try to improve it.

wy is the same as the residuals r in update and this is still in memory. So instead I think we can square r (in place) and then weight it with M (again in place). As r is not reused by the algorithm it's ok using the memory for this. The first step of the algorithm recomputes the new residuals in to r.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tried this but got numerical errors (effectively from dividing and multiplying by M). I will try to debug this again.


24 changes: 20 additions & 4 deletions Wrappers/Python/test/test_algorithms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -962,7 +962,10 @@ def test_SIRT_relaxation_parameter(self):
alg.run(20,verbose=0)
np.testing.assert_array_almost_equal(alg.x.array, self.b2.array)

np.testing.assert_almost_equal(0.5 * alg.D.array, alg._Dscaled.array)

np.testing.assert_almost_equal(0.5 *alg.D.array, alg._Dscaled.array)



def test_SIRT_nan_inf_values(self):
Aop_nan_inf = self.Aop
Expand Down Expand Up @@ -999,6 +1002,16 @@ def test_SIRT_remove_nan_or_inf_with_BlockDataContainer(self):

self.assertFalse(np.any(sirt.D == np.inf))


def test_SIRT_objective(self):
alg = SIRT(initial=self.initial, operator=self.Aop,
data=self.bop)

alg.run(10)
f=LeastSquares(self.Aop, self.bop, c=0.5, weight=alg.M )
np.testing.assert_allclose(alg.get_last_loss(), f(alg.x) , rtol=1e-5)


def test_SIRT_with_TV(self):
data = dataexample.SIMPLE_PHANTOM_2D.get(size=(128, 128))
ig = data.geometry
Expand All @@ -1012,7 +1025,9 @@ def test_SIRT_with_TV(self):
fista = FISTA(initial=initial, f=f, g=constraint, max_iteration=1000)
fista.run(100, verbose=0)
self.assertNumpyArrayAlmostEqual(fista.x.as_array(), sirt.x.as_array())

f=LeastSquares(A, data ,c=0.5, weight=sirt.M )
np.testing.assert_allclose(sirt.get_last_loss(), f(sirt.x), rtol=1e-5)

def test_SIRT_with_TV_warm_start(self):
data = dataexample.SIMPLE_PHANTOM_2D.get(size=(128, 128))
ig = data.geometry
Expand All @@ -1023,8 +1038,9 @@ def test_SIRT_with_TV_warm_start(self):
max_iteration=150, constraint=constraint)
sirt.run(25, verbose=0)

self.assertNumpyArrayAlmostEqual(
sirt.x.as_array(), ig.allocate(0.25).as_array(), 3)
self.assertNumpyArrayAlmostEqual(sirt.x.as_array(), ig.allocate(0.25).as_array(),3)
f=LeastSquares(A, data ,c=0.5, weight=sirt.M )
np.testing.assert_allclose(sirt.get_last_loss(), f(sirt.x), rtol=1e-5)


class TestSPDHG(unittest.TestCase):
Expand Down
54 changes: 29 additions & 25 deletions Wrappers/Python/test/test_preconditioners.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@

from cil.optimisation.algorithms import SIRT, GD, ISTA, FISTA
from cil.optimisation.functions import LeastSquares, IndicatorBox
from cil.framework import ImageGeometry, VectorGeometry
from cil.framework import ImageGeometry, VectorGeometry, VectorData
from cil.optimisation.operators import IdentityOperator, MatrixOperator

from cil.optimisation.utilities import Sensitivity, AdaptiveSensitivity, Preconditioner
Expand Down Expand Up @@ -72,9 +72,17 @@ def test_sensitivity_calculation(self):

def test_sensitivity_gd_against_sirt(self):

ig = ImageGeometry(12, 13, 14)
data = ig.allocate('random', seed=4)
A = IdentityOperator(ig)
n = 50
m = 500

A = np.random.uniform(0, 1, (m, n)).astype('float32')
b = (A.dot(np.random.randn(n)) + 0.1 *
np.random.randn(m)).astype('float32')

A = MatrixOperator(A)
data = VectorData(b)
ig=A.domain


sirt = SIRT(ig.allocate(0), A, data, update_objective_interval=1)
sirt.run(10)
Expand All @@ -91,26 +99,28 @@ def test_sensitivity_gd_against_sirt(self):
precond_pwls = GD(initial=ig.allocate(0), objective_function=f, preconditioner=preconditioner,
max_iteration=100, update_objective_interval=1, step_size=step_size)

def correct_update_objective(alg):
# SIRT computes |Ax_{k} - b|_2^2
# GD with weighted LeastSquares computes the objective included the weight, so we remove the weight
return 0.5*(alg.objective_function.A.direct(alg.x) - alg.objective_function.b).squared_norm()

precond_pwls.run(10)
np.testing.assert_allclose(
sirt.solution.array, precond_pwls.solution.array, atol=1e-4)
np.testing.assert_allclose(sirt.get_last_loss(
), correct_update_objective(precond_pwls), atol=1e-4)
np.testing.assert_allclose(f(sirt.x), precond_pwls.get_last_loss(), atol=1e-4)


def test_sensitivity_ista_against_sirt(self):

ig = ImageGeometry(12, 13, 14)
data = ig.allocate('random')
A = IdentityOperator(ig)
n = 50
m = 500

A = np.random.uniform(0, 1, (m, n)).astype('float32')
b = (A.dot(np.random.randn(n)) + 0.1 *
np.random.randn(m)).astype('float32')

A = MatrixOperator(A)
data = VectorData(b)
ig=A.domain

sirt = SIRT(ig.allocate(0), A, data, lower=0,
update_objective_interval=1)
sirt.run(10)
sirt.run(15)

M = 1./A.direct(A.domain_geometry().allocate(value=1.0))
f = LeastSquares(A=A, b=data, c=0.5, weight=M)
Expand All @@ -122,19 +132,13 @@ def test_sensitivity_ista_against_sirt(self):
max_iteration=100, update_objective_interval=1, step_size=step_size)
self.assertEqual(alg.preconditioner, None)

precond_pwls = GD(initial=ig.allocate(0), objective_function=f, preconditioner=preconditioner,
precond_pwls = ISTA(initial=ig.allocate(0), f=f, g=g, preconditioner=preconditioner,
max_iteration=100, update_objective_interval=1, step_size=step_size)

def correct_update_objective(alg):
# SIRT computes |Ax_{k} - b|_2^2
# GD with weighted LeastSquares computes the objective included the weight, so we remove the weight
return 0.5*(alg.objective_function.A.direct(alg.x) - alg.objective_function.b).squared_norm()

precond_pwls.run(10)

precond_pwls.run(15)
np.testing.assert_allclose(
sirt.solution.array, precond_pwls.solution.array, atol=1e-4)
np.testing.assert_allclose(sirt.get_last_loss(
), correct_update_objective(precond_pwls), atol=1e-4)
np.testing.assert_allclose(f(sirt.x), precond_pwls.get_last_loss(), atol=1e-4)



Expand Down
Loading