Use cholespy for CHOLMOD instead of scikit-sparse

.github/workflows/pytest.yml

@@ -45,7 +45,7 @@ jobs:
       - name: Install package
         run: |
           python -m pip install --progress-bar off --upgrade pip setuptools wheel
-          python -m pip install --progress-bar off .[test]
+          python -m pip install --progress-bar off .[chol,test]
       - name: Display system information
         run: lapy-sys_info --developer
       - name: Run pytest
@@ -81,7 +81,7 @@ jobs:
         run: |
           python -m pip install --progress-bar off uv
           uv pip install --system pytest pytest-timeout
-          uv pip install --system --resolution lowest-direct .
+          uv pip install --system --resolution lowest-direct ".[chol]"
       - name: Display system information
         run: lapy-sys_info --developer
       - name: Run pytest

README.md

@@ -62,8 +62,7 @@ python3 -m pip install --user --src /my/preferred/location --editable git+https:
 
 ### Dependencies and optional backends
 - Core: Python 3, NumPy, SciPy.
-- Optional (recommended): scikit-sparse (for CHOLMOD) to accelerate Cholesky sparse solves.
-Note: CHOLMOD (via scikit-sparse) is not currently installable via plain `pip` on all platforms; use `conda` when possible. If `use_cholmod=True` is requested, LaPy attempts to import CHOLMOD and will raise if it is unavailable. Install ordering: install `numpy` and `scipy` first, then `scikit-sparse`.
+- Optional: cholespy (for CHOLMOD) to accelerate sparse solves with the Cholesky decomposition. Install it with `pip install lapy[chol]`. If `use_cholmod=True` is requested without cholespy, LaPy raises an ImportError.
 
 ## Solver backends
 - Default: SciPy sparse LU/QR routines.

lapy/_cholmod.py

@@ -0,0 +1,102 @@
+"""Sparse Cholesky factorization with CHOLMOD, provided by the cholespy package."""
+
+import numpy as np
+from scipy import sparse
+
+from .utils._imports import import_optional_dependency
+
+_INT32_MAX = np.iinfo(np.int32).max
+
+
+def require_cholmod() -> None:
+    """Raise an ImportError if cholespy, which provides CHOLMOD, is not installed."""
+    import_optional_dependency("cholespy", extra="It is required for use_cholmod=True.", raise_error=True)
+
+
+class CholeskyFactor:
+    """Cholesky factorization of a sparse symmetric positive definite matrix.
+
+    CHOLMOD reads only one triangle of the matrix, so the matrix has to be
+    genuinely symmetric. The factorization and the solves run in double
+    precision; solutions are returned in ``float32`` for a ``float32`` matrix
+    and in ``float64`` otherwise.
+
+    Parameters
+    ----------
+    matrix : scipy.sparse.spmatrix
+        Sparse symmetric positive definite matrix of shape (n, n).
+
+    Raises
+    ------
+    ImportError
+        If cholespy is not installed.
+    ValueError
+        If the matrix is not square, too large for 32-bit indices, or not
+        positive definite.
+    """
+
+    def __init__(self, matrix) -> None:
+        cholespy = import_optional_dependency(
+            "cholespy", extra="It is required for use_cholmod=True.", raise_error=True
+        )
+        a = sparse.csc_matrix(matrix)
+        a.sum_duplicates()
+        a.sort_indices()
+        n = a.shape[0]
+        if a.shape != (n, n):
+            raise ValueError(f"Cholesky factorization needs a square matrix, got shape {a.shape}.")
+        # cholespy takes 32-bit indices; scipy switches to 64-bit ones beyond that range
+        if n > _INT32_MAX or a.nnz > _INT32_MAX:
+            raise ValueError(
+                f"Matrix with {n} rows and {a.nnz} non-zeros is too large for CHOLMOD's 32-bit indices."
+            )
+        # output precision; cholespy factorizes in double precision on the CPU either way
+        self.dtype = np.dtype(np.float32) if a.dtype == np.float32 else np.dtype(np.float64)
+        # Passed as CSR, which for a symmetric matrix holds the same arrays as CSC.
+        # cholespy reallocates the arrays of a CSC input in place, which corrupts the
+        # heap when NumPy owns them (it crashes on Windows), but copies a CSR input.
+        try:
+            self._solver = cholespy.CholeskySolverD(
+                n,
+                a.indptr.astype(np.int32),
+                a.indices.astype(np.int32),
+                a.data.astype(np.float64),
+                cholespy.MatrixType.CSR,
+            )
+        except ValueError as err:
+            if "positive definite" not in str(err):
+                raise
+            raise ValueError(
+                "Cholesky factorization failed because the matrix is not positive definite, "
+                "for example a Laplace operator without boundary conditions. "
+                "Use use_cholmod=False for the LU decomposition instead."
+            ) from err
+        self.shape = (n, n)
+
+    def solve(self, b) -> np.ndarray:
+        """Solve ``A x = b`` with the factorized matrix.
+
+        Parameters
+        ----------
+        b : array_like or scipy.sparse.spmatrix
+            Right hand side of shape (n,) or (n, n_rhs).
+
+        Returns
+        -------
+        np.ndarray
+            Solution of the same shape as ``b``, ``float32`` for a ``float32``
+            matrix and ``float64`` otherwise.
+
+        Raises
+        ------
+        ValueError
+            If the first dimension of ``b`` does not match the matrix.
+        """
+        if sparse.issparse(b):
+            b = b.toarray()
+        b = np.ascontiguousarray(np.asarray(b), dtype=np.float64)
+        if b.ndim not in (1, 2) or b.shape[0] != self.shape[0]:
+            raise ValueError(f"Right hand side of shape {b.shape} does not match a matrix of shape {self.shape}.")
+        x = np.empty_like(b)
+        self._solver.solve(b, x)
+        return x.astype(self.dtype, copy=False)

