DOCSEARCHbyAlgolia
Submit Docs

PyPardiso

repository·master·Indexed 41 days ago

https://github.com/haasad/pypardiso

A Python interface to the Intel oneAPI Math Kernel Library (MKL) PARDISO solver for solving large sparse linear systems of equations. It provides a high-performance, multi-threaded alternative to SciPy's default sparse solvers, featuring a PyPardisoSolver class for granular control and SciPy-compatible aliases like spsolve and factorized.

Tokens
2.9K
Snippets
9
Records
15
Agent score
52%

What's inside PyPardiso

  1. PyPardiso limitations and alternatives

    master

    Before using PyPardiso, be aware of the following constraints:

    • Complex Numbers: PyPardiso does not currently support complex numbers. For advanced use cases involving complex numbers, consider Pardiso.jl from the JuliaSparse project.
    • macOS (Apple Silicon): MKL is not available on Apple silicon. For fast sparse solving on macOS, it is recommended to use scikit-umfpack instead.
    • Solver Identity: PyPardiso is NOT a Python interface to the PARDISO 7.2 Solver Project; it specifically uses the Intel oneAPI MKL PARDISO solver.

  2. Install PyPardiso via conda or pip

    master

    PyPardiso can be installed using either conda (from the conda-forge channel) or pip. It is recommended to use a virtual environment for installation. Note that PyPardiso runs on Linux and Windows, but is not recommended for macOS users on Apple silicon (use scikit-umfpack instead).

    # Using conda
conda install -c conda-forge pypardiso

# Using pip
pip install pypardiso

  3. Use PyPardisoSolver to solve sparse linear systems

    master

    The PyPardisoSolver class provides a Python interface to the Intel MKL PARDISO library for solving large sparse linear systems of equations $Ax=b$.

    Basic Usage

    By default, the solver is configured for real (float64) and non-symmetric matrices (mtype=11).

    • Solve once: Use solve(A, b) where A is a sparse CSR or CSC matrix and b is a numpy array.
    • Solve multiple times with the same matrix: Call factorize(A) first. This performs the factorization step once, allowing subsequent calls to solve(A, b) with different right-hand sides b to be significantly faster as the factorization is reused automatically.

    Requirements

    • Matrix Format: A must be a scipy.sparse.csr_matrix or scipy.sparse.csc_matrix with dtype=np.float64.
    • Matrix Properties: A must be square and must not contain empty rows (for CSR) or empty columns (for CSC), as this will cause a segmentation fault.
    • Right-hand side: b should be a dense numpy array. If a sparse array is provided, it will be converted to dense automatically (with a warning).
    import numpy as np
import scipy.sparse as sp
from pypardiso import PyPardisoSolver

# Create a sparse CSR matrix
A = sp.csr_matrix([[3, 1], [1, 2]], dtype=np.float64)
b = np.array([9, 8], dtype=np.float64)

# Initialize solver
solver = PyPardisoSolver()

# Option 1: Direct solve
x = solver.solve(A, b)

# Option 2: Efficient solve for multiple b vectors
solver.factorize(A)
x1 = solver.solve(A, b)
x2 = solver.solve(A, np.array([1, 1], dtype=np.float64))

  5. Solve sparse linear systems with pypardiso.spsolve()

    master

    Use pypardiso.spsolve(A, b) to solve the sparse linear system $Ax=b$ for $x$.

    Requirements:

    • A must be a square, sparse matrix in CSR or CSC format.
    • b can be a vector or a matrix.

    PyPardiso provides functionality similar to scipy.sparse.linalg.spsolve but utilizes the Intel oneAPI Math Kernel Library (MKL) PARDISO solver, which is a shared-memory multiprocessing parallel direct sparse solver. It is often significantly faster than SciPy's default single-threaded SuperLU solver.

    import pypardiso
import numpy as np
import scipy.sparse as sp

# Create a square, sparse matrix A in CSR format
A = sp.rand(10, 10, density=0.5, format='csr')

# Create a vector b
b = np.random.rand(10)

# Solve Ax = b
x = pypardiso.spsolve(A, b)

  6. Set matrix type and solving phase

    master

    Advanced control over the solver's behavior can be achieved through these methods:

    • set_matrix_type(mtype): Changes the matrix type (e.g., for symmetric matrices). Note that complex matrix types are currently not supported; the library is primarily tested for mtype=11 (real and non-symmetric).
    • set_phase(phase): Sets the specific PARDISO solving phase (e.g., factorization, solve, or analysis).

  7. Manage memory and factorization in PyPardisoSolver

    master

    To prevent memory leaks or to free up resources after large computations, use the following methods:

    • remove_stored_factorization(): Deletes the wrapper's copy of matrix A. This frees Python-side memory, though the factorization remains accessible in the underlying MKL solver via a direct call with phase=33.
    • free_memory(everything=False): Releases the internal memory of the MKL solver.
      • If everything=False (default), it only releases the factorization (LU-decomposition) memory.
      • If everything=True, it releases all internal MKL memory.

  8. Solve linear systems with spsolve and factorized aliases

    master

    PyPardiso provides functional aliases that mimic the SciPy sparse linear algebra API, making it easy to drop into existing SciPy workflows.

    • spsolve: Solves the linear system $Ax = b$.
    • factorized: Returns a callable that solves the system using a pre-computed factorization of $A$.
    • ps: A short alias for pypardiso_solver.
    from pypardiso import spsolve, factorized, ps

