Adding power method and tests

Author: lorenzotomada

requirements.txt

@@ -15,10 +15,12 @@ iniconfig==2.0.0
 Jinja2==3.1.5
 kiwisolver==1.4.8
 line_profiler==4.2.0
+llvmlite==0.44.0
 MarkupSafe==3.0.2
 matplotlib==3.10.0
 mypy-extensions==1.0.0
-numpy==2.2.3
+numba==0.61.0
+numpy==2.1.3
 packaging==24.2
 pathspec==0.12.1
 pillow==11.1.0

src/pyclassify/utils.py

@@ -0,0 +1,55 @@
+import numpy as np
+import scipy.sparse as sp
+from numpy.linalg import eig, eigh
+
+
+
+def make_symmetric(A):
+    if not sp.isspmatrix(A):
+        raise TypeError("Input must be a scipy sparse matrix")
+    
+    A_sym = (A + A.T) / 2
+    return A_sym
+
+
+
+def eigenvalues_np(A, symmetric=True):
+    if A.shape[0] != A.shape[1]:
+        raise ValueError("Matrix must be square!")
+    if not isinstance(A, (np.ndarray, sp.spmatrix)):
+        raise TypeError("Input matrix must be a NumPy array or a SciPy sparse matrix!")
+    
+    eigenvalues, _ = eigh(A) if symmetric else eig(A)
+    return eigenvalues
+
+
+
+def eigenvalues_sp(A, symmetric=True):
+    if A.shape[0] != A.shape[1]:
+        raise ValueError("Matrix must be square!")
+    if not isinstance(A, (np.ndarray, sp.spmatrix)):
+        raise TypeError("Input matrix must be a NumPy array or a SciPy sparse matrix!")
+
+    eigenvalues, _ = sp.linalg.eigsh(A, k=A.shape[0] - 1) if symmetric else sp.linalg.eigs(A, k=A.shape[0] - 1)
+    return eigenvalues
+
+
+
+def power_method(A, max_iter=1000, tol=1e-4):
+    assert A.shape[0] == A.shape[1], "Matrix must be square"
+
+    x = np.random.rand(A.shape[0])
+    x /= np.linalg.norm(x)
+    x_old = x
+
+    iteration = 0
+    update_norm = tol + 1
+
+    while iteration < max_iter and update_norm > tol:
+        x = A @ x
+        x /= np.linalg.norm(x)
+        update_norm = np.linalg.norm(x - x_old)/np.linalg.norm(x_old)
+        x_old = x.copy()
+        iteration += 1
+
+    return x @ A @ x

test/test_.py

@@ -0,0 +1,53 @@
+import numpy as np
+import scipy.sparse as sp
+import pytest
+from pyclassify.utils import eigenvalues_np, eigenvalues_sp, power_method, make_symmetric
+
+
+np.random.seed(42)
+sizes = [10, 100, 1000]
+densities = [0.1, 0.5]
+
+
+
+@pytest.mark.parametrize("size", sizes)
+@pytest.mark.parametrize("density", densities)
+def test_make_symmetric(size, density):
+    matrix = sp.random(size, size, density=density, format="csr")
+    symmetric_matrix = make_symmetric(matrix)
+    
+    assert np.allclose(symmetric_matrix.toarray(), symmetric_matrix.toarray().T)
+    assert symmetric_matrix.shape == matrix.shape
+
+    with pytest.raises(TypeError):
+        _ = make_symmetric("banana")
+
+
+
+@pytest.mark.parametrize("size", sizes)
+@pytest.mark.parametrize("density", densities)
+def test_power_method(size, density):
+    matrix = sp.random(size, size, density=density, format="csr")
+    matrix = make_symmetric(matrix)
+
+    eigs_np = eigenvalues_np(matrix.toarray(), symmetric=True)
+    eigs_sp = eigenvalues_sp(matrix, symmetric=True)
+
+    index_np = np.argmax(np.abs(eigs_np))
+    index_sp = np.argmax(np.abs(eigs_sp))
+
+    biggest_eigenvalue_np = eigs_np[index_np]
+    biggest_eigenvalue_sp = eigs_sp[index_sp]
+
+    biggest_eigenvalue_pm = power_method(matrix)
+
+    assert np.isclose(biggest_eigenvalue_np, biggest_eigenvalue_sp, rtol=1e-4)
+    assert np.isclose(biggest_eigenvalue_pm, biggest_eigenvalue_sp, rtol=1e-4)
+
+    ugly_nonsquare_matrix = np.random.rand(5, 3)
+
+    with pytest.raises(ValueError):
+        _ = eigenvalues_np(eigenvalues_np(ugly_nonsquare_matrix, symmetric=True))
+
+    with pytest.raises(ValueError):
+        _ = eigenvalues_sp(eigenvalues_sp(ugly_nonsquare_matrix, symmetric=True))