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,6 +23,10 @@ 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

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

@@ -8,6 +8,7 @@
 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,21 +61,29 @@ 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
 
 
 @profile
 def eigenvalues_cp(A):
-    """
+    r"""
     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 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 +96,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
 
 

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:
@@ -72,7 +72,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!")
 
 

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
 
@@ -83,8 +80,10 @@ def test_implementations_power_method(size, density):
 
     biggest_eigenvalue_pm = power_method(matrix)
     biggest_eigenvalue_pm_numba = power_method_numba(matrix.toarray())
+    biggest_eigenvalue_pm_cp = power_method_cp(cp_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)
-    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)
+    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