From fb4b2988ba333a5ab5536b2e0e102e648594be21 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 21 Oct 2026 16:35:05 +0900 Subject: [PATCH] fix: take the sparse solver path only for a positive-definite K Shift-invert at zero returns the eigenvalues smallest in magnitude, which are the lowest only when K is positive definite. With an indefinite K, such as a tower past buckling under gravity, the sparse path dropped the strongly negative modes that the dense path keeps, so the result changed with mesh refinement alone. A banded Cholesky now gates the sparse path; anything not definite goes dense. The sparse path also builds its symmetrised matrices sparsely, which more than pays for the check. --- CHANGELOG.md | 24 ++++- VALIDATION.md | 2 + src/pybmodes/fem/solver.py | 86 +++++++++++++++--- tests/fem/test_sparse_solver.py | 149 +++++++++++++++++++++++++++++++- 4 files changed, 247 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 691b37a..db96b44 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/VALIDATION.md b/VALIDATION.md index 2f10b11..4a0589a 100644 --- a/VALIDATION.md +++ b/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 | diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 0815a7c..ff540d1 100644 --- a/src/pybmodes/fem/solver.py +++ b/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, diff --git a/tests/fem/test_sparse_solver.py b/tests/fem/test_sparse_solver.py index 09ab849..69238c8 100644 --- a/tests/fem/test_sparse_solver.py +++ b/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)