# Using spsolve
# x = spsolve(A, b)

# Using factorized
# solve = factorized(A)
# x = solve(b)

  9. Configure PyPardisoSolver via iparms

    master

    You can fine-tune the PARDISO solver behavior using the iparm array, which corresponds to the internal PARDISO control parameters.

    • Accessing parameters: Use get_iparm(i) to get the $i$-th parameter or get_iparms() to get a dictionary of all parameters (using 1-based indexing).
    • Setting parameters: Use set_iparm(i, value) to set the $i$-th parameter (using 1-based indexing).

    Warning: Be cautious with iparm[60] (Out-of-Core/OOC solver), as it may cause system instability if not managed correctly.

    Supported iparm indices for input include: 1, 2, 4, 5, 6, 8, 10, 11, 12, 13, 18, 19, 21, 24, 25, 27, 28, 31, 34, 35, 36, 37, 56, 60.

  10. Use factorized() for compatibility with SciPy

    master

    The factorized function mimics scipy.sparse.linalg.factorized. It returns a callable that, when passed an array b, returns the solution to $Ax = b$.

    Parameters:

    • A: Sparse square CSR or CSC matrix.
    • solver: An instance of PyPardisoSolver (default is the internal pypardiso_solver).

    Returns:

    • solve_b: A callable. Passing an array b to this callable returns the solution.

    Important Usage Notes:

    • Performance: Unlike SciPy, there is no performance gain in PyPardiso by using factorized instead of spsolve. It is recommended to use spsolve directly whenever possible.
    • Memory: The returned callable stores a copy of matrix A to ensure correctness. Do not use factorized for very large matrices as this increases memory consumption.
    • Factorization: The factorization is not stored within the callable; the first call to the returned function will take longer as it performs the factorization.
    import scipy.sparse as sp
import numpy as np
from pypardiso import factorized

A = sp.csr_matrix([[3, 1], [1, 2]])
b = np.array([9, 8])

# Get the callable
solve_b = factorized(A)

# Use the callable to solve
x = solve_b(b)
print(x)

  11. Solve linear systems with spsolve()

    master

    The spsolve function mimics scipy.sparse.linalg.spsolve but uses the Pardiso solver for improved performance. It solves the linear system $Ax = b$.

    Parameters:

    • A: Sparse square CSR or CSC matrix (scipy.sparse.csr_matrix or csc_matrix). Note: If a CSC matrix is provided, it is internally converted to CSR.
    • b: numpy.ndarray representing the right-hand side(s). b.shape[0] must match A.shape[0].
    • factorize: bool (default True). If True, the matrix $A$ is factorized so the factorization can be reused. Subsequent calls to spsolve with the same matrix $A$ will be significantly faster.
    • squeeze: bool (default True). Mimics SciPy's behavior by returning x.squeeze(). If False, the returned shape matches b.
    • solver: An instance of PyPardisoSolver (default is the internal pypardiso_solver).

    Returns:

    • x: numpy.ndarray containing the solution. If squeeze=True and b has shape (n, 1), the result will have shape (n,).

    Warning: Avoid using multiple instances of PyPardisoSolver in parallel (e.g., in different threads), as this can cause Python to crash on Windows.

    import numpy as np
import scipy.sparse as sp
from pypardiso import spsolve

# Create a sparse matrix A and right-hand side b
A = sp.csr_matrix([[3, 1], [1, 2]])
b = np.array([9, 8])

# Solve Ax = b
x = spsolve(A, b)
print(x)

  12. PyPardisoSolver API Reference

    master

    Class: PyPardisoSolver

    Initialization

    __init__(self, mtype=11, phase=13, size_limit_storage=5e7)

    • mtype (int): Matrix type (default 11 for real non-symmetric).
    • phase (int): Initial phase (default 13).
    • size_limit_storage (float): Threshold for using hashing instead of full matrix storage for factorization tracking.

    Methods

    • factorize(A): Factorizes sparse square CSR or CSC matrix A.
    • solve(A, b): Solves $Ax=b$ for $x$. Returns a numpy ndarray.
    • get_iparm(i): Returns the $i$-th iparm (1-based indexing).
    • get_iparms(): Returns a dictionary of all iparms (1-based indexing).
    • set_iparm(i, value): Sets the $i$-th iparm (1-based indexing).
    • set_matrix_type(mtype): Sets the matrix type.
    • set_phase(phase): Sets the solving phase.
    • set_statistical_info_on(): Enables printing of statistical info to stdout.
    • set_statistical_info_off(): Disables statistical info.
    • remove_stored_factorization(): Removes the stored factorization copy of matrix A.
    • free_memory(everything=False): Releases internal MKL memory.