Repository: jdmoorman/kaczmarz-algorithms

Issue:
Deal with inconsistent systems.
We have two options when the system of equations is inconsistent:
* Set a default finite value for `maxiter`.
* Raise a StopIteration whenever the error increases from one iteration to the next.
My preference is the former, since I frequently plot error vs iteration for inconsistent systems knowing that they will not converge.

Discussion:

Patch:
diff --git a/README.md b/README.md
index ef25c1e..e41acc2 100644
--- a/README.md
+++ b/README.md
@@ -111,7 +111,7 @@ For example, to implement a strategy which uses of the equations of your system
 ...         self.n_rows = A.shape[0]
 ...         self.row_index = None
 ...
-...     def select_row_index(self, xk):
+...     def _select_row_index(self, xk):
 ...         if self.row_index is None:
 ...             self.row_index = self.n_rows
 ...         self.row_index = (self.row_index - 1) % self.n_rows
diff --git a/setup.py b/setup.py
index d53d7c2..4334c63 100644
--- a/setup.py
+++ b/setup.py
@@ -8,7 +8,13 @@ from setuptools import find_packages, setup
 with open("README.md") as readme_file:
     readme = readme_file.read()
 
-test_requirements = ["pytest", "pytest-cov", "pytest-raises", "pytest-allclose"]
+test_requirements = [
+    "pytest",
+    "pytest-cov",
+    "pytest-raises",
+    "pytest-allclose",
+    "pytest-timeout",
+]
 
 docs_requirements = ["sphinx==1.8.5", "sphinx-rtd-theme", "m2r"]
 
diff --git a/src/kaczmarz/_abc.py b/src/kaczmarz/_abc.py
index bb65812..d655324 100644
--- a/src/kaczmarz/_abc.py
+++ b/src/kaczmarz/_abc.py
@@ -18,13 +18,14 @@ class Base(ABC):
         Tolerance for convergence, `norm(normalized_residual) <= tol`.
     maxiter : int or float, optional
         Maximum number of iterations.
+        Defaults to `log(tol^2)/log(1 - 1/(10*min(m, n)))`.
     callback : function, optional
         User-supplied function to call after each iteration.
         It is called as callback(xk), where xk is the current solution vector.
     """
 
     def __init__(
-        self, A, b, x0=None, tol=1e-5, maxiter=float("inf"), callback=None,
+        self, A, b, x0=None, tol=1e-5, maxiter=None, callback=None,
     ):
         row_norms = np.sqrt((A ** 2).sum(axis=1)).reshape(-1, 1)
         b = np.array(b)
@@ -42,6 +43,10 @@ class Base(ABC):
 
         self._x0 = np.array(x0, dtype="float64").ravel()
         self._tol = tol
+
+        if maxiter is None:
+            maxiter = 2 * np.log(tol) / np.log(1 - 1 / (10 * min(A.shape)))
+
         self._maxiter = maxiter
         if callback is None:
 
@@ -66,6 +71,37 @@ class Base(ABC):
         """
         return self._xk.copy().reshape(*self._iterate_shape)
 
+    @classmethod
+    def iterates(cls, *args, **kwargs):
+        """Get the Kaczmarz iterates.
+
+        TODO: Describe the invisible arguments.
+
+        Returns
+        -------
+        iterates : iterable((n,) or (n, 1) array)
+            An iterable of the Kaczmarz iterates.
+            The shapes will be inferred from the shape of `x0` if provided, or `b` otherwise.
+        """
+        return cls(*args, **kwargs)
+
+    @classmethod
+    def solve(cls, *args, **kwargs):
+        """Solve a linear system of equations using the Kaczmarz algorithm.
+
+        TODO: Describe the invisible arguments.
+
+        Returns
+        -------
+        x : (n,) or (n, 1) array
+            The solution to the system `Ax = b`.
+            The shape will be inferred from the shape of `x0` if provided, or `b` otherwise.
+        """
+        iterates = cls.iterates(*args, **kwargs)
+        for x in iterates:
+            pass
+        return x
+
     def __next__(self):
         """Perform an iteration of the Kaczmarz algorithm.
 
@@ -86,8 +122,8 @@ class Base(ABC):
             raise StopIteration
 
         self._k += 1
-        self._ik = self.select_row_index(self._xk)
-        self._xk = self.update_iterate(self._xk, self._ik)
+        self._ik = self._select_row_index(self._xk)
+        self._xk = self._update_iterate(self._xk, self._ik)
         self._callback(self.xk)
 
         return self.xk
@@ -96,7 +132,22 @@ class Base(ABC):
         """Iterator for iterates of the Kaczmarz algorithm."""
         return self
 
-    def update_iterate(self, xk, ik):
+    @abstractmethod
+    def _select_row_index(self, xk):
+        """Select a row to use for the next Kaczmarz update.
+
+        Parameters
+        ----------
+        xk : (n,) array
+            The current Kaczmarz iterate.
+
+        Returns
+        -------
+        ik : int
+            The index of the next row to use.
+        """
+
+    def _update_iterate(self, xk, ik):
         """Apply the Kaczmarz update.
 
         Parameters
