From c633a3047c1583f20ea0aa86d2f23b117912144e Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Sat, 29 Jan 2022 14:50:20 +0100 Subject: [PATCH 1/4] Correct the norm update for Power Method x_new should be divided by its norm, not by x_old_norm. --- modopt/math/matrix.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/modopt/math/matrix.py b/modopt/math/matrix.py index 939cf41f..a496ef13 100644 --- a/modopt/math/matrix.py +++ b/modopt/math/matrix.py @@ -348,17 +348,21 @@ def get_spec_rad(self, tolerance=1e-6, max_iter=20, extra_factor=1.0): # Set (or reset) values of x. x_old = self._set_initial_x() + xp = get_array_module(x_old) + x_old_norm = xp.linalg.norm(x_old) + + x_old = x_old / x_old_norm + # Iterate until the L2 norm of x converges. for i_elem in range(max_iter): - xp = get_array_module(x_old) - x_old_norm = xp.linalg.norm(x_old) - - x_new = self._operator(x_old) / x_old_norm + x_new = self._operator(x_old) x_new_norm = xp.linalg.norm(x_new) + x_new = x_new / x_new_norm + if (xp.abs(x_new_norm - x_old_norm) < tolerance): message = ( ' - Power Method converged after {0} iterations!' @@ -374,6 +378,7 @@ def get_spec_rad(self, tolerance=1e-6, max_iter=20, extra_factor=1.0): print(message.format(max_iter)) xp.copyto(x_old, x_new) + x_old_norm = x_new_norm self.spec_rad = x_new_norm * extra_factor self.inv_spec_rad = 1.0 / self.spec_rad From 6459291b3be2f1bddaeb0904a746f68bafd8dca7 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Sat, 29 Jan 2022 16:03:13 +0100 Subject: [PATCH 2/4] fix test value We are testing for eigen value of Identity. It should be one. --- modopt/tests/test_math.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/modopt/tests/test_math.py b/modopt/tests/test_math.py index 99908e02..b3720fed 100644 --- a/modopt/tests/test_math.py +++ b/modopt/tests/test_math.py @@ -234,13 +234,13 @@ def test_powermethod_converged(self): """Test PowerMethod converged.""" npt.assert_almost_equal( self.pmInstance1.spec_rad, - 0.90429242629600837, + 1.0, err_msg='Incorrect spectral radius: converged', ) npt.assert_almost_equal( self.pmInstance1.inv_spec_rad, - 1.1058369736612865, + 1.0, err_msg='Incorrect inverse spectral radius: converged', ) @@ -248,13 +248,13 @@ def test_powermethod_unconverged(self): """Test PowerMethod unconverged.""" npt.assert_almost_equal( self.pmInstance2.spec_rad, - 0.92048833577059219, + 1.0, err_msg='Incorrect spectral radius: unconverged', ) npt.assert_almost_equal( self.pmInstance2.inv_spec_rad, - 1.0863798715741946, + 1.0, err_msg='Incorrect inverse spectral radius: unconverged', ) From b271d30c962afeb1058f767453ad2a8e27af77d8 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Sat, 29 Jan 2022 16:05:52 +0100 Subject: [PATCH 3/4] fix WPS350 --- modopt/math/matrix.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modopt/math/matrix.py b/modopt/math/matrix.py index a496ef13..6bd225f6 100644 --- a/modopt/math/matrix.py +++ b/modopt/math/matrix.py @@ -351,7 +351,7 @@ def get_spec_rad(self, tolerance=1e-6, max_iter=20, extra_factor=1.0): xp = get_array_module(x_old) x_old_norm = xp.linalg.norm(x_old) - x_old = x_old / x_old_norm + x_old /= x_old_norm # Iterate until the L2 norm of x converges. for i_elem in range(max_iter): @@ -361,7 +361,7 @@ def get_spec_rad(self, tolerance=1e-6, max_iter=20, extra_factor=1.0): x_new_norm = xp.linalg.norm(x_new) - x_new = x_new / x_new_norm + x_new /= x_new_norm if (xp.abs(x_new_norm - x_old_norm) < tolerance): message = ( From d7d917f4788d3c889dc8feaec0aaa42c1ef4804e Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Sun, 30 Jan 2022 15:09:10 +0100 Subject: [PATCH 4/4] fix test value for unconverged case --- modopt/math/matrix.py | 4 ++-- modopt/tests/test_math.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/modopt/math/matrix.py b/modopt/math/matrix.py index 6bd225f6..8361531d 100644 --- a/modopt/math/matrix.py +++ b/modopt/math/matrix.py @@ -285,9 +285,9 @@ class PowerMethod(object): >>> np.random.seed(1) >>> pm = PowerMethod(lambda x: x.dot(x.T), (3, 3)) >>> np.around(pm.spec_rad, 6) - 0.904292 + 1.0 >>> np.around(pm.inv_spec_rad, 6) - 1.105837 + 1.0 Notes ----- diff --git a/modopt/tests/test_math.py b/modopt/tests/test_math.py index b3720fed..ba175ae6 100644 --- a/modopt/tests/test_math.py +++ b/modopt/tests/test_math.py @@ -248,13 +248,13 @@ def test_powermethod_unconverged(self): """Test PowerMethod unconverged.""" npt.assert_almost_equal( self.pmInstance2.spec_rad, - 1.0, + 0.8675467477372257, err_msg='Incorrect spectral radius: unconverged', ) npt.assert_almost_equal( self.pmInstance2.inv_spec_rad, - 1.0, + 1.152675636913221, err_msg='Incorrect inverse spectral radius: unconverged', )