lapy/conformal.py

@@ -16,7 +16,6 @@
 Closed Brain Surfaces."
 SIAM Journal on Imaging Sciences, vol. 8, no. 1, pp. 67-94, 2015.
 """
-import importlib
 import logging
 from typing import Any
 
@@ -25,7 +24,7 @@
 from scipy.optimize import minimize
 
 from . import Solver, TriaMesh
-from .utils._imports import import_optional_dependency
+from ._cholmod import CholeskyFactor
 
 logger = logging.getLogger(__name__)
 
@@ -178,8 +177,8 @@ def spherical_conformal_map(tria: TriaMesh, use_cholmod: bool = False) -> np.nda
     tria : TriaMesh
         A triangular mesh object representing a genus-0 closed surface.
     use_cholmod : bool, default=False
-        Which solver to use. If True, use Cholesky decomposition from
-        scikit-sparse cholmod. If False, use spsolve (LU decomposition).
+        Which solver to use. If True, use the sparse Cholesky decomposition of
+        CHOLMOD (requires ``cholespy``). If False, use the LU decomposition.
 
     Returns
     -------
@@ -194,7 +193,7 @@ def spherical_conformal_map(tria: TriaMesh, use_cholmod: bool = False) -> np.nda
         If edge lengths are degenerate.
         If projection contains NaN values.
     ImportError
-        If use_cholmod is True but scikit-sparse is not installed.
+        If use_cholmod is True but cholespy is not installed.
     """
     # Ensure the input mesh has genus-0 topology
     _ensure_genus_zero(tria)
