Working cupy methods (power method and eighs)

Author: lorenzotomada

.github/workflows/docs.yml

@@ -5,11 +5,13 @@ on:
     branches:
       - main
       - final_project
+      - gpu_branch
 
   pull_request:
     branches:
       - main
       - final_project
+      - gpu_branch
 
 jobs:
   docs:

.github/workflows/formatter.yml

@@ -5,11 +5,13 @@ on:
     branches:
       - main
       - final_project
+      - gpu_branch
 
   pull_request:
     branches:
       - main
       - final_project
+      - gpu_branch
 
 jobs:
   linter:

.github/workflows/tests.yml

@@ -23,13 +23,19 @@ jobs:
     steps:
     - uses: actions/checkout@v3
 
+    - name: Install system dependencies
+      run: |
+        apt-get update && apt-get install -y python3-venv python3-pip
+
     - name: Install additional dependencies
       run: |
         python3 -m venv venv
         source venv/bin/activate
         python -m pip install .
+      shell: bash -e {0}
 
     - name: Test with pytest
       run: |
         source venv/bin/activate
-        python -m pytest
+        python -m pytest
+      shell: bash -e {0}

experiments/config.yaml

@@ -1,3 +1,4 @@
-dimension: 1000
-
+dim: 1000
 density: 0.1
+tol: 0.001
+max_iter: 100

scripts/run.py

@@ -25,8 +25,9 @@
 
 
 kwargs = read_config(filename)
-dim = kwargs["dimension"]
+dim = kwargs["dim"]
 density = kwargs["density"]
+# also tol and max iter
 
 
 matrix = sp.random(dim, dim, density=density, format="csr")

shell/submit.sbatch

@@ -29,8 +29,6 @@ export SLURM_NTASKS_PER_NODE=1  # due to Ulysses's bug
 module load cuda/12.1
 conda init
 source ~/.bashrc
-conda activate ~/miniconda3/envs/devtools_scicomp
-pip install cupy-cuda12x
-python -m pip freeze > requirements.txt
-# Run the script
+conda activate /scratch/ltomada/miniconda3/envs/devtools_scicomp
+python -m pytest -v
 #python scripts/run.py fit --config=experiments/config.yaml

src/pyclassify/__init__.py

@@ -4,6 +4,7 @@
     "eigenvalues_cp",
     "power_method",
     "power_method_numba",
+    "power_method_cp",
 ]
 
 from .eigenvalues import (
@@ -12,4 +13,5 @@
     eigenvalues_cp,
     power_method,
     power_method_numba,
-)
+    power_method_cp,
+)

src/pyclassify/eigenvalues.py

@@ -3,11 +3,17 @@
 import scipy.linalg as spla
 from line_profiler import profile
 from numpy.linalg import eig, eigh
-from pyclassify.utils import check_square_matrix, check_symm_square, power_method_numba_helper
+from pyclassify.utils import (
+    check_square_matrix,
+    check_symm_square,
+    power_method_numba_helper,
+)
 import cupy as cp
+
 cp.cuda.Device(0)
 import cupy.linalg as cpla
 import cupyx.scipy.sparse as cpsp
+from cupyx.scipy.sparse.linalg._eigen import eigsh as cp_eigsh
 
 
 @profile
@@ -60,9 +66,9 @@ def eigenvalues_sp(A, symmetric=True):
     """
     check_square_matrix(A)
     eigenvalues, _ = (
-        spla.eigsh(A, k=A.shape[0] - 1)
+        sp.linalg.eigsh(A, k=A.shape[0] - 1)
         if symmetric
-        else spla.eigs(A, k=A.shape[0] - 1)
+        else sp.linalg.eigs(A, k=A.shape[0] - 1)
     )
     return eigenvalues
 
@@ -73,8 +79,16 @@ def eigenvalues_cp(A):
     Compute the eigenvalues of a sparse matrix using CuPy's `eigsh` function.
 
     This function checks if the input matrix is square and symmetric, then computes its eigenvalues using
-    CupY's sparse linear algebra solvers. For symmetric matrices, it uses `eigsh` for
-    more efficient computation.
+    CuPy's sparse linear algebra solvers. It uses `eigsh` for more efficient computation.
+    IMPORTANT: it important to underline a couple of things regarding this function:
+    - installing using the command
+    .. code-block:: shell
+            python -m pip install cupy-cuda12x
+      does not allow, for some reason, to import cupyx.scipy.sparse.linalg, and it necessary to import the
+      function manually form the source code. The problem can be observed using python3.10, while for
+      python3.12 things seem to work.
+    - the eighs function in this case does not allow to compute *all* the eigenvalues, but only a number
+    $m<n$, so here just a reduced portion is computed (starting form the ones which are greater in magnitude).
 
     Args:
         A (cpsp.spmatrix): A square sparse matrix whose eigenvalues are to be computed.
@@ -87,7 +101,8 @@ def eigenvalues_cp(A):
         ValueError: If number of rows != number of columns.
     """
     check_symm_square(A)
