Added tests for Lanczos and QR

Author: Gaspare99

src/pyclassify/__init__.py

@@ -3,11 +3,15 @@
     "eigenvalues_sp",
     "power_method",
     "power_method_numba",
+    "Lanczos_PRO",
+    "QR_method",
 ]
 
 from .eigenvalues import (
     eigenvalues_np,
     eigenvalues_sp,
     power_method,
     power_method_numba,
+    Lanczos_PRO,
+    QR_method,
 )

src/pyclassify/eigenvalues.py

@@ -145,6 +145,15 @@ def Lanczos_PRO(A, q, m=None, toll=np.sqrt(np.finfo(float).eps)):
     Raises:
         ValueError: If the input matrix A is not square or if m is greater than the size of A.
     """
+    if m == None:
+        m = A.shape[0]
+
+    if A.shape[0] != A.shape[1]:
+        raise ValueError("Input matrix A must be square.")
+
+    if A.shape[0] != q.shape[0]:
+        raise ValueError("Input vector q must have the same size as the matrix A.")
+
     q = q / np.linalg.norm(q)
     Q = np.array([q])
     r = A @ q
@@ -172,7 +181,6 @@ def Lanczos_PRO(A, q, m=None, toll=np.sqrt(np.finfo(float).eps)):
         if np.abs(beta[j]) < 1e-15:
             return Q, alpha, beta[:-1]
     return Q, alpha, beta[:-1]
-   
 
 
 def QR_method(A_copy, tol=1e-10, max_iter=100):
@@ -197,32 +205,43 @@ def QR_method(A_copy, tol=1e-10, max_iter=100):
         ValueError: If the input matrix A_copy is not square.
     """
 
-    A=A_copy.copy()
-    T=A.copy()  
-    A=np.array(A)
-    iter=0
-    
-    while np.linalg.norm(np.diag(A, -1), np.inf) > tol and iter<max_iter:
-        Matrix_trigonometry=np.array([]) 
-        QQ, RR=np.linalg.qr(T)
-        T=RR@QQ
-        for i in range(A.shape[0]-1):
-            c=A[i, i]/np.sqrt(A[i, i]**2+A[i+1, i]**2)
-            s=-A[i+1, i]/np.sqrt(A[i, i]**2+A[i+1, i]**2)
-            Matrix_trigonometry=np.vstack((Matrix_trigonometry, [c, s])) if Matrix_trigonometry.size else np.array([[c, s]])
-
-            R=np.array([[c, -s], [s, c]])
-            A[i:i+2, i:i+3]=R@A[i:i+2, i:i+3]
-            A[i+1, i]=0
-        Q=np.eye(A.shape[0])
-        i=0
-        Q[0:2, 0:2]=np.array([[Matrix_trigonometry[i, 0], Matrix_trigonometry[i, 1]], [-Matrix_trigonometry[i, 1], Matrix_trigonometry[i, 0]]])
-        for i in range(1, A.shape[0]-1):
-            R=np.eye(A.shape[0])
-            R[i: i+2, i:i+2]=np.array([[Matrix_trigonometry[i, 0], Matrix_trigonometry[i, 1]], [-Matrix_trigonometry[i, 1], Matrix_trigonometry[i, 0]]])
-            Q = Q@R
-        A=A@Q
-        iter+=1
-
-    
-    return np.diag(A), Q
+    A = A_copy.copy()
+    A = np.array(A)
+    iter = 0
+    Q = np.array([])
+
+    while np.linalg.norm(np.diag(A, -1), np.inf) > tol and iter < max_iter:
+        Matrix_trigonometry = np.array([])
+        for i in range(A.shape[0] - 1):
+            c = A[i, i] / np.sqrt(A[i, i] ** 2 + A[i + 1, i] ** 2)
+            s = -A[i + 1, i] / np.sqrt(A[i, i] ** 2 + A[i + 1, i] ** 2)
+            Matrix_trigonometry = (
+                np.vstack((Matrix_trigonometry, [c, s]))
+                if Matrix_trigonometry.size
+                else np.array([[c, s]])
+            )
+
+            R = np.array([[c, -s], [s, c]])
+            A[i : i + 2, i : i + 3] = R @ A[i : i + 2, i : i + 3]
+            A[i + 1, i] = 0
+        Q = np.eye(A.shape[0])
+        i = 0
+        Q[0:2, 0:2] = np.array(
+            [
+                [Matrix_trigonometry[i, 0], Matrix_trigonometry[i, 1]],
+                [-Matrix_trigonometry[i, 1], Matrix_trigonometry[i, 0]],
+            ]
+        )
+        for i in range(1, A.shape[0] - 1):
+            R = np.eye(A.shape[0])
+            R[i : i + 2, i : i + 2] = np.array(
+                [
+                    [Matrix_trigonometry[i, 0], Matrix_trigonometry[i, 1]],
+                    [-Matrix_trigonometry[i, 1], Matrix_trigonometry[i, 0]],
+                ]
+            )
+            Q = Q @ R
+        A = A @ Q
+        iter += 1
+
+    return np.diag(A), Q

