diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
index 6abf12a..1c3a299 100644
--- a/.github/workflows/pytest.yml
+++ b/.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
diff --git a/README.md b/README.md
index 12c9096..6121854 100644
--- a/README.md
+++ b/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.
diff --git a/lapy/_cholmod.py b/lapy/_cholmod.py
new file mode 100644
index 0000000..c2adc7b
--- /dev/null
+++ b/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)
diff --git a/lapy/conformal.py b/lapy/conformal.py
index ac5ef9b..2c8f75c 100644
--- a/lapy/conformal.py
+++ b/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)")
diff --git a/lapy/diffgeo.py b/lapy/diffgeo.py
index dfdc7fc..3eea007 100644
--- a/lapy/diffgeo.py
+++ b/lapy/diffgeo.py
@@ -7,7 +7,6 @@
 tet and tria versions.
 """
 
-import importlib
 import logging
 from typing import TYPE_CHECKING, Union
 
@@ -15,7 +14,7 @@
 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)")
diff --git a/lapy/heat.py b/lapy/heat.py
index 460805a..5d91459 100644
--- a/lapy/heat.py
+++ b/lapy/heat.py
@@ -6,12 +6,11 @@
 
 from __future__ import annotations
 
-import importlib
 import logging
 
 import numpy as np
 
-from .utils._imports import import_optional_dependency
+from ._cholmod import CholeskyFactor, require_cholmod
 
 logger = logging.getLogger(__name__)
 
@@ -158,8 +157,8 @@ def diffusion(
     aniso : int, default=None
         Number of smoothing iterations for curvature computation on vertices.
     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
     -------
@@ -172,13 +171,10 @@ def diffusion(
     ValueError
         If vids contains out-of-range vertex indices.
     ImportError
-        If use_cholmod is True but scikit-sparse is not installed.
+        If use_cholmod is True but cholespy is not installed.
     """
     if use_cholmod:
-        sksparse = import_optional_dependency("sksparse", raise_error=True)
-        importlib.import_module(".cholmod", sksparse.__name__)
-    else:
-        sksparse = None
+        require_cholmod()
     from . import Solver
 
     nv = len(geometry.v)
