fix: take the sparse solver path only for a positive-definite K

CHANGELOG.md

@@ -8,7 +8,29 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/).
 
 ## [Unreleased]
 
-(nothing yet)
+### Fixed
+
+- **The sparse solver path dropped negative modes.** Above 500 DOF a
+  symmetric solve went through `eigsh` shift-invert at zero, which
+  returns the eigenvalues smallest in *magnitude*. Those are the lowest
+  only when `K` is positive definite. With an indefinite `K`, such as a
+  tower loaded past buckling by `run(gravity=...)`, a strongly negative
+  mode was left out and the window backfilled with positive ones, while
+  the dense path below the threshold kept it. The answer therefore
+  changed with mesh refinement alone. `K = diag(−100, 1, …, 500)` returned
+  `[1, 2, 3]` for three modes instead of `[−100, 1, 2]`, and the bundled
+  NREL 5MW land tower at 1000 times its RNA mass returned four buckled
+  modes on the default mesh but two buckled and two positive ones once
+  refined to 120 nodes. `solve_modes` now takes the sparse path only
+  after a banded Cholesky confirms `K` is positive definite, and sends
+  anything else to the dense path, which returns the algebraically
+  lowest modes. The check costs about a millisecond of factorisation at
+  9000 DOF because the assembled stiffness has a half-bandwidth of 10.
+  The sparse path also now builds its symmetrised matrices sparsely
+  instead of densely, which more than pays for the check: a 1000-node
+  tower solves in 5.3 s instead of 8.2 s. Positive-definite results are
+  unchanged within the run-to-run scatter of the sparse solver itself
+  (about 1e-10 relative, from ARPACK's random starting vector).
 
 ## [1.19.0] — 2026-09-30
 

VALIDATION.md

@@ -242,6 +242,8 @@ is behavioural / contract-style.
 | Sparse path triggered above threshold | construction | log message announces "sparse shift-invert" | always | (within tol) | [`tests/fem/test_sparse_solver.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_sparse_solver.py) | no |
 | Sparse path skipped below threshold | construction | log message does not mention "sparse" | always | (within tol) | [`tests/fem/test_sparse_solver.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_sparse_solver.py) | no |
 | Asymmetric problems fall back to dense general `eig` | construction (asymmetric K) | log message announces "dense general"; sparse not invoked | always | (within tol) | [`tests/fem/test_sparse_solver.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_sparse_solver.py) | no |
+| Indefinite `K` keeps the algebraically lowest window on every path | construction (diagonal, coupled with scaled mass, zero modes, SPD control) + NREL 5MW land tower under gravity at 1000× RNA mass, default vs 120-node mesh | lowest-k eigenvalues; negative and zero modes kept; SPD still sparse; buckled tower identical on both meshes | `rtol = 1e-9` | exact | [`tests/fem/test_sparse_solver.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_sparse_solver.py) | no |
+| Positive-definiteness guard (banded Cholesky) agrees with a dense Cholesky | construction (random symmetric, tridiagonal, singular, full-width coupling, tiny asymmetry) | definite / not definite | always | (within tol) | [`tests/fem/test_sparse_solver.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_sparse_solver.py) | no |
 | `pybmodes patch --dry-run` writes nothing (mtime check) | construction | mtime of every staged file unchanged | exact match | (within tol) | [`tests/test_validate.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_validate.py) | yes |
 | `pybmodes patch --diff` PR format contains `×` improvement ratio | construction | stdout contains `×`, `RMS improvement:`, `better` literals | exact match | (within tol) | [`tests/test_validate.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_validate.py) | yes |
 | `pybmodes patch --output` leaves source byte-identical | construction | source-file sha256 unchanged | exact match | (within tol) | [`tests/test_validate.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_validate.py) | yes |

src/pybmodes/fem/solver.py

@@ -26,7 +26,12 @@
    ``ngd``). Routes through ``scipy.sparse.linalg.eigsh`` with
    ``sigma=0`` shift-invert; an order-of-magnitude faster than the
    dense LAPACK solve for the few-lowest-modes case on a 500+ DOF
-   tower mesh.
+   tower mesh. Taken only when ``K`` is positive definite, checked by
+   a banded Cholesky (:func:`_is_positive_definite`): shift-invert at
+   zero returns the eigenvalues *smallest in magnitude*, which are the
+   lowest only when none is negative. An indefinite ``K`` (a column
+   loaded past buckling) goes to the dense path, which returns the
+   algebraically lowest modes.
 2. **Dense symmetric** — ``scipy.linalg.eigh`` on the symmetrised
    matrices, applied to the *inverted* pencil ``M x = μ (K + s M) x``
    so that the lowest modes are resolved to relative rather than
