PyClassify Module Documentation

This module contains all the functions for the PyClassify package.

Functions:

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.

Submodule: utils

pyclassify.utils.check_square_matrix(A)

Checks if the input matrix is a square matrix of type NumPy ndarray or SciPy 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). 2. A square matrix.

Parameters:

A (np.ndarray or sp.spmatrix or cpsp.spmatrix) – The matrix to be checked.

Raises:

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

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

pyclassify.utils.check_symm_square(A)

Checks if the input matrix is a square symmetric matrix of type NumPy ndarray or SciPy sparse matrix. It uses check_square_matrix.

This function performs the following steps: 1. Verifies the input type is either np.ndarray or sp.spmatrix. 2. Checks that the matrix is square. 3. Validates that the matrix is symmetric.

Parameters:

A (np.ndarray or sp.spmatrix) – The matrix to be validated.

Raises:

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

• ValueError – If the matrix is not square or not symmetric.

pyclassify.utils.make_symmetric(A)

Ensures that the input matrix is symmetric by averaging it with its transpose.

This function performs the following steps: 1. Checks if the matrix is a square matrix using check_square_matrix. 2. Computes the symmetric version of the matrix using the formula: (A + A.T) / 2.

Parameters:

A (np.ndarray or sp.spmatrix) – The input matrix to be symmetrized.

Returns:

The symmetric version of the input matrix.

Return type:

np.ndarray or sp.spmatrix

Raises:

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

• ValueError – If the matrix is not square.

pyclassify.utils.max_iteration_warning()

Prints a warning message indicating that the maximum number of iterations has been reached.

This function is used to alert the user that the iterative method likely did not converge, and suggests that a lower tolerance or alternative approach may be needed.

pyclassify.utils.poisson_2d_structure(n, k=None)

Constructs a sparse matrix of size n x n with the same 5-diagonal structure as the 2D Poisson matrix (main, ±1, ±k). Used for testing structure, not actual PDEs: in real life, this function would return a matrix of size $n^2$.

Parameters:

• n (int) – Size of the square matrix.

• k (int or None) – The offset for the “long-range” diagonals. If None, uses k = int(sqrt(n)).

Returns:

A sparse matrix with 5 diagonals for testing.

Return type:

scipy.sparse.spmatrix

pyclassify.utils.profile_numpy_eigvals(A)

Profiles the memory usage of computing eigenvalues and eigenvectors using NumPy.

This function performs the following steps: 1. Measures memory before the eigendecomposition. 2. Computes eigenvalues and eigenvectors using np.linalg.eigh. 3. Measures memory after the computation. 4. Returns the difference in memory usage.

Parameters:

A (np.ndarray) – A symmetric NumPy array.

Returns:

Memory used in MB during the computation.

Return type:

float

pyclassify.utils.profile_scipy_eigvals(A)

Profiles the memory usage of computing eigenvalues and eigenvectors using SciPy.

This function performs the following steps: 1. Measures memory before the eigendecomposition. 2. Computes eigenvalues and eigenvectors using scipy.linalg.eigh. 3. Measures memory after the computation. 4. Returns the difference in memory usage.

Parameters:

A (np.ndarray) – A symmetric NumPy array.

Returns:

Memory used in MB during the computation.

Return type:

float

pyclassify.utils.read_config(file: str) → dict

Reads a YAML configuration file and loads its contents into a dictionary.

This function performs the following steps: 1. Constructs the absolute path to the YAML file by appending ‘.yaml’ to the given base name. 2. Opens and parses the file using yaml.safe_load.

Parameters:

file (str) – The base name of the YAML file (without extension).

Returns:

A dictionary containing the configuration parameters.

Return type:

dict

Submodule: helpers_secular

pyclassify.zero_finder.bisection(f, a, b, tol, max_iter)

Standard bisection method to find a root of the function f in the interval [a, b].

This implementation is used in compute_outer_zero to locate the outer eigenvalue.

Parameters: f (callable): A continuous function for which f(a) * f(b) < 0. a (float): Left endpoint of the interval. b (float): Right endpoint of the interval. tol (float): Tolerance for convergence. The method stops when the interval is smaller than tol or when f(c) is sufficiently small. max_iter (int): Maximum number of iterations before stopping.

Returns: float: Approximation of the root within the specified tolerance.

