Keeping up to date wrt final_project branch

Author: lorenzotomada

docs/Documentation.ipynb

@@ -0,0 +1,70 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Documentation \n",
+    "**Disclaimer**: *this document is not meant to be neither a formal nor an exhaustive description of the iterative method to solve the eigenvalue problem. The aim of this text is only to provide the necessary and minimal information needed for the reader to understand the code implementation. We will use these pages to document, explain and justify the choice made from a numerical and scientific computing standpoint*.\n",
+    "\n",
+    "## Problem statement\n",
+    "Given a matrix $\\boldsymbol{A} \\in \\mathbb{C}^{n, n}$, with $n \\in \\mathbb{N}$ the matrix dimension, the eigenvalue problem can be formulated as finding the eigenpair $\\{(\\lambda_i, \\boldsymbol{v}_i)\\}_{i=1} ^ n$, with $\\lambda_i \\in \\mathbb{C}$, $\\boldsymbol{v}_i \\in \\mathbb{C}^{n ,1}$ and $\\boldsymbol{v}_i \\ne \\boldsymbol{0}$ such that \n",
+    "$$\\boldsymbol{A} \\boldsymbol{v}_i = \\lambda_i \\boldsymbol{v}_i, \\quad i=1, 2, ..., n $$\n",
+    "\n",
+    "$\\lambda_i$ are called eigenvalue, while $\\boldsymbol{v}_i$ are the corresponding eigenvectors. Theoretically, finding the matrix's eigenvalue is possible by imposing the following condition:\n",
+    "$$ \\det(A-\\lambda I)=0 $$\n",
+    "which lead to the well known characteristic polynomial of degree $n$. The eigenvalue are the roots of the characteristic polynomial. Although correct, this apporach is not viable, due to the difficulties in both expressing the characteristic polynomial and find its roots, when the Matrix gets very large.\n",
+    "\n",
+    "Numerical methods takle the eigenproblem using a different strategy, most of the time being a iterative methods. Rather than aspire to get the exact solution, they seek after an approximation of the solution, which hopefully, under a set of conditions, converges to the exact solution.\n",
+    "Among all the method devolped to solve the eigenproblem, the power method, along with its variants suct that the inverse power method and power method with shift, and the QR are the widest spread and the most used.\n",
+    "\n",
+    "## Power method\n",
+    "Let $\\boldsymbol{A}$ being a matric and $\\lambda _i$ the eigenvalue already sorted accordingly to their module. If\n",
+    "$$ |\\lambda_1| > |\\lambda _i | \\quad i=2, 3, ..., n$$\n",
+    "then the power method allow to recover the eigenpair $(\\lambda_1, \\boldsymbol{v}_1)$ by applying iteratively, the following steps\n",
+    "\n",
+    "\n",
+    "$$\n",
+    "\\begin{array}{l}\n",
+    "\\textbf{Power Method Algorithm} \\\\\n",
+    "\\textbf{Input:} A \\in \\mathbb{C}^{n \\times n}, \\text{initial vector } x_0, \\text{ tolerance } \\text{tol}, \\text{ max iterations } \\text{max\\_iter} \\\\\n",
+    "\\textbf{Output:} \\lambda \\text{ (dominant eigenvalue approximation), } x \\text{ (eigenvector approximation)} \\\\\n",
+    "\n",
+    "1.\\ \\text{Initialize } x = x_0 \\text{ (random or chosen guess)} \\\\\n",
+    "2.\\ \\text{Normalize } x \\text{: } y^{(1)} = x/ \\| x\\| \\\\\n",
+    "3.\\ \\lambda_{\\text{old}} = 0 \\\\\n",
+    "\n",
+    "4.\\ \\textbf{for } k = 1 \\textbf{ to } \\text{max\\_iter} \\textbf{ do} \\\\\n",
+    "\\quad 5.\\ x^{(k+1)} = A \\cdot y^{(k)} \\\\\n",
+    "\\quad 6.\\  y^{(k+1)} = x^{(k+1)}/ \\| x^{(k+1)}\\|  \\\\\n",
+    "\\quad 7.\\  \\lambda^{(k+1)} =  (y^{(k+1)})^H \\boldsymbol{A} y^{(k+1)} \\\\\n",
+    "\\quad 8.\\ \\textbf{if } |\\lambda^{(k+1)} - \\lambda^{(k)}| < \\text{tol} \\textbf{ then} \\\\\n",
+    "\\quad \\quad \\text{Break (convergence reached)} \\\\\n",
+    "\n",
+    "9.\\ \\textbf{Return } \\lambda, x\n",
+    "\\end{array}\n",
+    "$$\n",
+    "\n",
+    "By exploiting Matrix properties, it is also possible to find the eigenvalue with the smallest norm, and its associated eigenvector, or either the the eigenvalue closest to a given number.\n",
+    "\n",
+    "## QR algorithm\n",
+    "Being the problem we are interested in to solve incolve a symmetric Matrix, we will start introduce the QR alforithm in a general way and specialize it for the case of symmetric matrices as soon as we have the chance.\n",
+    "\n",
+    "The QR algorithm can be divide in to two stages: during the first stage, the matrix is reduced, by means of similar trasformation, into a Hessemberg matrix, or, for Hermetian matrix, into a tridiagonal matrix. Although the QR algorithm can be applied directly to a full matrix, it converges faster if it is applied to Hessemberg Matrix or triagular matrices, with also a smaller computational cost.\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.12.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

