Add basic CI w/ GH actions

.github/workflows/test.yml

@@ -0,0 +1,21 @@
+name: Run Tests
+
+on:
+  push:
+    branches: [main]
+  pull_request:
+    branches: [main]
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: checkout code
+        uses: actions/checkout@v4
+      - name: Install UV
+        uses: astral-sh/setup-uv@v5
+      - name: Install the project
+        run: uv sync --group test
+      - name: Run tests
+        run: uv run pytest tests

pyproject.toml

@@ -8,6 +8,7 @@ authors = [
 ]
 dependencies = [
     "numpy>=2.0",
+    "scipy>=1.15.3",
 ]
 dynamic = ["version"]
 
@@ -21,8 +22,11 @@ path = "src/arnoldi/_version.py"
 [dependency-groups]
 dev = [
     "notebook>=7.4.2",
-    "scipy>=1.15.3",
 ]
 test = [
     "pytest>=8.3.5",
+    "pytest-cov>=6.2.1",
 ]
+
+[tool.pytest.ini_options]
+addopts = "--cov=arnoldi --cov-report=term-missing --cov-branch"

src/arnoldi/decomposition.py

@@ -1,4 +1,5 @@
 import numpy as np
+import scipy.linalg as slin
 import scipy.sparse as sp
 
 
@@ -13,16 +14,24 @@ def __init__(self, n, k, dtype=np.complex128):
         self.q = np.zeros((n, k+1), dtype=dtype)
         self.h = np.zeros((k+1, k), dtype=dtype)
 
-        self.dtype = dtype
+    @property
+    def _dtype(self):
+        return self.h.dtype
+
+    @property
+    def _atol(self):
+        # XXX: logic of sqrt copied from Julia's ArnoldiMethod.jl package
+        return np.sqrt(np.finfo(self._dtype).eps)
 
     def initialize(self, init_vector=None):
         if init_vector is None:
-            init_vector = np.random.randn(self.n).astype(self.dtype)
+            init_vector = np.random.randn(self.n).astype(self._dtype)
             init_vector /= np.linalg.norm(init_vector)
 
         self.q[:, 0] = init_vector
 
-    def iterate(self, a, tol=1e-12):
+    def iterate(self, a, tol=None):
+        tol = tol or self._atol
         for j in range(self.k):
             v = self.q[:, j+1]
             v[:] = a @ self.q[:, j]
@@ -34,7 +43,7 @@ def iterate(self, a, tol=1e-12):
 
             self.h[j + 1, j] = np.linalg.norm(v)
 
-            if self.h[j + 1, j] < tol:
+            if self.h[j + 1, j] < self._atol:
                 return j
             v /= self.h[j + 1, j]
 

tests/test_decomposition.py

@@ -4,18 +4,22 @@
 from arnoldi import Arnoldi
 
 
+ATOL = 1e-8
+RTOL = 1e-4
+
+
 def largest_eigvals(m, k):
     """ Compute top k eigenvalues of m, when sorted by amplitude.
     """
     r_eigvals, _ = np.linalg.eig(m)
     return np.array(sorted(r_eigvals, key=np.abs, reverse=True)[:k])
 
 
-RTOL = 1e-5
-ATOL= 1e-8
+class TestArnoldiExpansion:
+    def test_invariant_simple(self):
+        ## Test the invariant A * Q ~ Q * H, with H Hessenberg matrix and Q
+        ## orthonormal
 
-class TestArnoldiSimple:
-    def test_simple(self):
         # Given
         n = 10
         k = 6
@@ -41,27 +45,27 @@ def test_simple(self):
             a @ q[:, :n_iter], q @ h, rtol=RTOL, atol=ATOL
         )
 
-    def test_eigvals(self):
+    def test_eigvals_simple(self):
         # Given
         n = 20
         ## We use large krylov space since basic arnoldi does not converge
         ## quickly. Once we implement restarts, it should converge much faster
         k = n - 1
 
         a = sp.random(n, n, density=5 / n, dtype=np.complex128)
+        # We add ones on the diag to have a well conditioned matrix and eigen
+        # values not too far from 1
         a += sp.diags_array(np.ones(n))
 
         # When
         arnoldi = Arnoldi(n, k)
         arnoldi.initialize()
         n_iter = arnoldi.iterate(a)
 
-        q, h = arnoldi.q[:, :n_iter+1], arnoldi.h[:n_iter+1, :n_iter+1]
+        q, h = arnoldi.q[:, :n_iter+1], arnoldi.h[:n_iter, :n_iter]
 
         r_eigvals = largest_eigvals(a.toarray(), 3)
-        ritz_values = largest_eigvals(h[:-1, :], 3)
+        ritz_values = largest_eigvals(h, 3)
 
         # Then
-        np.testing.assert_allclose(
-            r_eigvals, ritz_values, rtol=1000 * RTOL, atol=ATOL
-        )
+        np.testing.assert_allclose(r_eigvals, ritz_values, rtol=RTOL, atol=ATOL)