@@ -138,13 +143,16 @@
 import logging
 import warnings
 from dataclasses import dataclass
-from typing import Literal, overload
+from typing import TYPE_CHECKING, Literal, overload
 
 import numpy as np
 from scipy.linalg import eig, eigh
 
 from pybmodes.options import DEFAULT_SOLVER_OPTIONS as _SOLVER_OPTIONS
 
+if TYPE_CHECKING:
+    from scipy.sparse import csc_matrix
+
 _log = logging.getLogger(__name__)
 
 # Above this reduced-system size the dense conditioning estimate
@@ -296,15 +304,30 @@ def solve_modes(
     eigvals: np.ndarray | None = None
     eigvecs: np.ndarray | None = None
 
-    # Sparse path — symmetric, big enough, small-subset request.
-    if (
+    # Sparse path — symmetric, big enough, small-subset request, and a
+    # positive-definite K. Shift-invert at zero returns the eigenvalues
+    # nearest zero in magnitude; with a negative one among them (a column
+    # past buckling) that is not the lowest window, and the unstable mode
+    # would be dropped silently. The dense path keeps it.
+    sparse_candidate = (
         sym
         and ngd > _SPARSE_NDOF_THRESHOLD
         and n_modes is not None
         and n_modes < ngd // 2
-    ):
+    )
+    k_sp = _symmetric_sparse(gk) if sparse_candidate else None
+    if k_sp is not None and not _is_positive_definite(k_sp):
+        k_sp = None
+        _log.info(
+            "solve_modes: K is not positive definite; using the dense "
+            "path so negative and zero modes stay in the window (ngd=%d)",
+            ngd,
+        )
+    if k_sp is not None and n_modes is not None:
         try:
-            eigvals, eigvecs = _solve_sparse_shift_invert(gk, gm, n_modes)
+            eigvals, eigvecs = _solve_sparse_shift_invert(
+                gk, gm, n_modes, k_sp=k_sp,
+            )
             path = "sparse_shift_invert"
             _log.info(
                 "solve_modes: sparse shift-invert path "
@@ -791,8 +814,47 @@ def _mass_matrix_cond(
 # Path implementations
 # ---------------------------------------------------------------------------
 
+def _symmetric_sparse(a: np.ndarray) -> csc_matrix:
+    """``0.5 (A + A^T)`` as a CSC matrix, without another dense copy.
+
+    Bit-identical to symmetrising densely first: the same pairwise sum,
+    and the halving is exact.
+    """
+    from scipy.sparse import csc_matrix
+
+    s = csc_matrix(a)
+    return csc_matrix((s + s.T) * 0.5)
+
+
+def _is_positive_definite(k_sp: csc_matrix) -> bool:
+    """Whether a symmetric sparse ``K`` admits a Cholesky factorisation.
+
+    The assembled stiffness is banded (half-bandwidth 10 on a tower mesh,
+    whatever its size), so the check runs on the band: ``O(n b²)``, about
+    a millisecond at 9000 DOF where a dense factorisation takes seconds.
+    The band is measured from ``K`` itself, so a model with wider coupling
+    is still checked exactly, only less cheaply.
+    """
+    from scipy.linalg import cholesky_banded
+    from scipy.sparse import triu
+
+    upper = triu(k_sp).tocoo()
+    if upper.nnz == 0:
+        return False
+    offset = upper.col - upper.row
+    b = int(offset.max())
+    ab = np.zeros((b + 1, k_sp.shape[0]))
+    np.add.at(ab, (b - offset, upper.col), upper.data)
+    try:
+        cholesky_banded(ab, lower=False, check_finite=False)
+    except np.linalg.LinAlgError:
+        return False
+    return True
+
+
 def _solve_sparse_shift_invert(
     gk: np.ndarray, gm: np.ndarray, n_modes: int,
+    *, k_sp: csc_matrix | None = None,
 ) -> tuple[np.ndarray, np.ndarray]:
     """Sparse symmetric generalised eigensolve via shift-invert near zero.
 
@@ -802,15 +864,17 @@ def _solve_sparse_shift_invert(
     ``OP = (K - σM)^-1 · M = K^-1 M`` under ``mode='normal'``. The
     eigenvalues of OP are ``1/λ``; ``which='LM'`` returns the largest,
     i.e. the smallest ``λ`` — exactly the modal-analysis ask.
+
+    Smallest in *magnitude*: that is the lowest window only when ``K`` is
+    positive definite, which :func:`solve_modes` checks before calling
+    this. ``k_sp`` passes in the symmetrised sparse ``K`` that check
+    already built.
     """
-    from scipy.sparse import csc_matrix
     from scipy.sparse.linalg import eigsh
 
     # Symmetrise to suppress sub-ULP scatter before factorisation.
-    gk_sym = 0.5 * (gk + gk.T)
-    gm_sym = 0.5 * (gm + gm.T)
-    K_sp = csc_matrix(gk_sym)
-    M_sp = csc_matrix(gm_sym)
+    K_sp = k_sp if k_sp is not None else _symmetric_sparse(gk)
+    M_sp = _symmetric_sparse(gm)
 
     eigvals, eigvecs = eigsh(
         K_sp,

tests/fem/test_sparse_solver.py

@@ -1,10 +1,16 @@
 """Tests for the sparse-solver dispatch in ``pybmodes.fem.solver``.
 
-The solver gates the sparse shift-invert path on three conditions:
+The solver gates the sparse shift-invert path on four conditions:
 
   ngd > _SPARSE_NDOF_THRESHOLD AND symmetric matrices AND a small
   subset of modes was requested (``n_modes is not None`` and
-  ``n_modes < ngd // 2``).
+  ``n_modes < ngd // 2``) AND a positive-definite K.
+
+The last one exists because shift-invert at zero returns the
+eigenvalues smallest in magnitude, which are the lowest only when none
+is negative. Without it an indefinite K (a column loaded past buckling)
+lost its unstable modes on the sparse path while the dense path kept
+them, so the answer changed with mesh refinement alone.
 
 These tests build small synthetic symmetric / asymmetric problems
 (no FEM pipeline needed) so they're fast and self-contained.
@@ -144,3 +150,142 @@ def test_dense_fallback_on_asymmetric(caplog: pytest.LogCaptureFixture) -> None:
     assert eigvals.size == 6
     assert eigvecs.shape == (n, 6)
     assert np.all(np.isfinite(eigvals))
+
+
+# ---------------------------------------------------------------------------
+# Indefinite K: the sparse path must not change the returned window
+# ---------------------------------------------------------------------------
+
+_N_BIG = _SPARSE_NDOF_THRESHOLD + 1
+
+
+def _rotated(d: np.ndarray, seed: int) -> np.ndarray:
+    """``Q diag(d) Q^T`` with a random orthogonal ``Q``: a fully coupled
+    symmetric matrix with a chosen spectrum."""
+    rng = np.random.default_rng(seed)
+    q, _r = np.linalg.qr(rng.standard_normal((d.size, d.size)))
+    k = q @ np.diag(d) @ q.T
+    return 0.5 * (k + k.T)
+
+
+def _lowest(gk: np.ndarray, gm: np.ndarray, n: int) -> np.ndarray:
+    from scipy.linalg import eigh
+
+    return eigh(gk, gm, eigvals_only=True, subset_by_index=(0, n - 1))
+
+
+class TestIndefiniteStiffnessKeepsTheLowestWindow:
+    def test_large_negative_mode_is_not_dropped(self):
+        """The reported case: shift-invert at zero picked [1, 2, 3] here
+        because |-100| is larger than 3."""
+        gk = np.diag(np.concatenate([[-100.0], np.arange(1.0, _N_BIG)]))
+        vals, _vecs, diag = solve_modes(
+            gk, np.eye(_N_BIG), 3, return_diagnostics=True,
+        )
+        np.testing.assert_allclose(vals, [-100.0, 1.0, 2.0], rtol=1e-12)
+        assert diag.path == "dense_symmetric"
+        assert diag.sparse_fallback is False
+
+    def test_coupled_stiffness_with_scaled_mass(self):
+        rng = np.random.default_rng(4)
+        d = np.concatenate([[-50.0, -7.0], rng.uniform(1.0, 400.0, _N_BIG - 2)])
+        gk = _rotated(d, seed=5)
+        gm = np.diag(rng.uniform(0.5, 20.0, _N_BIG))
+        vals, vecs = solve_modes(gk, gm, 4)
+        np.testing.assert_allclose(vals, _lowest(gk, gm, 4), rtol=1e-9)
+        assert vals[0] < 0.0 and vals[1] < 0.0
+        residual = gk @ vecs - (gm @ vecs) * vals
+        assert np.max(np.abs(residual)) < 1e-8 * np.max(np.abs(gk))
+
+    def test_zero_modes_stay_in_the_window(self):
+        d = np.concatenate([np.zeros(3), np.arange(1.0, _N_BIG - 2)])
+        gk = _rotated(d, seed=6)
+        vals, _vecs, diag = solve_modes(
+            gk, np.eye(_N_BIG), 5, return_diagnostics=True,
+        )
+        assert diag.path == "dense_symmetric"
+        assert np.all(np.abs(vals[:3]) < 1e-9 * float(np.max(d)))
+        np.testing.assert_allclose(vals[3:], [1.0, 2.0], rtol=1e-9)
+
+    def test_positive_definite_stiffness_still_goes_sparse(self):
+        gk, gm = _spd_problem(_N_BIG, seed=7)
+        vals, _vecs, diag = solve_modes(gk, gm, 4, return_diagnostics=True)
+        assert diag.path == "sparse_shift_invert"
+        np.testing.assert_allclose(vals, _lowest(gk, gm, 4), rtol=1e-8)
+
+
+class TestPositiveDefiniteCheck:
+    @staticmethod
+    def _check(a: np.ndarray) -> bool:
+        return solver._is_positive_definite(solver._symmetric_sparse(a))
+
+    def test_matches_a_dense_cholesky(self):
+        from scipy.linalg import cho_factor
+
+        rng = np.random.default_rng(8)
+        for shift in (-5.0, -0.5, 0.5, 5.0):
+            a = rng.standard_normal((60, 60))
+            k = 0.5 * (a + a.T) + shift * 10.0 * np.eye(60)
+            try:
+                cho_factor(k)
+                expected = True
+            except np.linalg.LinAlgError:
+                expected = False
+            assert self._check(k) is expected, shift
+
+    def test_banded_tridiagonal(self):
+        n = 200
+        k = 2.0 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
+        assert self._check(k)          # second-difference matrix, SPD
+        k[0, 0] = -1.0
+        assert not self._check(k)
+
+    def test_singular_and_empty_are_not_definite(self):
+        assert not self._check(np.diag([1.0, 0.0, 2.0]))
+        assert not self._check(np.zeros((4, 4)))
+
+    def test_wide_coupling_is_still_exact(self):
+        """A far off-diagonal coupling widens the band to the whole
+        matrix; the check must still see it."""
+        n = 40
+        k = np.eye(n)
+        k[0, n - 1] = k[n - 1, 0] = 2.0   # makes the 2x2 minor indefinite
+        assert not self._check(k)
+        k[0, n - 1] = k[n - 1, 0] = 0.5
+        assert self._check(k)
+
+    def test_symmetrises_tiny_asymmetry(self):
+        k = np.diag([3.0, 2.0, 1.0])
+        k[0, 1], k[1, 0] = 1.0, 1.0 + 1e-14
+        assert self._check(k)
+
+
+def test_buckled_tower_does_not_depend_on_the_mesh():
+    """A real model: the bundled NREL 5MW land tower under gravity with
+    1000 times its RNA mass is far past buckling. On the default mesh
+    (dense path) its lowest four eigenvalues are all negative, reported
+    as 0 Hz. Refined past the sparse threshold it used to return two
+    positive modes among them instead."""
+    import warnings
+
+    from pybmodes.cli import _resolve_examples_root
+    from pybmodes.models import Tower
+
+    bmi = (
+        _resolve_examples_root() / "sample_inputs" / "reference_turbines"
+        / "01_nrel5mw_land" / "01_nrel5mw_land_tower.bmi"
+    )
+
+    def run(nodes):
+        tower = Tower(bmi)
+        if nodes is not None:
+            tower.refine_mesh(nodes)
+        tower._bmi.tip_mass.mass *= 1000.0
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            return tower.run(n_modes=4, check_model=False, gravity=True)
+
+    coarse, fine = run(None), run(120)
+    assert fine.diagnostics.path == "dense_symmetric"
+    np.testing.assert_array_equal(coarse.frequencies, 0.0)
+    np.testing.assert_array_equal(fine.frequencies, 0.0)