test/test_.py

@@ -1,11 +1,14 @@
 import numpy as np
 import scipy.sparse as sp
+import scipy
 import pytest
 from pyclassify import (
     eigenvalues_np,
     eigenvalues_sp,
     power_method,
     power_method_numba,
+    Lanczos_PRO,
+    QR_method,
 )
 from pyclassify.utils import check_A_square_matrix, make_symmetric
 
@@ -60,3 +63,65 @@ def test_implementations_power_method(size, density):
     assert np.isclose(biggest_eigenvalue_np, biggest_eigenvalue_sp, rtol=1e-4)
     assert np.isclose(biggest_eigenvalue_pm, biggest_eigenvalue_sp, rtol=1e-4)
     assert np.isclose(biggest_eigenvalue_pm_numba, biggest_eigenvalue_sp, rtol=1e-4)
+
+
+@pytest.mark.parametrize("size", sizes)
+def test_Lanczos(size):
+    matrix = sp.random(size, size, density=0.3, format="csr")
+    matrix = make_symmetric(matrix)
+
+    random_vector = np.random.rand(size)
+
+    _, alpha, beta = Lanczos_PRO(matrix, random_vector)
+
+    T = np.diag(alpha) + np.diag(beta, 1) + np.diag(beta, -1)
+
+    EigVal_T = np.linalg.eig(T)[0]
+    EigVect_T = np.linalg.eig(T)[1]
+
+    EigVal_A = np.linalg.eig(matrix.toarray())[0]
+    EigVect_A = np.linalg.eig(matrix.toarray())[1]
+
+    EigVal_A = np.sort(EigVal_A)
+    EigVal_T = np.sort(EigVal_T)
+
+    assert np.allclose(EigVal_T, EigVal_A, rtol=1e-6)
+
+    with pytest.raises(ValueError):
+        random_matrix = np.random.rand(size, 2 * size)
+        _ = Lanczos_PRO(random_matrix, random_vector)
+
+    with pytest.raises(ValueError):
+        random_matrix = np.random.rand(size, size)
+        _ = Lanczos_PRO(random_matrix, np.random.rand(2 * size))
+
+
+@pytest.mark.parametrize("size", [10, 100])
+def test_QR_method(size):
+    eig = np.arange(1, size + 1)
+    A = np.diag(eig)
+    U = scipy.stats.ortho_group.rvs(size)
+    print(U)
+
+    A = U @ A @ U.T
+    print(A)
+    A = make_symmetric(A)
+    eig = np.linalg.eig(A)
+    print(eig.eigenvalues)
+    index = np.argsort(eig.eigenvalues)
+    eig = eig.eigenvalues
+    eig_vec = np.linalg.eig(A).eigenvectors
+    eig_vec = eig_vec[index]
+    eig = eig[index]
+    eig_vec = eig_vec / np.linalg.norm(eig_vec, axis=0)
+
+    random_vector = np.random.rand(size)
+    _, alpha, beta = Lanczos_PRO(A, random_vector)
+
+    T = np.diag(alpha) + np.diag(beta, 1) + np.diag(beta, -1)
+
+    eig_valQR, eig_vecQR = QR_method(T, max_iter=100 * size)
+    index = np.argsort(eig_valQR)
+    eig_valQR = eig_valQR[index]
+
+    assert np.allclose(eig, eig_valQR, rtol=1e-8)