From 504e0fe3327dce25fc4afb10b1dd3f623ca31b0a Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Sun, 15 Mar 2026 09:56:31 +0100 Subject: [PATCH 1/5] expose eigsh paramters --- lapy/solver.py | 64 ++++++++++++++++++++++++++++----- lapy/utils/tests/test_solver.py | 27 +++++++++++++- 2 files changed, 82 insertions(+), 9 deletions(-) diff --git a/lapy/solver.py b/lapy/solver.py index f4b02dc..d064f87 100644 --- a/lapy/solver.py +++ b/lapy/solver.py @@ -664,7 +664,15 @@ def _fem_voxels( b = sparse.csc_matrix((local_b, (i, j)), dtype=dtype) return a, b - def eigs(self, k: int = 10, sigma: float = -0.01) -> tuple[np.ndarray, np.ndarray]: + def eigs( + self, + k: int = 10, + sigma: float = -0.01, + which: str = "LM", + v0: np.ndarray | None = None, + mode: str = "normal", + rng: int | np.random.Generator | None = None, + ) -> tuple[np.ndarray, np.ndarray]: """Compute the linear finite-element method Laplace-Beltrami spectrum. Parameters @@ -675,15 +683,39 @@ def eigs(self, k: int = 10, sigma: float = -0.01) -> tuple[np.ndarray, np.ndarra compute all eigenvectors of a matrix. sigma : float, default=-0.01 Shift value for the shift-invert mode. The solver finds eigenvalues - near sigma. Negative values work well for finding smallest eigenvalues. - Adjust if convergence issues occur (typically small negative). + near ``sigma``. The default small negative value works well for + finding the smallest non-negative Laplacian eigenvalues. Adjust if + convergence issues occur. The returned eigenvalues are always sorted + in ascending order regardless of ``sigma``. + which : str, default="LM" + Which ``k`` eigenvalues to find. With shift-invert (``sigma`` set), + ``"LM"`` selects the ``k`` eigenvalues of the original problem + closest to ``sigma`` — this is the standard choice for computing + the smallest Laplacian eigenvalues. Other options accepted by + ``scipy.sparse.linalg.eigsh`` are ``"SM"``, ``"LA"``, ``"SA"``, + and ``"BE"``. + v0 : np.ndarray, shape (n_vertices,), optional + Starting vector for the ARPACK iteration. Providing a fixed vector + makes results reproducible. Takes precedence over ``rng``. If both + are ``None``, ARPACK uses its own random starting vector. + mode : str, default="normal" + Shift-invert mode passed to ``scipy.sparse.linalg.eigsh``. + ``"normal"`` is the standard choice for Laplacian spectra. + Other valid values are ``"buckling"`` and ``"cayley"``. + rng : int or numpy.random.Generator, optional + Seed or generator used to build a reproducible starting vector when + ``v0`` is ``None``. An integer seed is the most convenient choice, + e.g. ``rng=0``. Ignored when ``v0`` is supplied. Note that + ``scipy.sparse.linalg.eigsh`` has no native seed parameter, so + this is handled by generating ``v0`` internally via + ``numpy.random.default_rng(rng)``. Returns ------- eigenvalues : np.ndarray - Array of k eigenvalues, shape (k,). For closed meshes or Neumann - boundary condition, ``0`` will be the first eigenvalue (with - constant eigenvector). + Array of k eigenvalues, shape (k,), sorted in ascending order. + For closed meshes or Neumann boundary condition, ``0`` will be + the first eigenvalue (with constant eigenvector). eigenvectors : np.ndarray Array representing the k eigenvectors, shape (n_vertices, k). The column ``eigenvectors[:, i]`` is the eigenvector corresponding @@ -691,6 +723,12 @@ def eigs(self, k: int = 10, sigma: float = -0.01) -> tuple[np.ndarray, np.ndarra """ from scipy.sparse.linalg import LinearOperator, eigsh + # eigsh has no native seed parameter; derive v0 from rng when needed. + if v0 is None and rng is not None: + v0 = np.random.default_rng(rng).standard_normal( + self.stiffness.shape[0] + ).astype(self.stiffness.dtype) + if self.use_cholmod: logger.info("Solver: Cholesky decomposition from scikit-sparse cholmod ...") chol = self.sksparse.cholmod.cholesky(self.stiffness - sigma * self.mass) @@ -711,9 +749,19 @@ def eigs(self, k: int = 10, sigma: float = -0.01) -> tuple[np.ndarray, np.ndarra dtype=self.stiffness.dtype, ) eigenvalues, eigenvectors = eigsh( - self.stiffness, k, self.mass, sigma=sigma, OPinv=op_inv + self.stiffness, + k, + self.mass, + sigma=sigma, + which=which, + v0=v0, + mode=mode, + OPinv=op_inv, ) - return eigenvalues, eigenvectors + # eigsh with shift-invert does not guarantee sorted output (especially + # for sigma values far from 0). Always sort ascending for consistency. + sort_idx = np.argsort(eigenvalues) + return eigenvalues[sort_idx], eigenvectors[:, sort_idx] def poisson( self, diff --git a/lapy/utils/tests/test_solver.py b/lapy/utils/tests/test_solver.py index 55011c6..5b0bc53 100644 --- a/lapy/utils/tests/test_solver.py +++ b/lapy/utils/tests/test_solver.py @@ -1,4 +1,4 @@ -"""Tests for Solver.poisson — backward compatibility and 2-D rhs support.""" +"""Tests for Solver.eigs and Solver.poisson — parameters and 2-D rhs support.""" import numpy as np import pytest @@ -12,6 +12,31 @@ def tria_mesh(): return TriaMesh.read_off("data/square-mesh.off") +# --------------------------------------------------------------------------- +# eigs — new parameter and sorting tests +# --------------------------------------------------------------------------- + + +def test_eigs_rng_int_reproducible(tria_mesh): + """Two calls with the same integer rng seed must return identical results.""" + fem = Solver(tria_mesh, lump=True) + evals1, evecs1 = fem.eigs(k=4, rng=42) + evals2, evecs2 = fem.eigs(k=4, rng=42) + np.testing.assert_array_equal(evals1, evals2) + np.testing.assert_array_equal(evecs1, evecs2) + + +def test_eigs_v0_takes_precedence_over_rng(tria_mesh): + """Explicit v0 must take precedence over rng.""" + fem = Solver(tria_mesh, lump=True) + v0 = np.random.default_rng(0).standard_normal(len(tria_mesh.v)) + # Same v0 with different rng seeds must give identical results. + evals1, evecs1 = fem.eigs(k=4, v0=v0, rng=0) + evals2, evecs2 = fem.eigs(k=4, v0=v0, rng=99) + np.testing.assert_array_equal(evals1, evals2) + np.testing.assert_array_equal(evecs1, evecs2) + + def test_poisson_scalar_and_1d_return_1d(tria_mesh): """Scalar and 1-D rhs must return a 1-D array (backward compatibility).""" fem = Solver(tria_mesh, lump=True) From 22caf0f0dc538bcc26cf26f73b55919f1f0f4651 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Sun, 15 Mar 2026 10:07:14 +0100 Subject: [PATCH 2/5] fix shpinx --- lapy/solver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lapy/solver.py b/lapy/solver.py index d064f87..c975f83 100644 --- a/lapy/solver.py +++ b/lapy/solver.py @@ -694,7 +694,7 @@ def eigs( the smallest Laplacian eigenvalues. Other options accepted by ``scipy.sparse.linalg.eigsh`` are ``"SM"``, ``"LA"``, ``"SA"``, and ``"BE"``. - v0 : np.ndarray, shape (n_vertices,), optional + v0 : np.ndarray of shape (n_vertices,), default=None Starting vector for the ARPACK iteration. Providing a fixed vector makes results reproducible. Takes precedence over ``rng``. If both are ``None``, ARPACK uses its own random starting vector. @@ -702,7 +702,7 @@ def eigs( Shift-invert mode passed to ``scipy.sparse.linalg.eigsh``. ``"normal"`` is the standard choice for Laplacian spectra. Other valid values are ``"buckling"`` and ``"cayley"``. - rng : int or numpy.random.Generator, optional + rng : int or numpy.random.Generator, default=None Seed or generator used to build a reproducible starting vector when ``v0`` is ``None``. An integer seed is the most convenient choice, e.g. ``rng=0``. Ignored when ``v0`` is supplied. Note that From 245f2f59ecf49d486b3464661c87a27cb4f90a45 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Sun, 15 Mar 2026 10:13:49 +0100 Subject: [PATCH 3/5] fix dashes and remove pandas intersphinx --- doc/conf.py | 1 - lapy/utils/_imports.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/doc/conf.py b/doc/conf.py index f18653b..23f8571 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -118,7 +118,6 @@ "matplotlib": ("https://matplotlib.org/stable", None), "mne": ("https://mne.tools/stable/", None), "numpy": ("https://numpy.org/doc/stable", None), - "pandas": ("https://pandas.pydata.org/pandas-docs/stable/", None), "python": ("https://docs.python.org/3", None), "scipy": ("https://docs.scipy.org/doc/scipy", None), "sklearn": ("https://scikit-learn.org/stable/", None), diff --git a/lapy/utils/_imports.py b/lapy/utils/_imports.py index 637db89..198f644 100644 --- a/lapy/utils/_imports.py +++ b/lapy/utils/_imports.py @@ -41,7 +41,7 @@ def import_optional_dependency( False. Raises - ------- + ------ ImportError dependency not found; see raise_error """ From f33264dbf6b78d10cf7ff0244f545b9ed9ebb90c Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Mon, 16 Mar 2026 09:05:45 +0100 Subject: [PATCH 4/5] better v0 and rng treatment and docstring --- lapy/solver.py | 45 +++++++++++++++++++++++++++++---------------- 1 file changed, 29 insertions(+), 16 deletions(-) diff --git a/lapy/solver.py b/lapy/solver.py index c975f83..ee7216a 100644 --- a/lapy/solver.py +++ b/lapy/solver.py @@ -696,38 +696,51 @@ def eigs( and ``"BE"``. v0 : np.ndarray of shape (n_vertices,), default=None Starting vector for the ARPACK iteration. Providing a fixed vector - makes results reproducible. Takes precedence over ``rng``. If both - are ``None``, ARPACK uses its own random starting vector. + makes results reproducible. If ``None``, ARPACK uses a random + starting vector seeded by ``rng``. mode : str, default="normal" Shift-invert mode passed to ``scipy.sparse.linalg.eigsh``. ``"normal"`` is the standard choice for Laplacian spectra. Other valid values are ``"buckling"`` and ``"cayley"``. rng : int or numpy.random.Generator, default=None - Seed or generator used to build a reproducible starting vector when - ``v0`` is ``None``. An integer seed is the most convenient choice, - e.g. ``rng=0``. Ignored when ``v0`` is supplied. Note that - ``scipy.sparse.linalg.eigsh`` has no native seed parameter, so - this is handled by generating ``v0`` internally via - ``numpy.random.default_rng(rng)``. + Seed or generator used to produce a reproducible starting vector + when ``v0`` is ``None``. An integer seed is the most convenient + choice, e.g. ``rng=0``. + On scipy >= 1.17 both ``v0`` and ``rng`` are forwarded directly to + ``scipy.sparse.linalg.eigsh``, which owns all consistency checking. + On older scipy versions the behaviour is replicated locally: + ``numpy.random.default_rng(rng)`` is used to generate a starting + vector with ``uniform(-1, 1)`` (matching scipy's own + initialisation); ``rng`` is ignored when ``v0`` is provided. Returns ------- eigenvalues : np.ndarray Array of k eigenvalues, shape (k,), sorted in ascending order. - For closed meshes or Neumann boundary condition, ``0`` will be - the first eigenvalue (with constant eigenvector). + For closed meshes or Neumann boundary condition with the default + small negative ``sigma``, ``0`` will be the first eigenvalue + (with constant eigenvector). eigenvectors : np.ndarray Array representing the k eigenvectors, shape (n_vertices, k). The column ``eigenvectors[:, i]`` is the eigenvector corresponding to ``eigenvalues[i]``. """ + import inspect + from scipy.sparse.linalg import LinearOperator, eigsh - # eigsh has no native seed parameter; derive v0 from rng when needed. - if v0 is None and rng is not None: - v0 = np.random.default_rng(rng).standard_normal( - self.stiffness.shape[0] - ).astype(self.stiffness.dtype) + if "rng" in inspect.signature(eigsh).parameters: + # scipy >= 1.17: pass v0 and rng verbatim; eigsh owns all + # consistency checking (e.g. v0 takes precedence over rng). + start_kwargs: dict = {"v0": v0, "rng": rng} + else: + # scipy < 1.17: replicate eigsh's rng behaviour locally. + # Only generate v0 when it is not already provided. + if v0 is None and rng is not None: + v0 = np.random.default_rng(rng).uniform( + -1.0, 1.0, self.stiffness.shape[0] + ).astype(self.stiffness.dtype) + start_kwargs = {"v0": v0} if self.use_cholmod: logger.info("Solver: Cholesky decomposition from scikit-sparse cholmod ...") @@ -754,9 +767,9 @@ def eigs( self.mass, sigma=sigma, which=which, - v0=v0, mode=mode, OPinv=op_inv, + **start_kwargs, ) # eigsh with shift-invert does not guarantee sorted output (especially # for sigma values far from 0). Always sort ascending for consistency. From 7fa3d377c299e90977c9942dfed315497d838bfa Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Mon, 16 Mar 2026 13:48:57 +0100 Subject: [PATCH 5/5] fix test, doc and eigsh check at module level --- lapy/solver.py | 31 +++++++++++++++++-------------- lapy/utils/tests/test_solver.py | 26 ++++++++++++++++++++++---- 2 files changed, 39 insertions(+), 18 deletions(-) diff --git a/lapy/solver.py b/lapy/solver.py index ee7216a..daa61ca 100644 --- a/lapy/solver.py +++ b/lapy/solver.py @@ -1,9 +1,11 @@ import importlib +import inspect import logging import sys import numpy as np from scipy import sparse +from scipy.sparse.linalg import eigsh as _eigsh from .tet_mesh import TetMesh from .tria_mesh import TriaMesh @@ -11,6 +13,9 @@ logger = logging.getLogger(__name__) +# Evaluated once at import time; avoids repeated introspection on every eigs() call. +_EIGSH_SUPPORTS_RNG = "rng" in inspect.signature(_eigsh).parameters + class Solver: """FEM solver for Laplace Eigenvalue and Poisson Equation. @@ -706,11 +711,11 @@ def eigs( Seed or generator used to produce a reproducible starting vector when ``v0`` is ``None``. An integer seed is the most convenient choice, e.g. ``rng=0``. - On scipy >= 1.17 both ``v0`` and ``rng`` are forwarded directly to - ``scipy.sparse.linalg.eigsh``, which owns all consistency checking. - On older scipy versions the behaviour is replicated locally: - ``numpy.random.default_rng(rng)`` is used to generate a starting - vector with ``uniform(-1, 1)`` (matching scipy's own + When ``scipy.sparse.linalg.eigsh`` exposes an ``rng`` parameter, + both ``v0`` and ``rng`` are forwarded directly to it, which owns + all consistency checking. Otherwise the behaviour is replicated + locally: ``numpy.random.default_rng(rng)`` is used to generate a + starting vector with ``uniform(-1, 1)`` (matching scipy's own initialisation); ``rng`` is ignored when ``v0`` is provided. Returns @@ -725,17 +730,15 @@ def eigs( The column ``eigenvectors[:, i]`` is the eigenvector corresponding to ``eigenvalues[i]``. """ - import inspect - - from scipy.sparse.linalg import LinearOperator, eigsh + from scipy.sparse.linalg import LinearOperator - if "rng" in inspect.signature(eigsh).parameters: - # scipy >= 1.17: pass v0 and rng verbatim; eigsh owns all - # consistency checking (e.g. v0 takes precedence over rng). + if _EIGSH_SUPPORTS_RNG: + # eigsh supports rng natively: pass v0 and rng verbatim and let + # eigsh own all consistency checking (v0 takes precedence over rng). start_kwargs: dict = {"v0": v0, "rng": rng} else: - # scipy < 1.17: replicate eigsh's rng behaviour locally. - # Only generate v0 when it is not already provided. + # eigsh has no rng parameter: replicate its behaviour locally. + # Only generate v0 from rng when no explicit v0 is provided. if v0 is None and rng is not None: v0 = np.random.default_rng(rng).uniform( -1.0, 1.0, self.stiffness.shape[0] @@ -761,7 +764,7 @@ def eigs( shape=self.stiffness.shape, dtype=self.stiffness.dtype, ) - eigenvalues, eigenvectors = eigsh( + eigenvalues, eigenvectors = _eigsh( self.stiffness, k, self.mass, diff --git a/lapy/utils/tests/test_solver.py b/lapy/utils/tests/test_solver.py index 5b0bc53..906264a 100644 --- a/lapy/utils/tests/test_solver.py +++ b/lapy/utils/tests/test_solver.py @@ -12,6 +12,24 @@ def tria_mesh(): return TriaMesh.read_off("data/square-mesh.off") +# --------------------------------------------------------------------------- +# Helper +# --------------------------------------------------------------------------- + + +def _assert_evecs_allclose(evecs1, evecs2, **kwargs): + """Assert eigenvectors agree up to a per-column sign flip. + + The eigenspace returned by ARPACK/eigsh is unique only up to sign (and, + for degenerate eigenvalues, up to rotation within the subspace). For + non-degenerate eigenvalues the per-column dot product detects the sign and + aligns evecs2 before the numerical comparison. + """ + signs = np.sign(np.einsum("ij,ij->j", evecs1, evecs2)) + signs[signs == 0] = 1 # zero dot product: columns are already orthogonal + np.testing.assert_allclose(evecs1, evecs2 * signs, **kwargs) + + # --------------------------------------------------------------------------- # eigs — new parameter and sorting tests # --------------------------------------------------------------------------- @@ -22,8 +40,8 @@ def test_eigs_rng_int_reproducible(tria_mesh): fem = Solver(tria_mesh, lump=True) evals1, evecs1 = fem.eigs(k=4, rng=42) evals2, evecs2 = fem.eigs(k=4, rng=42) - np.testing.assert_array_equal(evals1, evals2) - np.testing.assert_array_equal(evecs1, evecs2) + np.testing.assert_allclose(evals1, evals2, rtol=1e-12) + _assert_evecs_allclose(evecs1, evecs2, rtol=1e-10, atol=1e-12) def test_eigs_v0_takes_precedence_over_rng(tria_mesh): @@ -33,8 +51,8 @@ def test_eigs_v0_takes_precedence_over_rng(tria_mesh): # Same v0 with different rng seeds must give identical results. evals1, evecs1 = fem.eigs(k=4, v0=v0, rng=0) evals2, evecs2 = fem.eigs(k=4, v0=v0, rng=99) - np.testing.assert_array_equal(evals1, evals2) - np.testing.assert_array_equal(evecs1, evecs2) + np.testing.assert_allclose(evals1, evals2, rtol=1e-12) + _assert_evecs_allclose(evecs1, evecs2, rtol=1e-10, atol=1e-12) def test_poisson_scalar_and_1d_return_1d(tria_mesh):