@@ -213,12 +209,11 @@ def diffusion(
     for k, v in enumerate(vids_list):
         b0[v, k] = 1.0
 
-    # Solve — both splu.solve and cholmod accept 2-D rhs natively
+    # Solve — both splu.solve and the Cholesky factor accept 2-D rhs
     logger.debug("Matrix Format: %s", hmat.getformat())
     if use_cholmod:
-        logger.info("Solver: Cholesky decomposition from scikit-sparse cholmod")
-        chol = sksparse.cholmod.cholesky(hmat)
-        vfunc = chol(b0)
+        logger.info("Solver: Cholesky decomposition (CHOLMOD)")
+        vfunc = CholeskyFactor(hmat).solve(b0)
     else:
         from scipy.sparse.linalg import splu
 
diff --git a/lapy/shapedna.py b/lapy/shapedna.py
index dee85c4..fce95a6 100644
--- a/lapy/shapedna.py
+++ b/lapy/shapedna.py
@@ -122,9 +122,9 @@ def compute_shapedna(
         Number of smoothing iterations for curvature computation on vertices.
         See `lapy.Solver` class.
     use_cholmod : bool, default=False
-        If True, attempts to use the Cholesky decomposition for improved
-        execution speed. Requires the ``scikit-sparse`` library. If it cannot
-        be found, an error will be thrown. If False, will use slower LU
+        If True, use the sparse Cholesky decomposition of CHOLMOD for improved
+        execution speed. Requires the ``cholespy`` package; if it is not
+        installed, an ImportError is raised. If False, use the slower LU
         decomposition.
 
     Returns
diff --git a/lapy/solver.py b/lapy/solver.py
index 2ff8eda..bd943cd 100644
--- a/lapy/solver.py
+++ b/lapy/solver.py
@@ -1,4 +1,3 @@
-import importlib
 import inspect
 import logging
 import sys
@@ -7,9 +6,9 @@
 from scipy import sparse
 from scipy.sparse.linalg import eigsh as _eigsh
 
+from ._cholmod import CholeskyFactor, require_cholmod
 from .tet_mesh import TetMesh
 from .tria_mesh import TriaMesh
-from .utils._imports import import_optional_dependency
 
 logger = logging.getLogger(__name__)
 
@@ -41,10 +40,10 @@ class Solver:
     aniso_smooth : int | None
         Number of smoothing iterations for curvature computation on vertices.
     use_cholmod : bool, default: False
-        If True, attempts to use the Cholesky decomposition for improved execution
-        speed. Requires the ``scikit-sparse`` library. If it can not be found, an error
-        will be thrown.
-        If False, will use slower LU decomposition.
+        If True, use the sparse Cholesky decomposition of CHOLMOD for improved
+        execution speed. Requires the ``cholespy`` package; if it is not installed,
+        an ImportError is raised.
+        If False, use the slower LU decomposition.
     dtype : numpy.dtype, default=np.float64
         Data type for the stiffness and mass matrices. Defaults to ``np.float64``
         for numerical stability. Use ``np.float32`` explicitly if memory is constrained
@@ -66,10 +65,7 @@ def __init__(
         dtype: np.dtype = np.float64,
     ) -> None:
         if use_cholmod:
-            self.sksparse = import_optional_dependency("sksparse", raise_error=True)
-            importlib.import_module(".cholmod", self.sksparse.__name__)
-        else:
-            self.sksparse = None
+            require_cholmod()
 
         # Ensure dtype is a numpy dtype object
         dtype = np.dtype(dtype)
@@ -746,10 +742,10 @@ def eigs(
             start_kwargs = {"v0": v0}
 
         if self.use_cholmod:
-            logger.info("Solver: Cholesky decomposition from scikit-sparse cholmod ...")
-            chol = self.sksparse.cholmod.cholesky(self.stiffness - sigma * self.mass)
+            logger.info("Solver: Cholesky decomposition (CHOLMOD) ...")
+            chol = CholeskyFactor(self.stiffness - sigma * self.mass)
             op_inv = LinearOperator(
-                matvec=chol,
+                matvec=chol.solve,
                 shape=self.stiffness.shape,
                 dtype=self.stiffness.dtype,
             )
@@ -936,12 +932,11 @@ def poisson(
                 raise ValueError("A matrix needs to be sparse CSC or CSR")
         else:
             a = self.stiffness
-        # solve A x = b  (both splu.solve and cholmod accept 2-D b natively)
+        # solve A x = b  (both splu.solve and the Cholesky factor accept 2-D b)
         logger.debug("Matrix format now: %s", a.getformat())
         if self.use_cholmod:
-            logger.info("Solver: Cholesky decomposition from scikit-sparse cholmod ...")
-            chol = self.sksparse.cholmod.cholesky(a)
-            x = chol(b.astype(dtype))
+            logger.info("Solver: Cholesky decomposition (CHOLMOD) ...")
+            x = CholeskyFactor(a).solve(b)
         else:
             from scipy.sparse.linalg import splu
 
diff --git a/lapy/utils/_imports.py b/lapy/utils/_imports.py
index fb53167..6e2da70 100644
--- a/lapy/utils/_imports.py
+++ b/lapy/utils/_imports.py
@@ -7,9 +7,7 @@
 
 # A mapping from import name to package name (on PyPI) when the package name
 # is different.
-INSTALL_MAPPING = {
-    "sksparse": "scikit-sparse",
-}
+INSTALL_MAPPING: dict[str, str] = {}
 
 
 def import_optional_dependency(
diff --git a/lapy/utils/tests/test_TetMesh_Geodesics.py b/lapy/utils/tests/test_TetMesh_Geodesics.py
index ad142b6..7141a88 100644
--- a/lapy/utils/tests/test_TetMesh_Geodesics.py
+++ b/lapy/utils/tests/test_TetMesh_Geodesics.py
@@ -121,26 +121,15 @@ def test_tetMesh_Geodesics_format(load_tet_mesh, loaded_data):
     X = np.nan_to_num(X)
     divx = compute_divergence(T, X)
 
-    useCholmod = True
-    try:
-        from sksparse.cholmod import cholesky
-    except ImportError:
-        useCholmod = False
-
     fem = Solver(T, lump=True)
     A, B = fem.stiffness, fem.mass
 
     H = A
     b0 = -divx
 
-    if useCholmod:
-        print("Solver: cholesky decomp - performance optimal ...")
-        chol = cholesky(H)
-        x = chol(b0)
-    else:
-        print("Solver: spsolve (LU decomp) - performance not optimal ...")
-        lu = splu(H)
-        x = lu.solve(b0)
+    # H is singular without boundary conditions, so a Cholesky factorization
+    # does not apply
+    x = splu(H).solve(b0)
 
     x = x - np.min(x)
 
@@ -156,7 +145,6 @@ def test_tetMesh_Geodesics_format(load_tet_mesh, loaded_data):
     ]
     assert H.getformat() == expected_matrix_format
     assert np.shape(x) == (9261,)
-    assert not useCholmod, "Solver: cholesky decomp - performance optimal ..."
     expected_max_x = loaded_data["expected_outcomes"]["test_TetMesh_Geodesics"][
         "max_distance"
     ]
diff --git a/lapy/utils/tests/test_TriaMesh_Geodesics.py b/lapy/utils/tests/test_TriaMesh_Geodesics.py
index f4fb355..83097ae 100644
--- a/lapy/utils/tests/test_TriaMesh_Geodesics.py
+++ b/lapy/utils/tests/test_TriaMesh_Geodesics.py
@@ -141,28 +141,14 @@ def test_Geodesics_format(loaded_data, load_square_mesh):
     divx = compute_divergence(T, X)
     # compute distance
 
-    useCholmod = True
-    try:
-        from sksparse.cholmod import cholesky
-    except ImportError:
-        useCholmod = False
-
     fem = Solver(T, lump=True)
     A, B = fem.stiffness, fem.mass
     H = -A
     b0 = divx
 
-    # solve H x = b0
-    # we don't need the B matrix here, as divx is the integrated divergence
-    print("Matrix Format now: " + H.getformat())
-    if useCholmod:
-        print("Solver: cholesky decomp - performance optimal ...")
-        chol = cholesky(H)
-        x = chol(b0)
-    else:
-        print("Solver: spsolve (LU decomp) - performance not optimal ...")
-        lu = splu(H)
-        x = lu.solve(b0)
+    # solve H x = b0; H is singular without boundary conditions, so a Cholesky
+    # factorization does not apply, and B is not needed as divx is integrated
+    x = splu(H).solve(b0)
 
     # remove shift
     x = x - min(x)
@@ -175,7 +161,6 @@ def test_Geodesics_format(loaded_data, load_square_mesh):
         "expected_matrix_format"
     ]
     assert H.getformat() == expected_matrix_format
-    assert not useCholmod, "Solver: cholesky decomp - performance optimal ..."
     expected_max_x = loaded_data["expected_outcomes"]["test_Geodesics_format"][
         "max_distance"
     ]
diff --git a/lapy/utils/tests/test_cholmod.py b/lapy/utils/tests/test_cholmod.py
new file mode 100644
index 0000000..4b19ef8
--- /dev/null
+++ b/lapy/utils/tests/test_cholmod.py
@@ -0,0 +1,76 @@
+"""The CHOLMOD backend (use_cholmod=True) must agree with the LU backend everywhere it is used."""
+
+import numpy as np
+import pytest
+
+from ...diffgeo import tria_mean_curvature_flow
+from ...heat import diffusion
+from ...solver import Solver
+from ...tria_mesh import TriaMesh
+
+pytest.importorskip("cholespy")
+
+
+@pytest.fixture
+def cube_tria():
+    """Closed triangle surface of a cube."""
+    return TriaMesh.read_vtk("data/cubeTria.vtk")
+
+
+@pytest.fixture
+def square():
+    """Open planar triangulation of the unit square."""
+    return TriaMesh.read_off("data/square-mesh.off")
+
+
+def test_eigs_match_lu(cube_tria):
+    """Shift-invert eigenvalues agree, with the factor used as a LinearOperator inside ARPACK."""
+    lu = Solver(cube_tria, use_cholmod=False).eigs(k=8, rng=0)[0]
+    chol = Solver(cube_tria, use_cholmod=True).eigs(k=8, rng=0)[0]
+    np.testing.assert_allclose(chol, lu, rtol=1e-8, atol=1e-10)
+
+
+def test_poisson_matches_lu(square):
+    """Poisson solutions with Dirichlet and Neumann conditions agree, for one and several right hand sides."""
+    rng = np.random.default_rng(0)
+    nv = len(square.v)
+    bnd = np.concatenate(square.boundary_loops())
+    dtup = (bnd[::2], np.sin(square.v[bnd[::2], 0]))
+    ntup = (bnd[1::2], np.cos(square.v[bnd[1::2], 1]))
+    for h in (rng.normal(size=nv), rng.normal(size=(nv, 3))):
+        lu = Solver(square, use_cholmod=False).poisson(h, dtup, ntup)
+        chol = Solver(square, use_cholmod=True).poisson(h, dtup, ntup)
+        assert chol.shape == lu.shape == h.shape
+        np.testing.assert_allclose(chol, lu, rtol=1e-8, atol=1e-10)
+
+
+def test_poisson_with_float32_solver(square):
+    """A float32 solver returns float32 output close to the double precision result."""
+    bnd = np.concatenate(square.boundary_loops())
+    dtup = (bnd, np.zeros(len(bnd)))
+    h = np.ones(len(square.v))
+    ref = Solver(square, use_cholmod=False).poisson(h, dtup)
+    chol = Solver(square, use_cholmod=True, dtype=np.float32).poisson(h, dtup)
+    assert chol.dtype == np.float32
+    np.testing.assert_allclose(chol, ref, rtol=1e-4, atol=1e-6)
+
+
+def test_diffusion_matches_lu(cube_tria):
+    """Heat diffusion from several seed sets agrees."""
+    seeds = [[0, 1], [5], [10, 20, 30]]
+    lu = diffusion(cube_tria, seeds, use_cholmod=False)
+    chol = diffusion(cube_tria, seeds, use_cholmod=True)
+    np.testing.assert_allclose(chol, lu, rtol=1e-8, atol=1e-10)
+
+
+def test_mean_curvature_flow_matches_lu(cube_tria):
+    """Mean curvature flow agrees step by step."""
+    lu = tria_mean_curvature_flow(cube_tria, max_iter=5, use_cholmod=False)
+    chol = tria_mean_curvature_flow(cube_tria, max_iter=5, use_cholmod=True)
+    np.testing.assert_allclose(chol.v, lu.v, rtol=1e-8, atol=1e-10)
+
+
+def test_not_positive_definite_is_reported(cube_tria):
+    """Without boundary conditions the Laplace operator of a closed surface is singular."""
+    with pytest.raises(ValueError, match="not positive definite"):
+        Solver(cube_tria, use_cholmod=True).poisson(np.ones(len(cube_tria.v)))
diff --git a/lapy/utils/tests/test_conformal.py b/lapy/utils/tests/test_conformal.py
index 68c1e61..a4cd582 100644
--- a/lapy/utils/tests/test_conformal.py
+++ b/lapy/utils/tests/test_conformal.py
@@ -238,7 +238,7 @@ def poisoned(z, tria, bigtri):
 
 def test_cholmod_matches_lu(sphere):
     """The Cholesky and the LU backend must agree."""
-    pytest.importorskip("sksparse.cholmod")
+    pytest.importorskip("cholespy")
 
     lu = spherical_conformal_map(sphere, use_cholmod=False)
     chol = spherical_conformal_map(sphere, use_cholmod=True)
diff --git a/pyproject.toml b/pyproject.toml
index 0f03088..36b3801 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -53,7 +53,7 @@ build = [
     'twine',
 ]
 chol = [
-    'scikit-sparse>=0.4.16',
+    'cholespy>=2.0',
 ]
 doc = [
     'furo!=2023.8.17',