@@ -436,8 +435,8 @@ def spherical_tutte_map(
     bigtri : int, default=0
         Index of the triangle to pin at the north pole.
     use_cholmod : bool, default=False
-        Which solver to use. If True, use Cholesky decomposition from
-        scikit-sparse cholmod. If False, use spsolve (LU decomposition).
+        Which solver to use. If True, use the sparse Cholesky decomposition of
+        CHOLMOD (requires ``cholespy``). If False, use the LU decomposition.
 
     Returns
     -------
@@ -449,7 +448,7 @@ def spherical_tutte_map(
     ValueError
         If mesh is not genus-0 (Euler characteristic != 2).
     ImportError
-        If use_cholmod is True but scikit-sparse is not installed.
+        If use_cholmod is True but cholespy is not installed.
     """
     _ensure_genus_zero(tria)
     return inverse_stereographic(
@@ -644,8 +643,8 @@ def linear_beltrami_solver(
     target : np.ndarray
         2D target positions for the landmark vertices, shape (n_landmarks, 2).
     use_cholmod : bool, default=False
-        Which solver to use. If True, use Cholesky decomposition from
-        scikit-sparse cholmod. If False, use spsolve (LU decomposition).
+        Which solver to use. If True, use the sparse Cholesky decomposition of
+        CHOLMOD (requires ``cholespy``). If False, use the LU decomposition.
 
     Returns
     -------
@@ -660,7 +659,7 @@ def linear_beltrami_solver(
         If triangle areas are degenerate.
         If Beltrami denominator is close to zero.
     ImportError
-        If use_cholmod is True but scikit-sparse is not installed.
+        If use_cholmod is True but cholespy is not installed.
     """
     # Ensure the triangulation is planar
     _ensure_planar_mesh(tria, "Linear Beltrami solver")
@@ -732,9 +731,8 @@ def _sparse_symmetric_solve(
 ) -> np.ndarray:
     """Solve the real sparse symmetric linear system of equations Ax = b.
 
-    Depending on the availability of the `scikit-sparse` package, it uses either:
-    - Cholesky decomposition (via scikit-sparse) for performance-optimal solving.
-    - LU decomposition (via SciPy) if scikit-sparse is not available.
+    Uses the sparse Cholesky decomposition of CHOLMOD (via ``cholespy``) if
+    ``use_cholmod`` is True, and the LU decomposition of SciPy otherwise.
 
     ``A`` has to be genuinely symmetric for the two branches to agree: the
     Cholesky solver reads only the lower triangular part of ``A``, whereas
@@ -748,8 +746,8 @@ def _sparse_symmetric_solve(
     b : np.ndarray
         Real right hand side of shape (n,) or (n, n_rhs).
     use_cholmod : bool, default=False
-        Which solver to use. If True, use Cholesky decomposition from
-        scikit-sparse cholmod. If False, use spsolve (LU decomposition).
+        Which solver to use. If True, use the sparse Cholesky decomposition of
+        CHOLMOD (requires ``cholespy``). If False, use the LU decomposition.
 
     Returns
     -------
@@ -761,12 +759,11 @@ def _sparse_symmetric_solve(
     ValueError
         If ``A`` or ``b`` is complex.
     ImportError
-        If use_cholmod is True but scikit-sparse is not installed.
+        If use_cholmod is True but cholespy is not installed.
     RuntimeError
         Propagated from ``splu`` when ``A`` is singular.
-    sksparse.cholmod.CholmodNotPositiveDefiniteError
-        Propagated from CHOLMOD when ``A`` is not positive definite and
-        use_cholmod is True.
+    ValueError
+        If use_cholmod is True and ``A`` is not positive definite.
     """
     if np.iscomplexobj(A) or np.iscomplexobj(b):
         raise ValueError(
@@ -776,11 +773,8 @@ def _sparse_symmetric_solve(
         )
     b = np.ascontiguousarray(b, dtype=np.float64)
     if use_cholmod:
-        sksparse = import_optional_dependency("sksparse", raise_error=True)
-        importlib.import_module(".cholmod", sksparse.__name__)
-        logger.info("Solver: Cholesky decomposition (scikit-sparse cholmod)")
-        chol = sksparse.cholmod.cholesky(sparse.csc_matrix(A))
-        x = chol(b)
+        logger.info("Solver: Cholesky decomposition (CHOLMOD)")
+        x = CholeskyFactor(A).solve(b)
     else:
         from scipy.sparse.linalg import splu
         logger.info("Solver: LU decomposition (spsolve)")

lapy/diffgeo.py

@@ -7,15 +7,14 @@
 tet and tria versions.
 """
 
-import importlib
 import logging
 from typing import TYPE_CHECKING, Union
 
 import numpy as np
 from scipy import sparse
 
 from . import Solver, TriaMesh
-from .utils._imports import import_optional_dependency
+from ._cholmod import CholeskyFactor, require_cholmod
 
 if TYPE_CHECKING:
     from .tet_mesh import TetMesh
@@ -132,8 +131,8 @@ def compute_geodesic_f(
         Scalar function at vertices, shape ``(n_vertices,)`` or
         ``(n_vertices, n_functions)``.
     use_cholmod : bool, default=False
-        If True, use Cholesky decomposition from scikit-sparse cholmod for the
-        linear solve. If False, use spsolve (LU decomposition).
+        If True, use the sparse Cholesky decomposition of CHOLMOD (requires
+        ``cholespy``) for the linear solve. If False, use the LU decomposition.
 
     Returns
     -------
@@ -185,8 +184,8 @@ def tria_compute_geodesic_f(
         Scalar function at vertices, shape ``(n_vertices,)`` or
         ``(n_vertices, n_functions)``.
     use_cholmod : bool, default=False
-        If True, use Cholesky decomposition from scikit-sparse cholmod for the
-        linear solve. If False, use spsolve (LU decomposition).
+        If True, use the sparse Cholesky decomposition of CHOLMOD (requires
+        ``cholespy``) for the linear solve. If False, use the LU decomposition.
 
     Returns
     -------
@@ -488,8 +487,8 @@ def tria_compute_rotated_f(
         Scalar function at vertices, shape ``(n_vertices,)`` or
         ``(n_vertices, n_functions)``.
     use_cholmod : bool, default=False
-        If True, use Cholesky decomposition from scikit-sparse cholmod for the
-        linear solve. If False, use spsolve (LU decomposition).
+        If True, use the sparse Cholesky decomposition of CHOLMOD (requires
+        ``cholespy``) for the linear solve. If False, use the LU decomposition.
 
     Returns
     -------
@@ -549,8 +548,8 @@ def tria_mean_curvature_flow(
     step : float, default=1.0
         Euler step size.
     use_cholmod : bool, default=False
-        Which solver to use. If True, use Cholesky decomposition from
-        scikit-sparse cholmod. If False, use spsolve (LU decomposition).
+        Which solver to use. If True, use the sparse Cholesky decomposition of
+        CHOLMOD (requires ``cholespy``). If False, use the LU decomposition.
 
     Returns
     -------
@@ -560,17 +559,14 @@ def tria_mean_curvature_flow(
     Raises
     ------
     ImportError
-        If use_cholmod is True but scikit-sparse is not installed.
+        If use_cholmod is True but cholespy is not installed.
 
     Notes
     -----
     Numexpr could speed up this functions if necessary.
     """
     if use_cholmod:
-        sksparse = import_optional_dependency("sksparse", raise_error=True)
-        importlib.import_module(".cholmod", sksparse.__name__)
-    else:
-        sksparse = None
+        require_cholmod()
     # pre-normalize
     trianorm = TriaMesh(tria.v, tria.t)
     trianorm.normalize_()
@@ -586,12 +582,9 @@ def tria_mean_curvature_flow(
         mass_v = mass.dot(trianorm.v)
         # solve (M + step*A) * v = Mv and update vertices
         if use_cholmod:
-            logger.info("Solver: Cholesky decomposition from scikit-sparse cholmod")
-            factor = sksparse.cholmod.cholesky(mass + step * a_mat)
-            trianorm.v = factor(mass_v)
+            logger.info("Solver: Cholesky decomposition (CHOLMOD)")
+            trianorm.v = CholeskyFactor(mass + step * a_mat).solve(mass_v)
         else:
-            # Note, it would be better to do sparse Cholesky (CHOLMOD)
-            # as it can be 5-6 times faster
             from scipy.sparse.linalg import spsolve
 
             logger.info("Solver: spsolve (LU decomposition)")