docs/Makefile

@@ -17,4 +17,5 @@ help:
 # Catch-all target: route all unknown targets to Sphinx using the new
 # "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
 %: Makefile
+	export LC_ALL=C
 	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)

requirements.txt

@@ -1,4 +1,5 @@
 alabaster==1.0.0
+asttokens==3.0.0
 babel==2.17.0
 black==25.1.0
 certifi==2025.1.31
@@ -7,27 +8,38 @@ click==8.1.8
 contourpy==1.3.1
 cupy-cuda12x==13.4.0
 cycler==0.12.1
+decorator==5.2.1
 docutils==0.21.2
 exceptiongroup==1.2.2
+executing==2.2.0
 fastrlock==0.8.3
 fonttools==4.56.0
 idna==3.10
 imagesize==1.4.1
 iniconfig==2.0.0
+ipython==9.0.1
+ipython_pygments_lexers==1.1.1
+jedi==0.19.2
 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
+matplotlib-inline==0.1.7
 mypy-extensions==1.0.0
 numba==0.61.0
 numpy==2.1.3
 packaging==24.2
+parso==0.8.4
 pathspec==0.12.1
+pexpect==4.9.0
 pillow==11.1.0
 platformdirs==4.3.6
 pluggy==1.5.0
+prompt_toolkit==3.0.50
+ptyprocess==0.7.0
+pure_eval==0.2.3
 Pygments==2.19.1
 pyparsing==3.2.1
 pytest==8.3.4
@@ -47,7 +59,11 @@ sphinxcontrib-jquery==4.1
 sphinxcontrib-jsmath==1.0.1
 sphinxcontrib-qthelp==2.0.0
 sphinxcontrib-serializinghtml==2.0.0
+stack-data==0.6.3
+tokenize_rt==6.1.0
 tomli==2.2.1
+traitlets==5.14.3
 typing_extensions==4.12.2
 urllib3==2.3.0
+wcwidth==0.2.13
 wheel==0.45.1

shell/submit.sh

@@ -0,0 +1,6 @@
+#!/bin/bash
+
+module load cuda/12.1
+module use cuda/12.1
+
+pytest -v > output_pytest.txt

src/pyclassify/__init__.py

@@ -5,6 +5,8 @@
     "power_method",
     "power_method_numba",
     "power_method_cp",
+    "Lanczos_PRO",
+    "QR_method",
 ]
 
 from .eigenvalues import (
@@ -14,4 +16,6 @@
     power_method,
     power_method_numba,
     power_method_cp,
+    Lanczos_PRO,
+    QR_method,
 )

src/pyclassify/eigenvalues.py

@@ -137,7 +137,7 @@ def power_method(A, max_iter=500, tol=1e-4, x=None):
 
 
 @profile