-    eigenvalues, _ = cpla.eigsh(A, k=A.shape[0] - 1)
+    k = 5 if A.shape[0] > 5 else A.shape[0] - 2
+    eigenvalues, _ = cp_eigsh(A, k=k)
     return eigenvalues
 
 
@@ -239,7 +254,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):
@@ -264,31 +278,45 @@ 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()
+    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

src/pyclassify/utils.py

@@ -13,9 +13,9 @@
 
 def check_square_matrix(A):
     """
-    Checks if the input matrix is a square matrix of type NumPy ndarray or SciPy sparse matrix.
+    Checks if the input matrix is a square matrix of type NumPy ndarray or SciPy or CuPy sparse matrix.
     This is done to ensure that the input matrix `A` is both:
-    1. Of type `np.ndarray` (NumPy array) or `scipy.sparse.spmatrix` (SciPy sparse matrix) or 'cupyx.scipy.sparse.spmatrix'.
+    1. Of type `np.ndarray` (NumPy array) or `scipy.sparse.spmatrix` (SciPy sparse matrix) or 'cupyx.scipy.sparse.spmatrix' (CuPy sparse matrix).
     2. A square matrix.
 
     Args:
@@ -26,7 +26,9 @@ def check_square_matrix(A):
         ValueError: If number of rows != number of columns.
     """
     if not isinstance(A, (np.ndarray, sp.spmatrix, cpsp.spmatrix)):
-        raise TypeError("Input matrix must be a NumPy array or a SciPy/CuPy sparse matrix!")
+        raise TypeError(
+            "Input matrix must be a NumPy array or a SciPy/CuPy sparse matrix!"
+        )
     if A.shape[0] != A.shape[1]:
         raise ValueError("Matrix must be square!")
 
@@ -72,7 +74,7 @@ def check_symm_square(A):
     check_square_matrix(A)
     if isinstance(A, sp.spmatrix) and not np.allclose(A.toarray(), A.toarray().T):
         raise ValueError("Matrix must be symmetric!")
-    elif isinstance(A, cpsp.spmatrix) and not cp.allclose(A.get(), A.get().T):
+    elif isinstance(A, cpsp.spmatrix) and not cp.allclose(A.toarray(), A.toarray().T):
         raise ValueError("Matrix must be symmetric!")
 
 
@@ -134,4 +136,4 @@ def read_config(file: str) -> dict:
     filepath = os.path.abspath(f"{file}.yaml")
     with open(filepath, "r") as stream:
         kwargs = yaml.safe_load(stream)
-    return kwargs
+    return kwargs

test/test_.py

@@ -1,9 +1,5 @@
 import numpy as np
 import scipy.sparse as sp
-import cupy as cp
-import os
-os.environ["CUPY_ACCELERATORS"] = "cub"
-cp.use_default_memory_pool = False
 import cupyx.scipy.sparse as cpsp
 import pytest
 from pyclassify import (
@@ -12,6 +8,7 @@
     eigenvalues_cp,
     power_method,
     power_method_numba,
+    power_method_cp,
 )
 from pyclassify.utils import check_square_matrix, make_symmetric, check_symm_square
 
@@ -26,7 +23,7 @@
 @pytest.mark.parametrize("size", sizes)
 @pytest.mark.parametrize("density", densities)
 def test_checks_on_A(size, density):
-    ugly_nonsquare_matrix = np.random.rand(size, size+1)
+    ugly_nonsquare_matrix = np.random.rand(size, size + 1)
 
     with pytest.raises(ValueError):
         check_square_matrix(ugly_nonsquare_matrix)
@@ -83,8 +80,20 @@ def test_implementations_power_method(size, density):
 
     biggest_eigenvalue_pm = power_method(matrix)
     biggest_eigenvalue_pm_numba = power_method_numba(matrix.toarray())
-
-    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_cp, biggest_eigenvalue_sp, rtol=1e-4)
-    assert np.isclose(biggest_eigenvalue_pm_numba, biggest_eigenvalue_sp, rtol=1e-4)
+    biggest_eigenvalue_pm_cp = power_method_cp(cp_matrix)
+
+    assert np.isclose(
+        biggest_eigenvalue_np, biggest_eigenvalue_sp, rtol=1e-4
+    )  # ensure numpy and scipy implementations are consistent
+    assert np.isclose(
+        biggest_eigenvalue_cp, biggest_eigenvalue_sp, rtol=1e-4
+    )  # ensure cupy and scipy implementations are consistent
+    assert np.isclose(
+        biggest_eigenvalue_pm, biggest_eigenvalue_sp, rtol=1e-4
+    )  # ensure power method and scipy implementation are consistent
+    assert np.isclose(
+        biggest_eigenvalue_pm_numba, biggest_eigenvalue_sp, rtol=1e-4
+    )  # ensure numba power method and scipy implementation are consistent
+    assert np.isclose(
+        biggest_eigenvalue_pm_cp, biggest_eigenvalue_sp, rtol=1e-4
+    )  # ensure cupy power method and scipy implementation are consistent