From 3e4f80e16e32c29e89d90fbf7bea64e40a78c9e4 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 18 Apr 2023 21:32:30 +0100 Subject: [PATCH 01/18] add symmetric flag to PowerMethod closes #1468 Signed-off-by: Edoardo Pasca --- .../cil/optimisation/operators/Operator.py | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index 764578d44c..98a65f8c21 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -141,13 +141,16 @@ def adjoint(self,x, out=None): raise NotImplementedError @staticmethod - def PowerMethod(operator, max_iteration=10, initial=None, tolerance = 1e-5, return_all=False): + def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False, symmetric=False): r"""Power method or Power iteration algorithm The Power method computes the largest (dominant) eigenvalue of a square matrix in magnitude, e.g., absolute value in the real case and module in the complex case. - For the non-square case, the power method is applied for the matrix :math: A^{T}*A. + Since we cannot estimate if the operator is diagonisable, the power method is applied on the diagonizable operator + :math:`A^{T}*A`. + + If the operator :math:`A` is diagonizable, the user can pass the flag `symmetric=True`. Parameters ---------- @@ -161,6 +164,8 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance = 1e-5, ret Stopping criterion for the Power method. Check if two consecutive eigenvalue evaluations are below the tolerance. return_all: boolean, default = False Toggles the verbosity of the return + symmetric: boolean, default False + Set this to True if the operator is symmetric Returns @@ -188,16 +193,6 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance = 1e-5, ret 2.0005647295658866 """ - - # Default case: non-symmetric - symmetric = False - try: - if operator.domain_geometry()==operator.range_geometry(): - symmetric = True - except AssertionError: - # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972 - pass - if initial is None: x0 = operator.domain_geometry().allocate('random') else: From 8ec43ac64147f95b5b2b9781636c9cbe77879ab2 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Wed, 23 Aug 2023 15:41:01 +0000 Subject: [PATCH 02/18] Changed symmetric to square and changed the documentation --- .../cil/optimisation/operators/Operator.py | 40 ++++++++++++------- Wrappers/Python/test/test_Operator.py | 1 + 2 files changed, 27 insertions(+), 14 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index 98a65f8c21..266bcdc381 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -141,17 +141,16 @@ def adjoint(self,x, out=None): raise NotImplementedError @staticmethod - def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False, symmetric=False): + def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False): r"""Power method or Power iteration algorithm The Power method computes the largest (dominant) eigenvalue of a square matrix in magnitude, e.g., - absolute value in the real case and module in the complex case. - Since we cannot estimate if the operator is diagonisable, the power method is applied on the diagonizable operator - :math:`A^{T}*A`. - - If the operator :math:`A` is diagonizable, the user can pass the flag `symmetric=True`. + absolute value in the real case and modulus in the complex case. + If the matrix is not square. The algorithm computes the largest (dominant) eigenvalue of :math: A^{T}*A, returning the square root of this value. + + Parameters ---------- @@ -164,8 +163,8 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur Stopping criterion for the Power method. Check if two consecutive eigenvalue evaluations are below the tolerance. return_all: boolean, default = False Toggles the verbosity of the return - symmetric: boolean, default False - Set this to True if the operator is symmetric + square: boolean, default False + Set this to True if the operator is square Returns @@ -193,6 +192,13 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur 2.0005647295658866 """ + square = False + try: + if operator.domain_geometry()==operator.range_geometry(): + square = True + except AssertionError: + # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972 + pass if initial is None: x0 = operator.domain_geometry().allocate('random') else: @@ -212,10 +218,15 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur i = 0 while (i < max_iteration and diff > tolerance): i+=1 - + if i==max_iteration: + warnings.warn('In {} iterations, the difference in eigenvalue between the last two iterations was {} over the tolerance of {}. \n \ + First try a larger iteration number and then a different initial vector. If that does not work consider the operator. For convergence, \n \ + we need the operator to have an eigenvalue that is strictly greater in magnitude than the magnitude of its other eigenvalues'.format(i, diff, tolerance)) + operator.direct(x0, out = y_tmp) - if symmetric: + + if square: #swap datacontainer references tmp = x0 x0 = y_tmp @@ -223,17 +234,18 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur else: operator.adjoint(y_tmp,out=x0) + # Get eigenvalue using Rayleigh quotient: denominator=1, due to normalization - x0_norm = x0.norm() + x0_norm=x0.norm() x0 /= x0_norm eig_new = numpy.abs(x0_norm) - if not symmetric: + if not square: eig_new = numpy.sqrt(eig_new) - diff = numpy.abs(eig_new - eig_old) eig_list.append(eig_new) - eig_old = eig_new + eig_old = eig_new + if return_all: return eig_new, i, x0, eig_list diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index 01988899c5..db16e76efa 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -33,6 +33,7 @@ from cil.utilities import dataexample import logging from testclass import CCPiTestClass +import numpy as np initialise_tests() From 03d2805b57456b608b4f67710d482be9b78c841b Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Thu, 24 Aug 2023 08:12:26 +0000 Subject: [PATCH 03/18] Unit tests and case of nil-potent matrix --- .../cil/optimisation/operators/Operator.py | 6 ++++-- Wrappers/Python/test/test_Operator.py | 19 +++++++++++++++++-- 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index 266bcdc381..ba6df45668 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -219,9 +219,9 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur while (i < max_iteration and diff > tolerance): i+=1 if i==max_iteration: - warnings.warn('In {} iterations, the difference in eigenvalue between the last two iterations was {} over the tolerance of {}. \n \ + warnings.warn('In {} iterations, the difference in eigenvalue between the last two iterations was {} larger than the tolerance of {}. \n \ First try a larger iteration number and then a different initial vector. If that does not work consider the operator. For convergence, \n \ - we need the operator to have an eigenvalue that is strictly greater in magnitude than the magnitude of its other eigenvalues'.format(i, diff, tolerance)) + we need the operator to have an eigenvalue that is strictly greater in magnitude than the magnitude of its other eigenvalues'.format(i, diff, tolerance)) operator.direct(x0, out = y_tmp) @@ -237,6 +237,8 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur # Get eigenvalue using Rayleigh quotient: denominator=1, due to normalization x0_norm=x0.norm() + if x0_norm Date: Thu, 24 Aug 2023 08:17:39 +0000 Subject: [PATCH 04/18] Spelling error --- Wrappers/Python/test/test_Operator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index 5deb45ccf7..738bf1f1ad 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -312,7 +312,7 @@ def test_PowerMethod(self): M1op = MatrixOperator(M1) res1 = M1op.PowerMethod(M1op,5) - # 2x2 nil + # 2x2 matrix, max absolute eigenvalue is not unique and initial vector chosen for non-convergence with self.assertWarns(Warning): M1=numpy.array([[2.,1.], [0.,-2.]]) M1op = MatrixOperator(M1) From 05aaa53374b5087d465c2dfc9f9acc5edf9b3fe4 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Fri, 25 Aug 2023 14:12:24 +0000 Subject: [PATCH 05/18] Playing with convergence --- Wrappers/Python/test/test_Operator.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index 738bf1f1ad..4d2cf5c1e4 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -318,6 +318,12 @@ def test_PowerMethod(self): M1op = MatrixOperator(M1) res1 = M1op.PowerMethod(M1op,100, initial=DataContainer(numpy.array([1.,1.]))) + # 2x2 matrix, max absolute eigenvalue is not unique and initial vector chosen for convergence + + M1=numpy.array([[2.,1.,0.],[0.,1.,1.], [0.,0.,1.]]) + M1op = MatrixOperator(M1) + res1 = M1op.PowerMethod(M1op,100) + numpy.testing.assert_almost_equal(res1,2., decimal=4) # Gradient Operator (float) From 0723002da86010731d21e6c660cdf45d49a58ac9 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Mon, 18 Sep 2023 14:07:15 +0000 Subject: [PATCH 06/18] Added option range_is_domain in the init --- .../cil/optimisation/operators/Operator.py | 47 +++++++++++-------- 1 file changed, 28 insertions(+), 19 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index ba6df45668..8da95c6e73 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -141,14 +141,14 @@ def adjoint(self,x, out=None): raise NotImplementedError @staticmethod - def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False): + def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False,range_is_domain=None ): r"""Power method or Power iteration algorithm The Power method computes the largest (dominant) eigenvalue of a square matrix in magnitude, e.g., absolute value in the real case and modulus in the complex case. - If the matrix is not square. The algorithm computes the largest (dominant) eigenvalue of :math: A^{T}*A, returning the square root of this value. + If the matrix is not square. The algorithm computes the largest (dominant) eigenvalue of :math: A^{T}*A :math:, returning the square root of this value. Parameters @@ -163,8 +163,9 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur Stopping criterion for the Power method. Check if two consecutive eigenvalue evaluations are below the tolerance. return_all: boolean, default = False Toggles the verbosity of the return - square: boolean, default False - Set this to True if the operator is square + range_is_domain: `boolean`, default None + Set this to `True` to apply the power method directly on the operator, :math: A :math:, and `False` to apply the power method to :math:A^TA:math: before taking the square root of the result. + Leave as default `None` to determine this from domain and range geometry. Returns @@ -176,6 +177,8 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur Corresponding eigenvector of the dominant eigenvalue. Only returned if return_all is True. list of eigenvalues: :obj:`list` List of eigenvalues. Only returned if return_all is True. + Convergence: `boolean' + Check on wether the difference between the last two iterations is less than tolerance. Only returned if return_all is True. Examples -------- @@ -192,13 +195,18 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur 2.0005647295658866 """ - square = False - try: - if operator.domain_geometry()==operator.range_geometry(): - square = True - except AssertionError: - # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972 - pass + convergence_check=True + if range_is_domain is None: + square = False + try: + if operator.domain_geometry()==operator.range_geometry(): + square = True + except AssertionError: + # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972 + pass + else: + square=range_is_domain + if initial is None: x0 = operator.domain_geometry().allocate('random') else: @@ -218,11 +226,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur i = 0 while (i < max_iteration and diff > tolerance): i+=1 - if i==max_iteration: - warnings.warn('In {} iterations, the difference in eigenvalue between the last two iterations was {} larger than the tolerance of {}. \n \ - First try a larger iteration number and then a different initial vector. If that does not work consider the operator. For convergence, \n \ - we need the operator to have an eigenvalue that is strictly greater in magnitude than the magnitude of its other eigenvalues'.format(i, diff, tolerance)) - + operator.direct(x0, out = y_tmp) @@ -238,7 +242,9 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur # Get eigenvalue using Rayleigh quotient: denominator=1, due to normalization x0_norm=x0.norm() if x0_norm Date: Mon, 18 Sep 2023 14:56:15 +0000 Subject: [PATCH 07/18] Some testing --- Wrappers/Python/test/test_BlockOperator.py | 6 ++--- Wrappers/Python/test/test_Operator.py | 31 +++++++++++++--------- 2 files changed, 22 insertions(+), 15 deletions(-) diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py index 3e81ab4cad..bff62cd5af 100644 --- a/Wrappers/Python/test/test_BlockOperator.py +++ b/Wrappers/Python/test/test_BlockOperator.py @@ -26,7 +26,7 @@ from cil.framework import ImageGeometry, ImageData import numpy from cil.optimisation.operators import FiniteDifferenceOperator - +from cil.optimisation.operators import LinearOperator initialise_tests() class TestBlockOperator(unittest.TestCase): @@ -187,9 +187,9 @@ def test_FiniteDiffOperator(self): G = FiniteDifferenceOperator(ig, direction=0, bnd_cond = 'Neumann') logging.info("{} {}".format(type(u), str(u.as_array()))) logging.info(str(G.direct(u).as_array())) - + LinearOperator.PowerMethod(G, range_is_domain=False) # Gradient Operator norm, for one direction should be close to 2 - numpy.testing.assert_allclose(G.norm(), numpy.sqrt(4), atol=0.1) + numpy.testing.assert_allclose(LinearOperator.PowerMethod(G, range_is_domain=False), numpy.sqrt(4), atol=0.1) M1, N1, K1 = 200, 300, 2 ig1 = ImageGeometry(voxel_num_x = M1, voxel_num_y = N1, channels = K1) diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index 4d2cf5c1e4..ea756157d6 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -306,25 +306,32 @@ def test_PowerMethod(self): numpy.testing.assert_almost_equal(res1,3.1624439599276974, decimal=4) # 2x2 non-diagonalisable nilpotent matrix - with self.assertRaises(ValueError): - - M1=numpy.array([[0.,1.], [0.,0.]]) - M1op = MatrixOperator(M1) - res1 = M1op.PowerMethod(M1op,5) - + M1=numpy.array([[0.,1.], [0.,0.]]) + M1op = MatrixOperator(M1) + res1 = M1op.PowerMethod(M1op,5) + numpy.testing.assert_almost_equal(res1,0, decimal=4) + + # 2x2 non-diagonalisable nilpotent matrix where range_is_domain is False + M1=numpy.array([[0.,1.], [0.,0.]]) + M1op = MatrixOperator(M1) + res1 = M1op.PowerMethod(M1op,5, range_is_domain=False) + numpy.testing.assert_almost_equal(res1,1, decimal=4) + + # 2x2 matrix, max absolute eigenvalue is not unique and initial vector chosen for non-convergence - with self.assertWarns(Warning): - M1=numpy.array([[2.,1.], [0.,-2.]]) - M1op = MatrixOperator(M1) - res1 = M1op.PowerMethod(M1op,100, initial=DataContainer(numpy.array([1.,1.]))) + + M1=numpy.array([[2.,1.], [0.,-2.]]) + M1op = MatrixOperator(M1) + _,_,_,_,convergence = M1op.PowerMethod(M1op,100, initial=DataContainer(numpy.array([1.,1.])), return_all=True) + numpy.testing.assert_equal(convergence,False) # 2x2 matrix, max absolute eigenvalue is not unique and initial vector chosen for convergence M1=numpy.array([[2.,1.,0.],[0.,1.,1.], [0.,0.,1.]]) M1op = MatrixOperator(M1) - res1 = M1op.PowerMethod(M1op,100) + res1,_,_,_,convergence = M1op.PowerMethod(M1op,100, return_all=True) numpy.testing.assert_almost_equal(res1,2., decimal=4) - + numpy.testing.assert_equal(convergence,True) # Gradient Operator (float) ig = ImageGeometry(30,30) From e5cc8c19a8313904f861c566e5d401c1ba4b9656 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Mon, 18 Sep 2023 15:25:14 +0000 Subject: [PATCH 08/18] Fixed the iterate bug --- Wrappers/Python/cil/optimisation/operators/Operator.py | 4 +--- Wrappers/Python/test/test_BlockOperator.py | 3 +-- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index 8da95c6e73..8944441809 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -225,11 +225,9 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur diff = numpy.finfo('d').max i = 0 while (i < max_iteration and diff > tolerance): - i+=1 operator.direct(x0, out = y_tmp) - - + if square: #swap datacontainer references tmp = x0 diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py index bff62cd5af..adfac71572 100644 --- a/Wrappers/Python/test/test_BlockOperator.py +++ b/Wrappers/Python/test/test_BlockOperator.py @@ -187,9 +187,8 @@ def test_FiniteDiffOperator(self): G = FiniteDifferenceOperator(ig, direction=0, bnd_cond = 'Neumann') logging.info("{} {}".format(type(u), str(u.as_array()))) logging.info(str(G.direct(u).as_array())) - LinearOperator.PowerMethod(G, range_is_domain=False) # Gradient Operator norm, for one direction should be close to 2 - numpy.testing.assert_allclose(LinearOperator.PowerMethod(G, range_is_domain=False), numpy.sqrt(4), atol=0.1) + numpy.testing.assert_allclose(G.norm(), numpy.sqrt(4), atol=0.1) M1, N1, K1 = 200, 300, 2 ig1 = ImageGeometry(voxel_num_x = M1, voxel_num_y = N1, channels = K1) From 7e7fa3d8997a1d8d029a2fff6e406c35662332b5 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Mon, 18 Sep 2023 15:41:55 +0000 Subject: [PATCH 09/18] Formatting issues --- .../cil/optimisation/operators/Operator.py | 275 +++++++++--------- 1 file changed, 138 insertions(+), 137 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index 8944441809..05a2d1fc46 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -22,6 +22,7 @@ import functools import warnings + class Operator(object): """ Operator that maps from a space X -> Y @@ -45,11 +46,11 @@ def __init__(self, domain_geometry, **kwargs): def is_linear(self): '''Returns if the operator is linear''' return False - def direct(self,x, out=None): + + def direct(self, x, out=None): '''Returns the application of the Operator on x''' raise NotImplementedError - def norm(self, **kwargs): '''Returns the norm of the Operator. On first call the norm will be calculated using the operator's calculate_norm method. Subsequent calls will return the cached norm. @@ -68,7 +69,7 @@ def norm(self, **kwargs): return self._norm - def set_norm(self,norm=None): + def set_norm(self, norm=None): '''Sets the norm of the operator to a custom value. ''' self._norm = norm @@ -76,44 +77,49 @@ def set_norm(self,norm=None): def calculate_norm(self): '''Calculates the norm of the Operator''' raise NotImplementedError + def range_geometry(self): '''Returns the range of the Operator: Y space''' return self._range_geometry + def domain_geometry(self): '''Returns the domain of the Operator: X space''' return self._domain_geometry + @property def domain(self): return self.domain_geometry() + @property def range(self): return self.range_geometry() + def __rmul__(self, scalar): '''Defines the multiplication by a scalar on the left returns a ScaledOperator''' return ScaledOperator(self, scalar) - + def compose(self, *other, **kwargs): - # TODO: check equality of domain and range of operators - #if self.operator2.range_geometry != self.operator1.domain_geometry: - # raise ValueError('Cannot compose operators, check domain geometry of {} and range geometry of {}'.format(self.operato1,self.operator2)) - - return CompositionOperator(self, *other, **kwargs) + # TODO: check equality of domain and range of operators + # if self.operator2.range_geometry != self.operator1.domain_geometry: + # raise ValueError('Cannot compose operators, check domain geometry of {} and range geometry of {}'.format(self.operato1,self.operator2)) + + return CompositionOperator(self, *other, **kwargs) def __add__(self, other): return SumOperator(self, other) def __mul__(self, scalar): - return self.__rmul__(scalar) - + return self.__rmul__(scalar) + def __neg__(self): """ Return -self """ - return -1 * self - + return -1 * self + def __sub__(self, other): """ Returns the subtraction of the operators.""" - return self + (-1) * other + return self + (-1) * other class LinearOperator(Operator): @@ -129,28 +135,30 @@ class LinearOperator(Operator): range_geometry : ImageGeometry or AcquisitionGeometry, optional, default None range of the operator """ + def __init__(self, domain_geometry, **kwargs): super(LinearOperator, self).__init__(domain_geometry, **kwargs) + def is_linear(self): '''Returns if the operator is linear''' return True - def adjoint(self,x, out=None): + + def adjoint(self, x, out=None): '''returns the adjoint/inverse operation - + only available to linear operators''' raise NotImplementedError - - @staticmethod - def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False,range_is_domain=None ): + @staticmethod + def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False, range_is_domain=None): r"""Power method or Power iteration algorithm - + The Power method computes the largest (dominant) eigenvalue of a square matrix in magnitude, e.g., absolute value in the real case and modulus in the complex case. If the matrix is not square. The algorithm computes the largest (dominant) eigenvalue of :math: A^{T}*A :math:, returning the square root of this value. - - + + Parameters ---------- @@ -161,12 +169,12 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur Starting point for the Power method. tolerance: positive:`float`, default = 1e-5 Stopping criterion for the Power method. Check if two consecutive eigenvalue evaluations are below the tolerance. - return_all: boolean, default = False + return_all: `boolean`, default = False Toggles the verbosity of the return range_is_domain: `boolean`, default None Set this to `True` to apply the power method directly on the operator, :math: A :math:, and `False` to apply the power method to :math:A^TA:math: before taking the square root of the result. Leave as default `None` to determine this from domain and range geometry. - + Returns ------- @@ -195,18 +203,18 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur 2.0005647295658866 """ - convergence_check=True + convergence_check = True if range_is_domain is None: square = False try: - if operator.domain_geometry()==operator.range_geometry(): + if operator.domain_geometry() == operator.range_geometry(): square = True except AssertionError: # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972 pass else: - square=range_is_domain - + square = range_is_domain + if initial is None: x0 = operator.domain_geometry().allocate('random') else: @@ -225,60 +233,52 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur diff = numpy.finfo('d').max i = 0 while (i < max_iteration and diff > tolerance): - - operator.direct(x0, out = y_tmp) - - if square: - #swap datacontainer references + operator.direct(x0, out=y_tmp) + + if square: + # swap datacontainer references tmp = x0 x0 = y_tmp y_tmp = tmp else: - operator.adjoint(y_tmp,out=x0) - - + operator.adjoint(y_tmp, out=x0) + # Get eigenvalue using Rayleigh quotient: denominator=1, due to normalization - x0_norm=x0.norm() - if x0_norm Date: Tue, 26 Sep 2023 11:14:38 +0000 Subject: [PATCH 10/18] Changed flag name and ensured norm was calculated on A^TA --- .../cil/optimisation/operators/Operator.py | 42 ++++++++++++------- 1 file changed, 26 insertions(+), 16 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index 05a2d1fc46..e3145b4f77 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -150,13 +150,25 @@ def adjoint(self, x, out=None): raise NotImplementedError @staticmethod - def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False, range_is_domain=None): + def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False, compose_with_adjoint=None): r"""Power method or Power iteration algorithm - The Power method computes the largest (dominant) eigenvalue of a square matrix in magnitude, e.g., - absolute value in the real case and modulus in the complex case. + Fpr an operator, :math:`A`, the power method computes the iterations + .. math:: + x_{k+1} = A (\frac{x_k}{\|x_{k}\|}) \; (*) + + initialised with a random vector :math:`x_0` and returning the largest (dominant) eigenvalue in magnitude given by :math:`\|x_k\|`. + + If the flag, `compose_with_adjoint` is set to `True`, the algorithm computes the largest (dominant) eigenvalue of :math:`A^{T}A` + returning the square root of this value, i.e. the iterations: + .. math:: + x_{k+1} = A^TA (\frac{x_k}{\|x_{k}\|}) + + and returning :math:`\sqrt{\|x_k\|}`. - If the matrix is not square. The algorithm computes the largest (dominant) eigenvalue of :math: A^{T}*A :math:, returning the square root of this value. + If the flag `compose_with_adjoint` is set to `False`, the iterates in equation :math:`*` are used. In the default case, `compose_with_adjoint=None`, + the code checks to see if the domain and range of the operator are equal and sets `compose_with_adjoint=True` in the case they aren't and + `compose_with_adjoint=False` if they are. Parameters @@ -171,8 +183,8 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur Stopping criterion for the Power method. Check if two consecutive eigenvalue evaluations are below the tolerance. return_all: `boolean`, default = False Toggles the verbosity of the return - range_is_domain: `boolean`, default None - Set this to `True` to apply the power method directly on the operator, :math: A :math:, and `False` to apply the power method to :math:A^TA:math: before taking the square root of the result. + compose_with_adjoint: `boolean`, default None + Set this to `False` to apply the power method directly on the operator, :math: A :math:, and `True` to apply the power method to :math:A^TA:math: before taking the square root of the result. Leave as default `None` to determine this from domain and range geometry. @@ -185,7 +197,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur Corresponding eigenvector of the dominant eigenvalue. Only returned if return_all is True. list of eigenvalues: :obj:`list` List of eigenvalues. Only returned if return_all is True. - Convergence: `boolean' + convergence: `boolean' Check on wether the difference between the last two iterations is less than tolerance. Only returned if return_all is True. Examples @@ -204,16 +216,14 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur """ convergence_check = True - if range_is_domain is None: - square = False + if compose_with_adjoint is None: try: if operator.domain_geometry() == operator.range_geometry(): - square = True + compose_with_adjoint= False except AssertionError: # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972 - pass - else: - square = range_is_domain + compose_with_adjoint=True + if initial is None: x0 = operator.domain_geometry().allocate('random') @@ -235,7 +245,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur while (i < max_iteration and diff > tolerance): operator.direct(x0, out=y_tmp) - if square: + if not compose_with_adjoint: # swap datacontainer references tmp = x0 x0 = y_tmp @@ -253,7 +263,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur x0 /= x0_norm eig_new = numpy.abs(x0_norm) - if not square: + if compose_with_adjoint: eig_new = numpy.sqrt(eig_new) diff = numpy.abs(eig_new - eig_old) eig_list.append(eig_new) @@ -271,7 +281,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur def calculate_norm(self): r""" Returns the norm of the LinearOperator calculated by the PowerMethod with default values. """ - return LinearOperator.PowerMethod(self) + return LinearOperator.PowerMethod(self, compose_with_adjoint=True) @staticmethod def dot_test(operator, domain_init=None, range_init=None, tolerance=1e-6, **kwargs): From 439567f83fa6bd7e602042e97fbafc452d0384a9 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Tue, 26 Sep 2023 12:57:17 +0000 Subject: [PATCH 11/18] Debugging and tests --- Wrappers/Python/cil/optimisation/operators/Operator.py | 3 +++ Wrappers/Python/test/test_Operator.py | 9 +++++---- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index e3145b4f77..cef879e930 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -216,10 +216,13 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur """ convergence_check = True + if compose_with_adjoint is None: try: if operator.domain_geometry() == operator.range_geometry(): compose_with_adjoint= False + else: + compose_with_adjoint=True except AssertionError: # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972 compose_with_adjoint=True diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index ea756157d6..7117dc23a7 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -284,6 +284,7 @@ def test_PowerMethod(self): # Test with the norm res2 = M1op.norm() + res1 = M1op.PowerMethod(M1op,100, compose_with_adjoint=True) numpy.testing.assert_almost_equal(res1,res2, decimal=4) @@ -302,8 +303,8 @@ def test_PowerMethod(self): # 3x3 complex matrix, (real+complex eigenvalue), dominant eigenvalue = 3.1624439599276974 M1 = numpy.array([[2,0,0],[1,2j,1j],[3, 3-1j,3]]) M1op = MatrixOperator(M1) - res1 = M1op.PowerMethod(M1op,100) - numpy.testing.assert_almost_equal(res1,3.1624439599276974, decimal=4) + res1 = M1op.PowerMethod(M1op,120) + numpy.testing.assert_almost_equal(res1,3.1624439599276974, decimal=3) # 2x2 non-diagonalisable nilpotent matrix M1=numpy.array([[0.,1.], [0.,0.]]) @@ -311,10 +312,10 @@ def test_PowerMethod(self): res1 = M1op.PowerMethod(M1op,5) numpy.testing.assert_almost_equal(res1,0, decimal=4) - # 2x2 non-diagonalisable nilpotent matrix where range_is_domain is False + # 2x2 non-diagonalisable nilpotent matrix where compose_with_adjoint=True M1=numpy.array([[0.,1.], [0.,0.]]) M1op = MatrixOperator(M1) - res1 = M1op.PowerMethod(M1op,5, range_is_domain=False) + res1 = M1op.PowerMethod(M1op,5, compose_with_adjoint=True) numpy.testing.assert_almost_equal(res1,1, decimal=4) From f85c98c081cbbd51a74d6f717cd3d6c0fb038008 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Tue, 26 Sep 2023 13:06:03 +0000 Subject: [PATCH 12/18] Updated CHANGELOG [skip ci] --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0c6794ead4..ba01af88fc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,7 +9,7 @@ - Fix for show_geometry bug for 2D data - Added warmstart capability to proximal evaluation of the CIL TotalVariation function. - FBP split processing bug fix - now respects panel origin - + - Bug fix in the LinearOperator norm with an additional flag for the algorithm linearOperator.PowerMethod * 23.0.1 - Fix bug with NikonReader requiring ROI to be set in constructor. From 47cd42c4c500e689d1a8d1d57ceda2e4d52a373c Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Tue, 26 Sep 2023 15:41:18 +0000 Subject: [PATCH 13/18] Naming changes - discussed with Gemma --- .../cil/optimisation/operators/Operator.py | 36 ++++++++++--------- Wrappers/Python/test/test_BlockOperator.py | 1 - Wrappers/Python/test/test_Operator.py | 25 ++++++++++--- 3 files changed, 41 insertions(+), 21 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index 800d8d137d..753309dd93 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -150,25 +150,26 @@ def adjoint(self, x, out=None): raise NotImplementedError @staticmethod - def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False, compose_with_adjoint=None): + def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False, method='auto'): r"""Power method or Power iteration algorithm - Fpr an operator, :math:`A`, the power method computes the iterations + The power method contains two different algorithms chosen by the `method` flag. + + In the case `method="direct_only"`, for operator, :math:`A`, the power method computes the iterations .. math:: x_{k+1} = A (\frac{x_k}{\|x_{k}\|}) \; (*) initialised with a random vector :math:`x_0` and returning the largest (dominant) eigenvalue in magnitude given by :math:`\|x_k\|`. - If the flag, `compose_with_adjoint` is set to `True`, the algorithm computes the largest (dominant) eigenvalue of :math:`A^{T}A` + In the case `method="composed_with_adjoint"`, the algorithm computes the largest (dominant) eigenvalue of :math:`A^{T}A` returning the square root of this value, i.e. the iterations: .. math:: - x_{k+1} = A^TA (\frac{x_k}{\|x_{k}\|}) + x_{k+1} = A^*A (\frac{x_k}{\|x_{k}\|}) and returning :math:`\sqrt{\|x_k\|}`. - If the flag `compose_with_adjoint` is set to `False`, the iterates in equation :math:`*` are used. In the default case, `compose_with_adjoint=None`, - the code checks to see if the domain and range of the operator are equal and sets `compose_with_adjoint=True` in the case they aren't and - `compose_with_adjoint=False` if they are. + The default flag is `method="auto"`, the algorithm checks to see if the `operator.domain_geometry() == operator.range_geometry()` and if so + uses the method "direct_only" and if not the method "composed_with_adjoint". Parameters @@ -217,16 +218,19 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur """ convergence_check = True - if compose_with_adjoint is None: + if method not in ["auto","direct_only","composed_with_adjoint"]: + ValueError(" The flag 'method' should take one of 'auto','direct_only','composed_with_adjoint'") + + if method=="auto": try: if operator.domain_geometry() == operator.range_geometry(): - compose_with_adjoint= False + method="direct_only" else: - compose_with_adjoint=True + method="composed_with_adjoint" except AssertionError: # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972 - compose_with_adjoint=True - + method="composed_with_adjoint" + if initial is None: x0 = operator.domain_geometry().allocate('random') @@ -248,12 +252,12 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur while (i < max_iteration and diff > tolerance): operator.direct(x0, out=y_tmp) - if not compose_with_adjoint: + if method=="direct_only": # swap datacontainer references tmp = x0 x0 = y_tmp y_tmp = tmp - else: + elif method=="composed_with_adjoint": operator.adjoint(y_tmp, out=x0) # Get eigenvalue using Rayleigh quotient: denominator=1, due to normalization @@ -266,7 +270,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur x0 /= x0_norm eig_new = numpy.abs(x0_norm) - if compose_with_adjoint: + if method=="composed_with_adjoint": eig_new = numpy.sqrt(eig_new) diff = numpy.abs(eig_new - eig_old) eig_list.append(eig_new) @@ -284,7 +288,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur def calculate_norm(self): r""" Returns the norm of the LinearOperator calculated by the PowerMethod with default values. """ - return LinearOperator.PowerMethod(self, compose_with_adjoint=True) + return LinearOperator.PowerMethod(self, method="composed_with_adjoint") @staticmethod def dot_test(operator, domain_init=None, range_init=None, tolerance=1e-6, **kwargs): diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py index adfac71572..7f417bb382 100644 --- a/Wrappers/Python/test/test_BlockOperator.py +++ b/Wrappers/Python/test/test_BlockOperator.py @@ -26,7 +26,6 @@ from cil.framework import ImageGeometry, ImageData import numpy from cil.optimisation.operators import FiniteDifferenceOperator -from cil.optimisation.operators import LinearOperator initialise_tests() class TestBlockOperator(unittest.TestCase): diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index 15c4904711..6a96610fca 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -285,7 +285,7 @@ def test_PowerMethod(self): # Test with the norm res2 = M1op.norm() - res1 = M1op.PowerMethod(M1op,100, compose_with_adjoint=True) + res1 = M1op.PowerMethod(M1op,100, method="composed_with_adjoint") numpy.testing.assert_almost_equal(res1,res2, decimal=4) @@ -313,10 +313,10 @@ def test_PowerMethod(self): res1 = M1op.PowerMethod(M1op,5) numpy.testing.assert_almost_equal(res1,0, decimal=4) - # 2x2 non-diagonalisable nilpotent matrix where compose_with_adjoint=True + # 2x2 non-diagonalisable nilpotent matrix where method="composed_with_adjoint" M1=numpy.array([[0.,1.], [0.,0.]]) M1op = MatrixOperator(M1) - res1 = M1op.PowerMethod(M1op,5, compose_with_adjoint=True) + res1 = M1op.PowerMethod(M1op,5, method="composed_with_adjoint") numpy.testing.assert_almost_equal(res1,1, decimal=4) @@ -350,7 +350,13 @@ def test_PowerMethod(self): # Identity Operator Id = IdentityOperator(ig) res1 = Id.PowerMethod(Id,100) - numpy.testing.assert_almost_equal(res1,1.0, decimal=4) + numpy.testing.assert_almost_equal(res1,1.0, decimal=4) + + # Test errors produced if not a valid method + try: + res1 = Id.PowerMethod(Id,100, method='gobledigook') + except ValueError: + pass def test_Norm(self): @@ -377,6 +383,17 @@ def test_Norm(self): #recalculates norm self.assertAlmostEqual(G.norm(), numpy.sqrt(8), 2) + # 2x2 real matrix, dominant eigenvalue = 2. Check norm uses the right flag for power method + M1 = numpy.array([[1,0],[1,2]], dtype=float) + M1op = MatrixOperator(M1) + res1 = M1op.norm() + res2 = M1op.PowerMethod(M1op,100) + res3 = M1op.PowerMethod(M1op,100, method="composed_with_adjoint") + res4 = M1op.PowerMethod(M1op,100, method="direct_only") + numpy.testing.assert_almost_equal(res1,res3, decimal=4) + self.assertNotEqual(res1, res2) + self.assertNotEqual(res1,res4) + def test_ProjectionMap(self): # Check if direct is correct From 0fcfeddf4540149bca504b98b41f68ba49c44aae Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Wed, 27 Sep 2023 10:56:37 +0000 Subject: [PATCH 14/18] Reponse to Gemma's review --- .../cil/optimisation/operators/Operator.py | 88 +++++++++++-------- 1 file changed, 50 insertions(+), 38 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index 753309dd93..ab99f0675a 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -153,25 +153,9 @@ def adjoint(self, x, out=None): def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False, method='auto'): r"""Power method or Power iteration algorithm - The power method contains two different algorithms chosen by the `method` flag. - - In the case `method="direct_only"`, for operator, :math:`A`, the power method computes the iterations - .. math:: - x_{k+1} = A (\frac{x_k}{\|x_{k}\|}) \; (*) - - initialised with a random vector :math:`x_0` and returning the largest (dominant) eigenvalue in magnitude given by :math:`\|x_k\|`. - - In the case `method="composed_with_adjoint"`, the algorithm computes the largest (dominant) eigenvalue of :math:`A^{T}A` - returning the square root of this value, i.e. the iterations: - .. math:: - x_{k+1} = A^*A (\frac{x_k}{\|x_{k}\|}) - - and returning :math:`\sqrt{\|x_k\|}`. - - The default flag is `method="auto"`, the algorithm checks to see if the `operator.domain_geometry() == operator.range_geometry()` and if so - uses the method "direct_only" and if not the method "composed_with_adjoint". - - + The Power method computes the largest (dominant) eigenvalue of a matrix in magnitude, e.g., + absolute value in the real case and modulus in the complex case. + Parameters ---------- @@ -184,10 +168,9 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur Stopping criterion for the Power method. Check if two consecutive eigenvalue evaluations are below the tolerance. return_all: `boolean`, default = False Toggles the verbosity of the return - compose_with_adjoint: `boolean`, default None - Set this to `False` to apply the power method directly on the operator, :math: A :math:, and `True` to apply the power method to :math:A^TA:math: before taking the square root of the result. - Leave as default `None` to determine this from domain and range geometry. - + method: `string` one of `auto`, `composed_with_adjoint` and `direct_only`, default = `auto` + The default `auto` lets the code choose the method, this can be specified with `direct_only` or `compose_with_adjoint` + Returns ------- @@ -200,6 +183,27 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur List of eigenvalues. Only returned if return_all is True. convergence: `boolean' Check on wether the difference between the last two iterations is less than tolerance. Only returned if return_all is True. + + + Note + ----- + The power method contains two different algorithms chosen by the `method` flag. + + In the case `method="direct_only"`, for operator, :math:`A`, the power method computes the iterations + .. math:: + x_{k+1} = A (\frac{x_k}{\|x_{k}\|}) \; (*) + + initialised with a random vector :math:`x_0` and returning the largest (dominant) eigenvalue in magnitude given by :math:`\|x_k\|`. + + In the case `method="composed_with_adjoint"`, the algorithm computes the largest (dominant) eigenvalue of :math:`A^{T}A` + returning the square root of this value, i.e. the iterations: + .. math:: + x_{k+1} = A^*A (\frac{x_k}{\|x_{k}\|}) + + and returning :math:`\sqrt{\|x_k\|}`. + + The default flag is `method="auto"`, the algorithm checks to see if the `operator.domain_geometry() == operator.range_geometry()` and if so + uses the method "direct_only" and if not the method "composed_with_adjoint". Examples -------- @@ -216,20 +220,26 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur 2.0005647295658866 """ - convergence_check = True - if method not in ["auto","direct_only","composed_with_adjoint"]: - ValueError(" The flag 'method' should take one of 'auto','direct_only','composed_with_adjoint'") + + allowed_methods = ["auto","direct_only","composed_with_adjoint"] + + if method not in allowed_methods: + raise ValueError("The argument 'method' can be set to one of {0} got {1}".format(allowed_methods, method)) + apply_adjoint=True + if method == "direct_only": + apply_adjoint=False if method=="auto": try: - if operator.domain_geometry() == operator.range_geometry(): - method="direct_only" - else: - method="composed_with_adjoint" + geometries_match = operator.domain_geometry() == operator.range_geometry() + except AssertionError: # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972 - method="composed_with_adjoint" + pass + else: + if geometries_match: + apply_adjoint=False if initial is None: @@ -245,19 +255,20 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur # initial guess for dominant eigenvalue eig_old = 1. - - eig_list = [] + if return_all: + eig_list = [] + convergence_check = True diff = numpy.finfo('d').max i = 0 while (i < max_iteration and diff > tolerance): operator.direct(x0, out=y_tmp) - if method=="direct_only": + if not apply_adjoint: # swap datacontainer references tmp = x0 x0 = y_tmp y_tmp = tmp - elif method=="composed_with_adjoint": + else: operator.adjoint(y_tmp, out=x0) # Get eigenvalue using Rayleigh quotient: denominator=1, due to normalization @@ -270,14 +281,15 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur x0 /= x0_norm eig_new = numpy.abs(x0_norm) - if method=="composed_with_adjoint": + if apply_adjoint: eig_new = numpy.sqrt(eig_new) diff = numpy.abs(eig_new - eig_old) - eig_list.append(eig_new) + if return_all: + eig_list.append(eig_new) eig_old = eig_new i += 1 - if i == max_iteration: + if return_all and i == max_iteration: convergence_check = False if return_all: From 39b1c62f4c0daf2637339380e8bdcfeb2127d549 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Wed, 27 Sep 2023 11:02:19 +0000 Subject: [PATCH 15/18] Documentation changes --- .../cil/optimisation/operators/Operator.py | 38 +++++++++---------- 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index ab99f0675a..b78ef4581a 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -155,7 +155,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur The Power method computes the largest (dominant) eigenvalue of a matrix in magnitude, e.g., absolute value in the real case and modulus in the complex case. - + Parameters ---------- @@ -168,9 +168,9 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur Stopping criterion for the Power method. Check if two consecutive eigenvalue evaluations are below the tolerance. return_all: `boolean`, default = False Toggles the verbosity of the return - method: `string` one of `auto`, `composed_with_adjoint` and `direct_only`, default = `auto` - The default `auto` lets the code choose the method, this can be specified with `direct_only` or `compose_with_adjoint` - + method: `string` one of `"auto"`, `"composed_with_adjoint"` and `"direct_only"`, default = `"auto"` + The default `auto` lets the code choose the method, this can be specified with `"direct_only"` or `"composed_with_adjoint"` + Returns ------- @@ -183,8 +183,8 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur List of eigenvalues. Only returned if return_all is True. convergence: `boolean' Check on wether the difference between the last two iterations is less than tolerance. Only returned if return_all is True. - - + + Note ----- The power method contains two different algorithms chosen by the `method` flag. @@ -220,27 +220,26 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur 2.0005647295658866 """ - - - allowed_methods = ["auto","direct_only","composed_with_adjoint"] - + + allowed_methods = ["auto", "direct_only", "composed_with_adjoint"] + if method not in allowed_methods: - raise ValueError("The argument 'method' can be set to one of {0} got {1}".format(allowed_methods, method)) + raise ValueError("The argument 'method' can be set to one of {0} got {1}".format( + allowed_methods, method)) - apply_adjoint=True + apply_adjoint = True if method == "direct_only": - apply_adjoint=False - if method=="auto": + apply_adjoint = False + if method == "auto": try: geometries_match = operator.domain_geometry() == operator.range_geometry() except AssertionError: # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972 - pass + pass else: if geometries_match: - apply_adjoint=False - + apply_adjoint = False if initial is None: x0 = operator.domain_geometry().allocate('random') @@ -594,13 +593,12 @@ def adjoint(self, x, out=None): raise ValueError('No adjoint operation with non-linear operators') def is_linear(self): - return self.linear_flag - + return self.linear_flag def calculate_norm(self): '''Returns the norm of the CompositionOperator, that is the product of the norms of its operators.''' norm = 1. for operator in self.operators: - norm *= operator.norm() + norm *= operator.norm() return norm From efa58b40d6ffcce05e4dd9e8eb388bb0f1deb321 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Wed, 27 Sep 2023 12:54:07 +0000 Subject: [PATCH 16/18] Changed warnings to logging.warning --- Wrappers/Python/cil/optimisation/operators/Operator.py | 6 +++--- Wrappers/Python/test/test_BlockOperator.py | 1 + 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index b78ef4581a..a0047ecd0f 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -20,7 +20,7 @@ from numbers import Number import numpy import functools -import warnings +import logging class Operator(object): @@ -61,7 +61,7 @@ def norm(self, **kwargs): ''' if len(kwargs) != 0: - warnings.warn('norm: the norm method does not use any parameters.\n\ + logging.warning('norm: the norm method does not use any parameters.\n\ For LinearOperators you can use PowerMethod to calculate the norm with non-default parameters and use set_norm to set it') if self._norm is None: @@ -273,7 +273,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur # Get eigenvalue using Rayleigh quotient: denominator=1, due to normalization x0_norm = x0.norm() if x0_norm < tolerance: - Warning( + logging.warning( 'The operator has at least one zero eigenvector and is likely to be nilpotent') eig_new = 0. break diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py index 7f417bb382..93e5bbd3e2 100644 --- a/Wrappers/Python/test/test_BlockOperator.py +++ b/Wrappers/Python/test/test_BlockOperator.py @@ -26,6 +26,7 @@ from cil.framework import ImageGeometry, ImageData import numpy from cil.optimisation.operators import FiniteDifferenceOperator + initialise_tests() class TestBlockOperator(unittest.TestCase): From dbf21524e44b5631347e87f08b2dd7b65889d2f6 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Thu, 28 Sep 2023 08:52:32 +0000 Subject: [PATCH 17/18] Documentation changes --- .../Python/cil/optimisation/operators/Operator.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index a0047ecd0f..df2b6e643b 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -190,17 +190,11 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur The power method contains two different algorithms chosen by the `method` flag. In the case `method="direct_only"`, for operator, :math:`A`, the power method computes the iterations - .. math:: - x_{k+1} = A (\frac{x_k}{\|x_{k}\|}) \; (*) - - initialised with a random vector :math:`x_0` and returning the largest (dominant) eigenvalue in magnitude given by :math:`\|x_k\|`. + :math:`x_{k+1} = A (x_k/\|x_{k}\|)` initialised with a random vector :math:`x_0` and returning the largest (dominant) eigenvalue in magnitude given by :math:`\|x_k\|`. In the case `method="composed_with_adjoint"`, the algorithm computes the largest (dominant) eigenvalue of :math:`A^{T}A` - returning the square root of this value, i.e. the iterations: - .. math:: - x_{k+1} = A^*A (\frac{x_k}{\|x_{k}\|}) - - and returning :math:`\sqrt{\|x_k\|}`. + returning the square root of this value, i.e. the iterations: + :math:`x_{k+1} = A^TA (x_k/\|x_{k}\|)` and returning :math:`\sqrt{\|x_k\|}`. The default flag is `method="auto"`, the algorithm checks to see if the `operator.domain_geometry() == operator.range_geometry()` and if so uses the method "direct_only" and if not the method "composed_with_adjoint". From 4fb1099a17ea35d5af583f4fe47763062054dce8 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Thu, 28 Sep 2023 09:58:08 +0000 Subject: [PATCH 18/18] Minor documentation change --- Wrappers/Python/cil/optimisation/operators/Operator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py index df2b6e643b..adf89b9fba 100644 --- a/Wrappers/Python/cil/optimisation/operators/Operator.py +++ b/Wrappers/Python/cil/optimisation/operators/Operator.py @@ -181,7 +181,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur Corresponding eigenvector of the dominant eigenvalue. Only returned if return_all is True. list of eigenvalues: :obj:`list` List of eigenvalues. Only returned if return_all is True. - convergence: `boolean' + convergence: `boolean` Check on wether the difference between the last two iterations is less than tolerance. Only returned if return_all is True.