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

Random seed for power method #1585

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
+18 −5
Open

Random seed for power method #1585

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
68eee1b
Initial solution
MargaretDuff Nov 22, 2023
fe49005
Merge branch 'master' into powermethod-random
MargaretDuff Nov 22, 2023
b505be6
Merge branch 'master' into powermethod-random
MargaretDuff Nov 24, 2023
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
10 changes: 6 additions & 4 deletions Wrappers/Python/cil/optimisation/operators/Operator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,7 @@ class LinearOperator(Operator):
"""

def __init__(self, domain_geometry, **kwargs):
self._random_number_seed=kwargs.get('seed', 0)
super(LinearOperator, self).__init__(domain_geometry, **kwargs)

def is_linear(self):
Expand All @@ -171,7 +172,7 @@ def adjoint(self, x, out=None):
raise NotImplementedError

@staticmethod
def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False, method='auto'):
def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, return_all=False, method='auto', seed=None):
r"""Power method or Power iteration algorithm

The Power method computes the largest (dominant) eigenvalue of a matrix in magnitude, e.g.,
Expand All @@ -191,7 +192,8 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur
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 `"composed_with_adjoint"`

seed: `int` or default = None
The random seed used by the random number generator to create the `initial` starting point, if it is not provided.

Returns
-------
Expand Down Expand Up @@ -257,7 +259,7 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5, retur
apply_adjoint = False

if initial is None:
x0 = operator.domain_geometry().allocate('random')
x0 = operator.domain_geometry().allocate('random', seed=seed)
else:
x0 = initial.copy()

Expand Down Expand Up @@ -314,7 +316,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, method="composed_with_adjoint")
return LinearOperator.PowerMethod(self, method="composed_with_adjoint", seed=self._random_number_seed)

@staticmethod
def dot_test(operator, domain_init=None, range_init=None, tolerance=1e-6, **kwargs):
Expand Down
13 changes: 12 additions & 1 deletion Wrappers/Python/test/test_Operator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -288,7 +288,18 @@ def test_PowerMethod(self):
res1 = M1op.PowerMethod(M1op,100, method="composed_with_adjoint")
numpy.testing.assert_almost_equal(res1,res2, decimal=4)


#Test random seed
numpy.testing.assert_equal(M1op._random_number_seed,0)
res1 = M1op.PowerMethod(M1op,1, method="composed_with_adjoint", seed=2)
res2 = M1op.PowerMethod(M1op,1, method="composed_with_adjoint", seed=2)
res3 = M1op.PowerMethod(M1op,1, method="composed_with_adjoint", seed=3)
numpy.testing.assert_almost_equal(res1,res2, decimal=4)
numpy.testing.assert_raises(AssertionError, numpy.testing.assert_almost_equal, res1,res3)
res4=M1op.norm()
res5=M1op.norm()
numpy.testing.assert_equal(res4,res5)


# 2x3 real matrix, dominant eigenvalue = 4.711479432297657
M1 = numpy.array([[1.,0.,3],[1,2.,3]])
M1op = MatrixOperator(M1)
Expand Down