@@ -140,49 +191,3 @@ class Base(ABC):
             return True
 
         return False
-
-    @classmethod
-    def iterates(cls, *args, **kwargs):
-        """Get the Kaczmarz iterates.
-
-        TODO: Describe the invisible arguments.
-
-        Returns
-        -------
-        iterates : iterable((n,) or (n, 1) array)
-            An iterable of the Kaczmarz iterates.
-            The shapes will be inferred from the shape of `x0` if provided, or `b` otherwise.
-        """
-        return cls(*args, **kwargs)
-
-    @classmethod
-    def solve(cls, *args, **kwargs):
-        """Solve a linear system of equations using the Kaczmarz algorithm.
-
-        TODO: Describe the invisible arguments.
-
-        Returns
-        -------
-        x : (n,) or (n, 1) array
-            The solution to the system `Ax = b`.
-            The shape will be inferred from the shape of `x0` if provided, or `b` otherwise.
-        """
-        iterates = cls.iterates(*args, **kwargs)
-        for x in iterates:
-            pass
-        return x
-
-    @abstractmethod
-    def select_row_index(self, xk):
-        """Select a row to use for the next Kaczmarz update.
-
-        Parameters
-        ----------
-        xk : (n,) array
-            The current Kaczmarz iterate.
-
-        Returns
-        -------
-        ik : int
-            The index of the next row to use.
-        """
diff --git a/src/kaczmarz/_variants.py b/src/kaczmarz/_variants.py
index 137ee80..a653cb8 100644
--- a/src/kaczmarz/_variants.py
+++ b/src/kaczmarz/_variants.py
@@ -21,7 +21,7 @@ class Cyclic(kaczmarz.Base):
         self.n_rows = A.shape[0]
         self.row_index = -1
 
-    def select_row_index(self, xk):
+    def _select_row_index(self, xk):
         """Get the next row index, in cyclic order.
 
         Parameters


Test patch:
diff --git a/tests/test_abc.py b/tests/test_abc.py
index 322fb6a..3438e9d 100644
--- a/tests/test_abc.py
+++ b/tests/test_abc.py
@@ -7,7 +7,7 @@ import kaczmarz
 @pytest.fixture()
 def DummyStrategy():
     class _DummyStrategy(kaczmarz.Base):
-        def select_row_index(self, xk):
+        def _select_row_index(self, xk):
             return 0
 
     return _DummyStrategy
@@ -37,6 +37,14 @@ def test_undefined_abstract_method(eye23, ones2, DummyStrategy, NonStrategy):
     DummyStrategy(eye23, ones2)
 
 
+@pytest.mark.timeout(1)
+def test_inconsistent_system_terminates(eye23, ones2, DummyStrategy, NonStrategy):
+    """Make sure inconsistent systems do not run forever."""
+    A = np.array([[1], [2]])
+    b = np.array([1, 1])
+    DummyStrategy.solve(A, b)
+
+
 def test_single_row_matrix(DummyStrategy, allclose):
     A = np.array([[0, 0, 1, 1]])
     b = np.array([1])
diff --git a/tests/test_cyclic.py b/tests/test_cyclic.py
index 591ef5d..63981f6 100644
--- a/tests/test_cyclic.py
+++ b/tests/test_cyclic.py
@@ -7,10 +7,10 @@ import kaczmarz
 def test_row_indexes(eye23, ones2):
     x = np.zeros(3)
     cyclic = kaczmarz.Cyclic(eye23, ones2)
-    assert 0 == cyclic.select_row_index(x)
-    assert 1 == cyclic.select_row_index(x)
-    assert 0 == cyclic.select_row_index(x)
-    assert 1 == cyclic.select_row_index(x)
+    assert 0 == cyclic._select_row_index(x)
+    assert 1 == cyclic._select_row_index(x)
+    assert 0 == cyclic._select_row_index(x)
+    assert 1 == cyclic._select_row_index(x)
 
     x0 = np.zeros(3)
     iterates = kaczmarz.Cyclic.iterates(eye23, ones2, x0)
@@ -28,7 +28,14 @@ def test_solve_identity(eye33, ones3):
 
 
 def test_solve_non_orthogonal_matrix():
-    A = np.array([[1, 0, 0], [2, 1, 0], [1, 2, 1], [0, 1, 2], [0, 0, 1]])
+    Alol = [
+        [1, 0, 0],
+        [2, 1, 0],
+        [1, 2, 1],
+        [0, 1, 2],
+        [0, 0, 1],
+    ]
+    A = np.array(Alol)
     row_norms = np.linalg.norm(A, axis=1)
     x_exact = np.ones(3)
     b = A @ x_exact