final_project documentation

Welcome to the documentation for the final project of the course in Development Tools for Scientific Computing. Our goal is to compute an efficient solver for eigenvalue problems.

For theoretical details, please check the file docs/Documentation.ipynb To install the package, please read the README.md file.

Table of Contents

Contents:

• PyClassify Module Documentation
• Functions:
  • EigenSolver
  • Lanczos_PRO()
  • eigenvalues_np()
  • eigenvalues_sp()
  • power_method()
  • power_method_numba()
• Submodule: utils
  • check_square_matrix()
  • check_symm_square()
  • make_symmetric()
  • max_iteration_warning()
  • poisson_2d_structure()
  • profile_numpy_eigvals()
  • profile_scipy_eigvals()
  • read_config()
• Submodule: helpers_secular
  • bisection()
  • compute_Psi()
  • compute_outer_zero()
  • find_root()
  • secular_solver_python()
• Submodule: parallel_tridiag_eigen
  • check_column_directions()
  • find_interval_extreme()
  • parallel_eigen()
  • parallel_tridiag_eigen()

Module Documentation

class pyclassify.EigenSolver(A: ndarray, max_iter=5000, tol=1e-08, tol_deflation=1e-12)

Bases: object

Class solving the eigenvalue problem for a given symmetric matrix.

Two different building blocks are present: Lanczos_PRO (used to transform the matrix to a tridiagonal one), and another function written in C++. The latter can be either Eigen_value_calculator (if the user is only interested in the computation of the eigenvalues) or QR_algorithm, if eigenvectors are needed as well.

We refer the interested reader to their implementation in C++ for further details.

Lanczos_PRO(A=None, q=None, m=None, tol=np.float64(1.4901161193847656e-08))

Perform the Lanczos algorithm for symmetric matrices.

This function computes an orthogonal matrix Q and tridiagonal matrix T such that .. math:: A approx Q T Q^T, where A is a symmetric matrix. The algorithm is useful for finding a few eigenvalues and eigenvectors of large symmetric matrices.

Parameters:

• A (np.ndarray) – A symmetric square matrix of size n x n.

• q (np.ndarray, optional) – Initial vector of size n. Default value is None (a random one is created).

• m (int, optional) – Number of eigenvalues to compute. Must be less than or equal to n. If None, defaults to the size of A.

• tol (float, optional) – Tolerance for orthogonality checks (default is sqrt(machine epsilon)).

Returns:

A tuple (Q, alpha, beta) where:

• Q (np.ndarray): Orthogonal matrix of size n x m.

• alpha (np.ndarray): Vector of size m containing the diagonal elements of the tridiagonal matrix.

• beta (np.ndarray): Vector of size m-1 containing the off-diagonal elements of the tridiagonal matrix.

Return type:

tuple

Raises:

• TypeError – If the input is not a NumPy array or SciPy/CuPy sparse matrix.

• ValueError – If number of rows != number of columns or the matrix is not symmetric or it m is greater than the size of A.

compute_eigenval(diag=None, off_diag=None)

Compute (only) the eigenvalues of a symmetric triangular matrix, passed as argument in the form of diagonal and off-diagonal.

This function relies on the Eigen_value_calculator function, written in C++. :param diag: Diagonal of the matrix. Default value is None. If no value is passed, the one used

is the one resulting from the Lanczos decomposition of the matrix passed to the constructor.

Parameters:

off_diag (np.ndarray, optional) – Off-iagonal of the matrix. Default value is None. If no value is passed, the one used is the one resulting from the Lanczos decomposition of the matrix passed to the constructor.

Returns:

an np.array containing the eigenvalues of the matrix.

Return type:

np.array

Raises:

ValueError – If there is a mismatch between the diagonal and off diagonal size.

eig(diag=None, off_diag=None)

Compute the eigenvalues and eigenvectors of a symmetric triangular matrix, passed as argument in the form of diagonal and off-diagonal.

This function relies on the QR_algorithm function, written in C++. :param diag: Diagonal of the matrix. Default value is None. If no value is passed, the one used

is the one resulting from the Lanczos decomposition of the matrix passed to the constructor.

Parameters:

off_diag (np.ndarray, optional) – Off-iagonal of the matrix. Default value is None. If no value is passed, the one used is the one resulting from the Lanczos decomposition of the matrix passed to the constructor.

Returns:

a tuple containing two arrays, the eigenvalues and the eigenvectors.

Return type:

tuple

Raises:

ValueError – If there is a mismatch between the diagonal and off diagonal size.

property initial_guess

pyclassify.Lanczos_PRO(A, q, m=None, tol=np.float64(1.4901161193847656e-08))

Perform the Lanczos algorithm for symmetric matrices.