-def power_method_numba(A):
+def power_method_numba(A, max_iter=500, tol=1e-4, x=None):
     """
     Compute the dominant eigenvalue of a matrix using the power method, with Numba optimization.
 
@@ -154,7 +154,7 @@ def power_method_numba(A):
     Returns:
         float: The approximated dominant eigenvalue of the matrix `A`.
     """
-    return power_method_numba_helper(A)
+    return power_method_numba_helper(A, max_iter, tol, x)
 
 
 @profile
@@ -219,6 +219,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
@@ -271,14 +280,12 @@ def QR_method(A_copy, tol=1e-10, max_iter=100):
     """
 
     A = A_copy.copy()
-    T = A.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([])
-        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)

test/test_.py

@@ -1,6 +1,8 @@
 import numpy as np
 import scipy.sparse as sp
 import cupyx.scipy.sparse as cpsp
+import cupy as cp
+import scipy
 import pytest
 from pyclassify import (
     eigenvalues_np,
@@ -9,14 +11,19 @@
     power_method,
     power_method_numba,
     power_method_cp,
+    Lanczos_PRO,
+    QR_method,
 )
 from pyclassify.utils import check_square_matrix, make_symmetric, check_symm_square
 
 
-np.random.seed(104)
+@pytest.fixture(autouse=True)
+def set_random_seed():
+    seed = 1422
+    np.random.seed(seed)
 
 
-sizes = [10, 100, 1000]
+sizes = [10, 100]
 densities = [0.1, 0.3]
 
 
@@ -34,8 +41,7 @@ def test_checks_on_A(size, density):
     symmetric_matrix = make_symmetric(matrix)
     check_square_matrix(symmetric_matrix)
 
-    # cp_symm_matrix = cpsp.csr_matrix(symmetric_matrix)
-    # check_square_matrix(cp_symm_matrix)
+    # regarding the CuPy implementation: see below!
 
     not_so_symmetric_matrix = np.random.rand(5, 5)
     if not_so_symmetric_matrix[1, 2] == not_so_symmetric_matrix[2, 1]:
@@ -64,36 +70,112 @@ def test_make_symmetric(size, density):
 def test_implementations_power_method(size, density):
     matrix = sp.random(size, size, density=density, format="csr")
     matrix = make_symmetric(matrix)
-    # cp_matrix = cpsp.csr_matrix(matrix)
 
     eigs_np = eigenvalues_np(matrix.toarray(), symmetric=True)
     eigs_sp = eigenvalues_sp(matrix, symmetric=True)
-    # eigs_cp = eigenvalues_cp(cp_matrix)
 
     index_np = np.argmax(np.abs(eigs_np))
     index_sp = np.argmax(np.abs(eigs_sp))
-    # index_cp = np.argmax(np.abs(eigs_cp))
 
     biggest_eigenvalue_np = eigs_np[index_np]
     biggest_eigenvalue_sp = eigs_sp[index_sp]
-    # biggest_eigenvalue_cp = eigs_cp[index_cp]
 
     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
     )  # 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
+
+
+@pytest.mark.parametrize("size", sizes)
+@pytest.mark.parametrize("density", densities)
+def test_cupy(size, density):
+    try:
+        if not cp.cuda.is_available():
+            pytest.skip("Skipping test because CUDA is not available")
+    except cp.cuda.runtime.CUDARuntimeError as e:
+        pytest.skip(f"Skipping test due to CUDA driver issues: {str(e)}")
+
+    matrix = sp.random(size, size, density=density, format="csr")
+    symmetric_matrix = make_symmetric(matrix)
+    cp_symm_matrix = cpsp.csr_matrix(symmetric_matrix)
+    check_square_matrix(cp_symm_matrix)
+
+    eigs_cp = eigenvalues_cp(cp_symm_matrix)
+
+    index_cp = np.argmax(np.abs(eigs_cp))
+    biggest_eigenvalue_cp = eigs_cp[index_cp]
+    biggest_eigenvalue_pm_cp = power_method_cp(cp_symm_matrix, max_iter=1000, tol=1e-5)
+
+    assert np.isclose(
+        biggest_eigenvalue_cp, biggest_eigenvalue_pm_cp, rtol=1e-4
+    )  # ensure cupy native and cupy power method implementations are consistent
+
+
+@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)