pyclassify.zero_finder.compute_Psi(i, v, d, rho)

This function computes the functions Psi1, and Psi2 and their derivative, accordingly to the equation described by the notes chapter 5-6 contained in the shared folders. The equation utilized are the equation 5.27, 5.25 and 5.29. d is assumed to be already a vector whose elements are already sorted. Inputs:

-i: It is the index of the i-th smallest eigenvalue that is being computed -v: Rank one correction -d: Diagonal element of the Diagonal matrix described in equation 5.7 -rho: off diagonal element used for the splitting.

Output:

Psi1: Lamoba function that returns the value of Psi1 at the point x Psi2: Lamoba function that returns the value of Psi2 at the point x dPsi1: Lamoba function that returns the value of the derivative dPsi1/dx at the point x dPsi2: Lamoba function that returns the value of the derivative dPsi2/dx at the point x

pyclassify.zero_finder.compute_outer_zero(v, d, rho, interval_end, tol=1e-12, max_iter=1000)

Computes the outer eigenvalue (lambda[0] if rho < 0, lambda[n-1] if rho > 0) of a rank-one modified diagonal matrix.

The secular function behaves such that:

• If rho > 0, the outer eigenvalue lies in (d[n-1], infty), and f is increasing in this interval.

• If rho < 0, the outer eigenvalue lies in (-infty, d[0]), and f is decreasing in this interval.

This function: 1. Determines the direction to search based on the sign of rho. 2. Finds an upper bound (or lower bound for rho < 0) where the secular function changes sign. 3. Uses the bisection method to find the root in the determined interval.

Returns: float: Approximation of the outer eigenvalue.

pyclassify.zero_finder.find_root(i, left_center, v, d, rho, lam_0, tol=1e-15, maxiter=100)

Find the roots of the secular equation contained inside the interval min(d)=d[0] and max(d)=d[-1]. Inputs:

-i: It is the index of the i-th smallest eigenvalue that is being computed -left_center: boolen variable. True if the left extreme is the center od the new - shifted reference system. -v: Rank one correction -d: Diagonal element of the Diagonal matrix described in equation 5.7 -rho: off diagonal element used for the splitting. -lam_0: initial guess of the i-th root of the secular equation -tol: absolute or relative tollerence. For lamba value below 1e-6, tols refers to the absolute tolerance (avoiding division by zero). For

lambda value greater than 1e-6, tol is the absolute tolerance.

-maxiter: Maximum number of iteration of the iterative solver.

Outputs:

-shift + lam_0: i-th smallest eigenvalue -lam_0: difference between the i-th eigenvalue and the nearest diagonal element.

pyclassify.zero_finder.secular_solver_python(rho, d, v)

Computes all the roots of the secular equation. Inputs:

-v: Rank one correction -d: Diagonal element of the Diagonal matrix described in equation 5.7 -rho: off diagonal element used for the splitting.

Outputs:

-eig_val: vector of the roots of the secular equation -index: vector of the index of the center of the shifted frame of reference associated to the i-ith root -delta: vector of difference between the i-th eigenvalue and the nearest diagonal element.

Submodule: parallel_tridiag_eigen

pyclassify.parallel_tridiag_eigen.check_column_directions(A, B)

pyclassify.parallel_tridiag_eigen.find_interval_extreme(total_dimension, n_processor)

Computes the intervals for vector for being scattered.

Input:

-total_dimension: the dimension of the vector that has to be splitted -n_processor: the number of processor to which the scatter vector has to be sent

pyclassify.parallel_tridiag_eigen.parallel_eigen(main_diag, off_diag, tol_QR=1e-15, max_iterQR=5000, tol_deflation=1e-15)

pyclassify.parallel_tridiag_eigen.parallel_tridiag_eigen(diag, off, comm=None, tol_factor=1e-16, min_size=1, depth=0, profiler=None, tol_QR=1e-08, max_iterQR=5000)

Computes eigenvalues and eigenvectors of a symmetric tridiagonal matrix.

Input: -diag: diagonal part of the tridiagonal matrix -off: off diagonal part of the tridiagonal matrix -comm: MPI communicator -tol_factor: tollerance for the deflating step

Output: -final_eigvals: return the eigenvalues of the tridiagonal matrix -final_eigvecs: return the eigenvectors of the tridiagonal matrix