This function computes an orthogonal matrix Q and tridiagonal matrix T such that .. math:: A approx Q T Q^T, where A is a symmetric matrix. The algorithm is useful for finding a few eigenvalues and eigenvectors of large symmetric matrices.

Parameters:

• A (np.ndarray) – A symmetric square matrix of size n x n.

• q (np.ndarray, optional) – Initial vector of size n. Default value is None (a random one is created).

• m (int, optional) – Number of eigenvalues to compute. Must be less than or equal to n. If None, defaults to the size of A.

• tol (float, optional) – Tolerance for orthogonality checks (default is sqrt(machine epsilon)).

Returns:

A tuple (Q, alpha, beta) where:

• Q (np.ndarray): Orthogonal matrix of size n x m.

• alpha (np.ndarray): Vector of size m containing the diagonal elements of the tridiagonal matrix.

• beta (np.ndarray): Vector of size m-1 containing the off-diagonal elements of the tridiagonal matrix.

Return type:

tuple

Raises:

• TypeError – If the input is not a NumPy array or SciPy/CuPy sparse matrix.

• ValueError – If number of rows != number of columns or the matrix is not symmetric or it m is greater than the size of A.

pyclassify.eigenvalues_np(A, symmetric=True)

Wrapper for the np eigenvalue solver. This function is only used in tests for better readability. Compute the eigenvalues of a square matrix using NumPy’s eig or eigh function.

This function checks if the input matrix is square (and is actually a matrix) using ‘check_square_matrix’, and then computes its eigenvalues.

If the matrix is symmetric, it uses eigh (which is more efficient for symmetric matrices). Otherwise, it uses eig.

Parameters:

• A (np.ndarray) – A square matrix whose eigenvalues are to be computed.

• symmetric (bool, optional) – If True, assumes the matrix is symmetric and uses eigh for faster computation. If False, uses eig for general matrices (default is True).

Returns:

An array containing the eigenvalues of the matrix A.

Return type:

np.ndarray

Raises:

• TypeError – If the input is not a NumPy array or a SciPy/CuPy sparse matrix.

• ValueError – If number of rows != number of columns.

pyclassify.eigenvalues_sp(A, symmetric=True)

Compute the eigenvalues of a sparse matrix using SciPy’s eigsh or eigs function.

This function checks if the input matrix is square, then computes its eigenvalues using SciPy’s sparse linear algebra solvers. For symmetric matrices, it uses eigsh for more efficient computation. For non-symmetric matrices, it uses eigs.

Parameters:

• A (sp.spmatrix) – A square sparse matrix whose eigenvalues are to be computed.

• symmetric (bool, optional) – If True, assumes the matrix is symmetric and uses eigsh. If False, uses eigs for general matrices (default is True).

Returns:

An array containing the eigenvalues of the sparse matrix A.

Return type:

np.ndarray

Raises:

• TypeError – If the input is not a NumPy array or a SciPy/CuPy sparse matrix.

• ValueError – If number of rows != number of columns.

pyclassify.power_method(A, max_iter=500, tol=1e-07, x=None)

Compute the dominant eigenvalue of a square matrix using the power method.

Parameters:

• A (np.ndarray or sp.spmatrix or cpsp.spmatrix) – A square matrix whose dominant eigenvalue is to be computed.

• max_iter (int, optional) – Maximum number of iterations to perform (default is 500).

• tol (float, optional) – Tolerance for convergence based on the relative change between iterations (default is 1e-4).

• x (np.ndarray, optional) – Initial guess for the eigenvector. If None, a random vector is generated.

Returns:

The approximated dominant eigenvalue of the matrix A, computed as the Rayleigh quotient x @ A @ x.

Return type:

float

Raises:

• TypeError – If the input is not a NumPy array or a SciPy/CuPy sparse matrix.

• ValueError – If number of rows != number of columns.

pyclassify.power_method_numba(A, max_iter=500, tol=1e-07, x=None)

Compute the dominant eigenvalue of a matrix using the power method, with Numba optimization.

This function and applies Numba’s Just-In-Time (JIT) compilation to optimize the performance of the power method for large matrices.

Remark that numba does not support scipy sparse matrices, so the input matrix must be a NumPy array. he function is optimized with Numba using the ‘njit’ decorator with nogil and parallel options. We have re-written the code due to the fact that using numba we cannot use the helper function check_square_matrix.

Parameters:

• A (np.ndarray) – A square matrix.

• max_iter (int, optional) – Maximum number of iterations to perform (default is 500).

• tol (float, optional) – Tolerance for convergence based on the relative change between iterations (default is 1e-4).

• x (np.ndarray, optional) – Initial guess for the eigenvector. If None, a random vector is generated.

Returns:

The approximated dominant eigenvalue of the matrix A.

Return type:

float

Raises:

ValueError – If the input matrix A is not square. The check is not done using ‘check_square_matrix’ because of numba technicalities.