From 0d290011446d4384cf9067b60df7b81f15052b33 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 17:09:17 +0900 Subject: [PATCH 01/28] fix: catch the silent symmetric-eigensolver breakdown on a near-singular M Both scipy.linalg.eigh and scipy.sparse.linalg.eigsh reduce K x = lam M x through a Cholesky factor of the mass matrix, and that reduction loses accuracy when a very light beam carries a very heavy lump. LAPACK does not raise there, it returns wrong frequencies. A 100 m cantilever with a 4000:1 lump-to-beam ratio reported 0.103 Hz against a true 0.0436 Hz, and the answer wandered non-monotonically with mesh density. solve_modes now measures the backward error of every symmetric solve and redoes it on the general dense path, which factorises neither matrix. Two conditions gate the swap and both are load-bearing. The general result is taken only on a decisive win. A real deck can carry a large backward error without being broken: the bundled NREL 5MW land tower sits at ~2e-2 because its adapter leaves M at cond ~4e10, and there the general path is 1.4x better while splitting a degenerate fore-aft / side-side pair the symmetric solver resolves exactly, which the FA/SS classifier depends on. A genuine breakdown improves by nine orders. Rigid-body modes are excluded from the measurement. For those K x and lam are both ~0, so the relative residual is a ratio of two near-zero quantities and reads ~1 however exact the eigenpair is. Judging a free-free floating solve by the raw maximum condemned it, and acting on that verdict was worse than doing nothing, since the general path filters non-positive eigenvalues and would delete the zero-frequency mode. No existing result changes: the full suite passes with no test edited except one whose workaround for this defect is no longer needed. --- CHANGELOG.md | 33 ++- VALIDATION.md | 2 + src/pybmodes/fem/solver.py | 118 +++++++++- src/pybmodes/options.py | 26 +++ tests/fem/test_gravity_and_point_mass.py | 13 +- tests/fem/test_ill_conditioned_mass.py | 274 +++++++++++++++++++++++ 6 files changed, 458 insertions(+), 8 deletions(-) create mode 100644 tests/fem/test_ill_conditioned_mass.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 4687948..ee2f4f2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,38 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ## [Unreleased] -(nothing yet) +### Fixed + +- **The symmetric eigensolvers could return confidently wrong low modes + on a near-singular mass matrix, silently.** Both `scipy.linalg.eigh` + and `scipy.sparse.linalg.eigsh` reduce `K x = λ M x` through a Cholesky + factor of the mass matrix, and that reduction loses accuracy when a + very light beam carries a very heavy lump. LAPACK does not raise there + — it returns wrong frequencies. On a 100 m cantilever with a 4000:1 + lump-to-beam mass ratio the reported fundamental was 0.103 Hz against a + true 0.0436 Hz, a factor of 2.4, and the answer wandered + non-monotonically with mesh density. + + `solve_modes` now checks the backward error `||K x - λ M x|| / ||K x||` + of every symmetric solve and, when it is large, redoes it through the + general dense path, which factorises neither matrix. The retried result + is taken only when it is better by an order of magnitude, and a + `RuntimeWarning` names the swap. `SolverDiagnostics` gains + `residual_fallback` recording it. + + **No existing result changes.** The decisive-improvement condition is + what guarantees that: a real deck can carry a large backward error + without being broken, and on the bundled NREL 5MW land tower (whose + adapter leaves the mass matrix at cond ~4e10) the general path is only + 1.4× better while *splitting* a degenerate fore-aft / side-side pair + the symmetric solver resolves exactly. Rigid-body modes are excluded + from the check for the same reason — their relative residual is a ratio + of two near-zero quantities and is ~1 however exact the eigenpair is, + so judging a floating solve by the raw maximum would have condemned it + and then deleted its zero-frequency mode. + +- `SolverOptions` gains `residual_retry_threshold` and + `residual_retry_improvement` for the two conditions above. ## [1.18.0] — 2026-08-12 diff --git a/VALIDATION.md b/VALIDATION.md index 2e45755..1f18530 100644 --- a/VALIDATION.md +++ b/VALIDATION.md @@ -102,6 +102,8 @@ metrics: | Discrete mid-span point mass (issue #35) | cantilever carrying one lump at station $a$: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1) | 1st frequency; mesh-position independence | < 0.5 %; coarse-vs-fine < 0.2 % | (within tol) | [`tests/fem/test_gravity_and_point_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_gravity_and_point_mass.py) | no | | Distributed Winkler soil bed vs the lumped mudline condensation (issue #118) | Psaroudakis et al. (2021) / Yu & Amdahl (2023) Eq. 25 is the exact static condensation of a constant-EI pile on a bed of rate $k = D_P E_{SO}$ | coupled 1st frequency, distributed bed vs lumped springs | < 1 % | 0.4 % (the embedded pile inertia the condensed form drops) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | | Distributed Winkler bed converges on the rigid mudline clamp (issue #118) | limit $E_{SO} \to \infty$; residual compliance scales as the elastic length $(4EI/k)^{1/4}$ | 1st frequency at $E_{SO} \times 10^8$ vs the clamped model | < 1 % | (within tol) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | +| Near-singular mass matrix does not silently corrupt the low modes | cantilever carrying one lump: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1), at a 4e5:1 lump-to-beam mass ratio | 1st frequency, dense and sparse dispatch sizes | < 0.5 % | (within tol; the unguarded symmetric solve is 137 % out) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no | +| Rigid-body modes are not mistaken for a solver breakdown | construction (rank-deficient $K$, well-conditioned $M$) | zero mode retained; no retry triggered | exact | (within tol) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no | | Deck-reader `n_nodes` refinement (issue #58) | Euler-Bernoulli closed form, modes 1-3, plus self-convergence $n{=}100$ vs $200$ | bending frequencies | < 1 %; self-convergence < 0.2 % | (within tol) | [`tests/test_refine_mesh.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_refine_mesh.py) | no | **Citations** (full author / year forms used in the table above). diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 3d04605..00733fd 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -32,7 +32,33 @@ when the sparse path fails to converge (logged as a warning). 3. **Dense general** — ``scipy.linalg.eig`` for genuinely asymmetric systems (offshore decks where the rigid-arm transformation makes - the platform-support block non-symmetric). Matches BModes JJ. + the platform-support block non-symmetric). Matches BModes JJ. Also + the retry path when a symmetric solve comes back with a large + backward error — see below. + +Both symmetric paths reduce ``K x = λ M x`` through a Cholesky factor +of the mass matrix, and that reduction degrades once ``M`` is nearly +singular, which a very light beam carrying a very heavy lump produces. +The failure mode is silent: LAPACK returns confidently wrong low modes +rather than raising. :func:`solve_modes` therefore checks the backward +error of the **elastic** modes of every symmetric solve (a rigid-body +mode has ``K x ≈ 0`` and ``λ ≈ 0``, so its relative residual is a ratio +of two near-zero quantities and carries no information — see +:func:`_max_elastic_residual`) and, when it exceeds +:attr:`~pybmodes.options.SolverOptions.residual_retry_threshold`, tries +the general path as well — taking its result only if it is better by +:attr:`~pybmodes.options.SolverOptions.residual_retry_improvement`, and +warning when it does. + +That second condition is the load-bearing one. A real deck can sit above +the threshold without being broken (the bundled NREL 5MW land tower +reaches ~2e-2, its adapter leaving ``M`` at cond ~4e10), and there the +general path is only marginally better while *splitting* the degenerate +fore-aft / side-side pair the symmetric solver resolves exactly — which +the FA / SS classifier downstream depends on. A true breakdown is not +marginal: it improves by nine orders of magnitude. Demanding a decisive +win keeps every validated frequency untouched and still catches the +failure this guard exists for. Note on the user-spec mode choice: ``eigsh(..., sigma=0, mode='buckling')`` reduces to ``OP = K^-1 K = I`` for ``sigma=0``, @@ -94,6 +120,14 @@ class SolverDiagnostics: attempted and failed, so the result came from the dense fallback. fallback_reason : the repr of the exception that triggered the fallback, or ``None`` when no fallback happened. + residual_fallback : ``True`` when a symmetric path returned modes whose + backward error exceeded + :attr:`~pybmodes.options.SolverOptions.residual_retry_threshold` + and the result was redone through the general dense path. The + symmetric routines factorise the mass matrix, which a very light + beam carrying a very heavy lump makes nearly singular; there they + return wrong low modes rather than failing, so the residual is + what catches it. max_residual : the largest per-mode relative residual ``||K x - λ M x|| / ||K x||`` over the returned modes (``0.0`` when no modes were returned). A healthy modal solve sits near @@ -114,6 +148,7 @@ class SolverDiagnostics: max_residual: float residuals: tuple[float, ...] matrix_cond: float | None + residual_fallback: bool = False # Sparse path activates once the reduced system has more than this # many DOFs and the caller asked for a small subset of modes. Below @@ -223,6 +258,47 @@ def solve_modes( _normalize_columns_l2(eigvecs) + # Accuracy guarantee for the symmetric paths. Both ``eigh`` and + # ``eigsh`` reduce ``K x = λ M x`` through a Cholesky factor of one of + # the matrices, and that reduction degrades once the factored matrix + # is nearly singular — a very light beam carrying a very heavy lump + # does exactly that to ``M``. The failure is silent: LAPACK returns + # confidently wrong low modes rather than raising. The backward error + # catches it (healthy solves sit at ~1e-4 or below, degraded ones + # above 1), and the general path, which factorises neither matrix, + # stays exact there. + residual_fallback = False + if sym: + worst = _max_elastic_residual(gk, gm, eigvals, eigvecs) + if worst > _SOLVER_OPTIONS.residual_retry_threshold: + alt_vals, alt_vecs = _solve_dense_general(gk, gm, n_modes) + _normalize_columns_l2(alt_vecs) + alt_worst = _max_elastic_residual(gk, gm, alt_vals, alt_vecs) + # Only take the general result on a decisive win. A marginal + # one is not a breakdown, and switching for it would churn + # validated frequencies and break the degenerate fore-aft / + # side-side pairs the symmetric solver resolves exactly — the + # bundled NREL 5MW land deck does exactly that. A genuine + # breakdown improves by many orders, not by a factor. + if alt_worst < _SOLVER_OPTIONS.residual_retry_improvement * worst: + warnings.warn( + f"the symmetric eigensolver returned modes with a " + f"backward error of {worst:.2e}, so the eigenpairs do " + f"not satisfy K x = lambda M x. Its Cholesky reduction " + f"of the mass matrix loses accuracy when that matrix is " + f"nearly singular, which a very light beam carrying a " + f"very heavy lump produces. Redone through the general " + f"dense path, which factorises neither matrix " + f"(backward error {alt_worst:.2e}); the returned modes " + f"come from that solve. Worth checking the mass " + f"distribution is the one you intended.", + RuntimeWarning, + stacklevel=2, + ) + eigvals, eigvecs = alt_vals, alt_vecs + path = "dense_general" + residual_fallback = True + # Mode-count guarantee: the general path filters complex / non- # positive eigenvalues, so it can return fewer modes than requested. # Surface that rather than letting it pass silently (a downstream @@ -258,6 +334,7 @@ def solve_modes( gk, gm, eigvals, eigvecs, path=path, symmetric=sym, n_requested=n_modes, sparse_fallback=sparse_fallback, fallback_reason=fallback_reason, + residual_fallback=residual_fallback, ) return eigvals, eigvecs, diagnostics @@ -273,6 +350,7 @@ def _build_diagnostics( n_requested: int | None, sparse_fallback: bool, fallback_reason: str | None, + residual_fallback: bool = False, ) -> SolverDiagnostics: """Assemble a :class:`SolverDiagnostics` for a completed solve.""" residuals = _modal_residuals(gk, gm, eigvals, eigvecs) @@ -287,9 +365,47 @@ def _build_diagnostics( max_residual=float(residuals.max()) if residuals.size else 0.0, residuals=tuple(float(r) for r in residuals), matrix_cond=cond, + residual_fallback=residual_fallback, ) +# A mode whose eigenvalue is below this fraction of the largest returned +# one is a rigid-body mode: a free-free floating platform has up to six, +# and an unrestrained DOF (a symmetric column's yaw) gives an exactly +# zero one. +_RIGID_BODY_EIGVAL_RTOL = 1.0e-8 + + +def _max_elastic_residual( + gk: np.ndarray, gm: np.ndarray, eigvals: np.ndarray, eigvecs: np.ndarray, +) -> float: + """Largest backward error over the **elastic** modes only. + + The relative residual ``||K x - λ M x|| / ||K x||`` is undefined for a + rigid-body mode: there ``K x ≈ 0`` and ``λ ≈ 0``, so it is a ratio of + two near-zero quantities and evaluates to ≈ 1 no matter how exact the + eigenpair is. A free-free floating model legitimately has up to six of + them, so judging a solve by the raw maximum would condemn every + floating result — and the general path drops zero eigenvalues + entirely, so acting on that verdict would delete a physically real + mode rather than improve anything. + + Rigid-body modes are therefore excluded before taking the maximum. + Returns ``0.0`` when every returned mode is rigid-body, i.e. when + there is nothing the metric can speak to. + """ + if eigvals.size == 0: + return 0.0 + scale = float(np.max(np.abs(eigvals))) + if scale <= 0.0: + return 0.0 + elastic = np.abs(eigvals) > _RIGID_BODY_EIGVAL_RTOL * scale + if not elastic.any(): + return 0.0 + r = _modal_residuals(gk, gm, eigvals, eigvecs) + return float(r[elastic].max()) + + def _modal_residuals( gk: np.ndarray, gm: np.ndarray, eigvals: np.ndarray, eigvecs: np.ndarray, ) -> np.ndarray: diff --git a/src/pybmodes/options.py b/src/pybmodes/options.py index 075d0e4..ea981ef 100644 --- a/src/pybmodes/options.py +++ b/src/pybmodes/options.py @@ -70,10 +70,36 @@ class SolverOptions: :func:`scipy.linalg.eig` instead of the symmetric :func:`scipy.linalg.eigh`. The OC3 Hywind cross-coupled ``hydro_K + mooring_K`` exercises this branch. + residual_retry_threshold : float, default 0.1 + Largest per-mode relative residual + ``||K x - λ M x|| / ||K x||`` a symmetric solve may return before + the general dense path is tried as well. ``scipy.linalg.eigh`` + reduces the generalised problem through a Cholesky factor of the + **mass** matrix, which loses accuracy once that matrix is nearly + singular — a very light beam carrying a very heavy lump — and + returns confidently wrong low modes rather than failing. + + The band this sits in is narrower than it looks. Ordinary solves + land at or below ~1e-3; a real deck whose adapter leaves ``M`` + genuinely ill-conditioned (the bundled NREL 5MW land tower, cond + ~4e10) reaches ~2e-2 and is *not* meant to trigger; the degraded + regime starts around 0.7. The default splits the last two gaps + with roughly 5x either side. + residual_retry_improvement : float, default 0.1 + How much better the general path's backward error must be before + its result is taken. The second, and more important, guard: on + that same land deck the general path is only ~1.4x better while + *breaking* a physically real degenerate fore-aft / side-side pair + the symmetric solver resolves exactly, which the downstream FA / + SS classifier depends on. A genuine breakdown is not marginal — + it improves by nine orders — so requiring a decisive win keeps + validated results untouched and still catches the real failure. """ sparse_ndof_threshold: int = 500 symmetry_rtol: float = 1.0e-12 + residual_retry_threshold: float = 0.1 + residual_retry_improvement: float = 0.1 @dataclass(frozen=True) diff --git a/tests/fem/test_gravity_and_point_mass.py b/tests/fem/test_gravity_and_point_mass.py index 20aae85..4f4f5fe 100644 --- a/tests/fem/test_gravity_and_point_mass.py +++ b/tests/fem/test_gravity_and_point_mass.py @@ -277,13 +277,14 @@ def test_lump_matches_static_cantilever_stiffness(self, tmp_path): a = 60.0 m_lump = 5.0e5 # A beam mass this far below the lump makes the generalised mass - # matrix badly conditioned, which the dense LAPACK subset path - # does not handle well; 101 elements puts the solve on the sparse - # shift-invert path, where it is clean. Nothing about the lump - # placement needs the fine mesh. - tower = _synthetic_tower(tmp_path, mass_den=1.0e-2, n_elements=101) + # matrix nearly singular, which the symmetric eigensolvers handle + # badly. The solver detects that from the backward error and + # redoes the solve on the general path, so a coarse mesh is fine + # here — see tests/fem/test_ill_conditioned_mass.py. + tower = _synthetic_tower(tmp_path, mass_den=1.0e-2, n_elements=25) tower.add_point_mass(a, m_lump) - f = tower.run(6, check_model=False).frequencies[0] + with pytest.warns(RuntimeWarning, match="backward error"): + f = tower.run(6, check_model=False).frequencies[0] f_ref = np.sqrt(3.0 * EI_PHYS / (m_lump * a**3)) / (2.0 * np.pi) assert f == pytest.approx(f_ref, rel=5.0e-3) diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py new file mode 100644 index 0000000..cef416e --- /dev/null +++ b/tests/fem/test_ill_conditioned_mass.py @@ -0,0 +1,274 @@ +"""The symmetric eigensolvers degrade silently on a near-singular mass +matrix, and the solver has to notice. + +Both ``scipy.linalg.eigh`` and ``scipy.sparse.linalg.eigsh`` reduce +``K x = lambda M x`` through a Cholesky factor of the mass matrix. When +that matrix is nearly singular — a very light beam carrying a very heavy +lump — the reduction loses accuracy, and LAPACK returns confidently +wrong low modes rather than raising. On the case pinned below the dense +symmetric path reported 0.103 Hz against a true 0.0436 Hz, a factor of +2.4, with no error and no warning. + +The guard is the backward error ``||K x - lambda M x|| / ||K x||``. It +has two conditions, and the second matters more than the first: the +error must exceed the retry threshold, *and* the general path must beat +it decisively. Being above the threshold alone is not evidence of a +breakdown — the bundled NREL 5MW land deck sits at ~2e-2 because its +adapter leaves ``M`` at cond ~4e10, and there the general path is only +1.4x better while splitting a degenerate fore-aft / side-side pair the +symmetric solver resolves exactly. Swapping for that would churn a +validated frequency by 0.84 % and break the FA / SS classifier +downstream. A real breakdown improves by nine orders, not by a factor. + +Analytical reference: a cantilever whose beam mass is negligible next to +a tip lump behaves as a spring-mass oscillator on the static tip +stiffness ``3 EI / L^3``, giving ``f = sqrt(3 EI / (m L^3)) / 2 pi`` +(Blevins 1979, Table 8-1). +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from pybmodes.fem.assembly import assemble +from pybmodes.fem.nondim import RM, ROMG, make_params, nondim_tip_mass +from pybmodes.fem.solver import eigvals_to_hz, solve_modes +from pybmodes.io.bmi import TipMassProps + +L = 100.0 +EI = 1.0e10 +M_TIP = 4.0e5 + +# The lump dominates the beam by ~4e5:1 here, which drives cond(M) past +# 1e12 and is what breaks the Cholesky reduction. +LIGHT = 1.0e-2 +REALISTIC = 1.0e3 + + +def _analytic() -> float: + return float(np.sqrt(3.0 * EI / (M_TIP * L**3)) / (2.0 * np.pi)) + + +def _cantilever_with_tip_lump( + nselt: int, mass_den: float, +) -> tuple[np.ndarray, np.ndarray]: + nd = make_params(radius=L, hub_rad=0.0, rot_rpm=0.0) + eiy = EI / nd.ref4 + eli = 1.0 / nselt + el = np.full(nselt, eli) + xb = np.array([1.0 - (i + 1) * eli for i in range(nselt)]) + tip = nondim_tip_mass( + TipMassProps(mass=M_TIP, cm_offset=0.0, cm_axial=0.0, ixx=0.0, + iyy=0.0, izz=0.0, ixy=0.0, izx=0.0, iyz=0.0), + nd, beam_type=2, id_form=1, hub_conn=1, + ) + gk, gm, _ = assemble( + nselt=nselt, el=el, xb=xb, cfe=np.zeros(nselt), + eiy=np.full(nselt, eiy), eiz=np.full(nselt, eiy), + gj=np.full(nselt, 1.0e3 * eiy), eac=np.full(nselt, 100.0), + rmas=np.full(nselt, mass_den / RM), + skm1=np.full(nselt, 1.0e-5), skm2=np.full(nselt, 1.0e-5), + eg=np.zeros(nselt), ea=np.zeros(nselt), omega2=0.0, + sec_loc=np.array([0.0, 1.0]), str_tw=np.zeros(2), hub_conn=1, + tip_mass=tip, + ) + return gk, gm + + +def _first_frequency(nselt: int, mass_den: float): + gk, gm = _cantilever_with_tip_lump(nselt, mass_den) + with pytest.warns(RuntimeWarning, match="backward error"): + eigvals, _vecs, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + return float(eigvals_to_hz(eigvals, ROMG)[0]), diag + + +class TestIllConditionedMassIsCaught: + """Sizes chosen to straddle the dense / sparse dispatch threshold.""" + + @pytest.mark.parametrize("nselt", [13, 27, 53]) + def test_dense_path_recovers_the_analytic_frequency(self, nselt): + f, diag = _first_frequency(nselt, LIGHT) + assert f == pytest.approx(_analytic(), rel=5.0e-3) + assert diag.residual_fallback is True + assert diag.path == "dense_general" + + @pytest.mark.parametrize("nselt", [13, 27, 53]) + def test_the_symmetric_result_would_have_been_wrong(self, nselt): + """Without the guard the answer is not merely imprecise.""" + from scipy.linalg import eigh + + gk, gm = _cantilever_with_tip_lump(nselt, LIGHT) + raw = eigh( + 0.5 * (gk + gk.T), 0.5 * (gm + gm.T), subset_by_index=(0, 3), + )[0] + f_raw = float(eigvals_to_hz(raw, ROMG)[0]) + assert abs(f_raw - _analytic()) / _analytic() > 0.25 + + def test_backward_error_drops_after_the_retry(self): + _f, diag = _first_frequency(27, LIGHT) + assert diag.max_residual < 1.0e-6 + + +class TestHealthyProblemsAreUntouched: + """The guard must not fire on a mass distribution any real structure + would have — the standing rule for a numerical gate.""" + + @pytest.mark.parametrize("nselt", [13, 27, 53]) + def test_realistic_mass_ratio_stays_on_the_symmetric_path(self, nselt): + import warnings as _w + + gk, gm = _cantilever_with_tip_lump(nselt, REALISTIC) + with _w.catch_warnings(): + _w.simplefilter("error") + eigvals, _vecs, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.symmetric is True + f = float(eigvals_to_hz(eigvals, ROMG)[0]) + # The beam mass now matters, so this sits a little below the + # massless-beam closed form rather than on it. + assert 0.9 * _analytic() < f < _analytic() + + def test_healthy_residuals_keep_real_margin_under_the_threshold(self): + """A 400:1 lump-to-beam ratio is already a demanding model and is + the worst healthy backward error in the suite, at ~1.3e-3. The + threshold sits well above that and an order below the + degraded regime, so neither side is close. + + A false positive would cost time rather than accuracy in any + case: the retry keeps whichever of the two solves has the lower + backward error. + """ + from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt + + worst = 0.0 + for nselt in (13, 27, 53): + gk, gm = _cantilever_with_tip_lump(nselt, REALISTIC) + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + worst = max(worst, diag.max_residual) + assert worst < 0.05 * opt.residual_retry_threshold + + +class TestMarginalImprovementIsRefused: + """The decisive-win condition, pinned directly. + + A synthetic pair where the symmetric solve is above the threshold but + the general path cannot do materially better must keep the symmetric + result — the behaviour that protects the bundled land deck's + degenerate fore-aft / side-side pair. + """ + + def test_marginal_gain_keeps_the_symmetric_result(self, monkeypatch): + import dataclasses + + import pybmodes.fem.solver as solvermod + + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + + # Force the "general path is barely better" situation by making + # the improvement bar unreachable, leaving the threshold tripped. + monkeypatch.setattr( + solvermod, "_SOLVER_OPTIONS", + dataclasses.replace( + solvermod._SOLVER_OPTIONS, residual_retry_improvement=0.0, + ), + ) + import warnings as _w + + with _w.catch_warnings(): + _w.simplefilter("error") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" + + def test_degenerate_pair_survives_a_symmetric_solve(self): + """The property the improvement bar exists to protect: a tower + with EI_FA == EI_SS returns its bending modes as an exactly + degenerate pair, which the general path would split.""" + gk, gm = _cantilever_with_tip_lump(27, REALISTIC) + eigvals, _v = solve_modes(gk, gm, n_modes=4) + f = eigvals_to_hz(eigvals, ROMG) + assert f[0] == pytest.approx(f[1], rel=1.0e-9) + + +class TestRigidBodyModesAreNotMistakenForBreakdown: + """A free-free model's zero-frequency modes must not trip the guard. + + For a rigid-body mode ``K x ~ 0`` and ``lambda ~ 0``, so the relative + residual is a ratio of two near-zero quantities and evaluates to ~1 + however exact the eigenpair is. Judging a solve by the raw maximum + would condemn every floating result, and acting on that verdict is + worse than doing nothing: the general path filters out non-positive + eigenvalues, so it would delete the zero mode and shift every index + after it. + """ + + def _free_free_with_a_zero_mode(self): + """A symmetric platform with no yaw restoring: one exactly zero + eigenvalue among otherwise well-conditioned elastic modes.""" + n = 12 + rng = np.random.default_rng(7) + a = rng.normal(size=(n, n)) + m = a @ a.T + n * np.eye(n) # SPD, well conditioned + b = rng.normal(size=(n, n - 1)) + k = b @ b.T # rank n-1: one zero mode + return 0.5 * (k + k.T), 0.5 * (m + m.T) + + def test_zero_mode_present_and_no_retry(self): + import warnings as _w + + gk, gm = self._free_free_with_a_zero_mode() + with _w.catch_warnings(): + _w.simplefilter("error") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=6, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" + # The zero mode survives rather than being filtered away. + assert abs(eigvals[0]) < 1.0e-8 * abs(eigvals).max() + assert eigvals.size == 6 + + def test_raw_maximum_residual_would_have_condemned_it(self): + """The metric really is ~1 on the rigid-body mode, so excluding + it is what makes the guard usable rather than a nicety.""" + from pybmodes.fem.solver import _max_elastic_residual, _modal_residuals + + gk, gm = self._free_free_with_a_zero_mode() + eigvals, eigvecs = solve_modes(gk, gm, n_modes=6) + raw = _modal_residuals(gk, gm, eigvals, eigvecs) + assert raw.max() > 0.1 + assert _max_elastic_residual(gk, gm, eigvals, eigvecs) < 1.0e-8 + + def test_all_rigid_body_modes_report_zero(self): + """Nothing to judge means no verdict, not a bad one.""" + from pybmodes.fem.solver import _max_elastic_residual + + gk = np.zeros((4, 4)) + gm = np.eye(4) + vals = np.zeros(4) + vecs = np.eye(4) + assert _max_elastic_residual(gk, gm, vals, vecs) == 0.0 + + +class TestDiagnosticsContract: + def test_residual_fallback_defaults_to_false(self): + gk, gm = _cantilever_with_tip_lump(13, REALISTIC) + _v, _x, diag = solve_modes(gk, gm, n_modes=4, return_diagnostics=True) + assert diag.residual_fallback is False + + def test_plain_two_tuple_return_still_works(self): + """The retry must not change the historical return shape.""" + gk, gm = _cantilever_with_tip_lump(13, LIGHT) + with pytest.warns(RuntimeWarning): + out = solve_modes(gk, gm, n_modes=4) + assert isinstance(out, tuple) + assert len(out) == 2 From 7873c56bff6b9e87e7543a8118a9c7c0ff9fa694 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 17:22:48 +0900 Subject: [PATCH 02/28] fix: let the decisive-win rule handle rigid-body modes, not a classifier Codex P1 on #140. My rigid-body exclusion keyed its scale off the returned eigenvalues, so a subset containing only rigid-body modes took its own numerical noise as the reference and re-admitted them as elastic, reading their ~1 residuals as a breakdown and retrying into the one path that discards them. Rewriting the classifier to key off strain instead was worse. A genuinely soft mode on a stiff structure carries little strain too: the 0.084 Hz lump mode of a 1e10 N.m^2 beam came in below the cutoff, so the guard stopped firing on the very case it was written for and the answer went back to being 538 percent wrong, silently. Nothing here reliably separates a rigid-body mode from a soft elastic one, so the classification is gone. Both candidate solves are measured the same way over every mode, and the retry still needs a decisive win, so a mode the metric cannot speak to says the same nothing twice and cannot tip the decision. What actually made a false positive dangerous was the fallback dropping zero eigenvalues, changing the mode set rather than just its precision. The retry now runs with keep_rigid_body=True, which clamps rounding-level negatives to zero and keeps them. The asymmetric production path is untouched and keeps the BModes-matching filter it is validated against. A rigid-only subset now returns the modes it was asked for at every size, and the point-mass case is exact again. --- CHANGELOG.md | 15 ++- VALIDATION.md | 2 +- src/pybmodes/fem/solver.py | 131 +++++++++++++++---------- tests/fem/test_ill_conditioned_mass.py | 123 ++++++++++++++++------- 4 files changed, 181 insertions(+), 90 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ee2f4f2..a86eb71 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,11 +32,16 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). without being broken, and on the bundled NREL 5MW land tower (whose adapter leaves the mass matrix at cond ~4e10) the general path is only 1.4× better while *splitting* a degenerate fore-aft / side-side pair - the symmetric solver resolves exactly. Rigid-body modes are excluded - from the check for the same reason — their relative residual is a ratio - of two near-zero quantities and is ~1 however exact the eigenpair is, - so judging a floating solve by the raw maximum would have condemned it - and then deleted its zero-frequency mode. + the symmetric solver resolves exactly. + + The same condition disposes of rigid-body modes, on which the residual + is a ratio of two near-zero quantities and reads ~1 however exact the + eigenpair is. They are not identified and excluded — nothing separates + them reliably from a genuinely soft mode — instead both candidates are + measured identically, so a mode the metric cannot speak to says the + same thing twice and cannot tip the decision. The retry also preserves + zero eigenvalues, so a spurious trigger on a free-free model costs a + little time rather than deleting a physical mode. - `SolverOptions` gains `residual_retry_threshold` and `residual_retry_improvement` for the two conditions above. diff --git a/VALIDATION.md b/VALIDATION.md index 1f18530..ff2a0c9 100644 --- a/VALIDATION.md +++ b/VALIDATION.md @@ -103,7 +103,7 @@ metrics: | Distributed Winkler soil bed vs the lumped mudline condensation (issue #118) | Psaroudakis et al. (2021) / Yu & Amdahl (2023) Eq. 25 is the exact static condensation of a constant-EI pile on a bed of rate $k = D_P E_{SO}$ | coupled 1st frequency, distributed bed vs lumped springs | < 1 % | 0.4 % (the embedded pile inertia the condensed form drops) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | | Distributed Winkler bed converges on the rigid mudline clamp (issue #118) | limit $E_{SO} \to \infty$; residual compliance scales as the elastic length $(4EI/k)^{1/4}$ | 1st frequency at $E_{SO} \times 10^8$ vs the clamped model | < 1 % | (within tol) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | | Near-singular mass matrix does not silently corrupt the low modes | cantilever carrying one lump: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1), at a 4e5:1 lump-to-beam mass ratio | 1st frequency, dense and sparse dispatch sizes | < 0.5 % | (within tol; the unguarded symmetric solve is 137 % out) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no | -| Rigid-body modes are not mistaken for a solver breakdown | construction (rank-deficient $K$, well-conditioned $M$) | zero mode retained; no retry triggered | exact | (within tol) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no | +| Rigid-body modes are not lost to the solver guard | construction (rank-deficient $K$, well-conditioned $M$; requested subsets of 1, 3, 6 and 10 modes) | zero modes retained, requested mode count returned | exact | (within tol) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no | | Deck-reader `n_nodes` refinement (issue #58) | Euler-Bernoulli closed form, modes 1-3, plus self-convergence $n{=}100$ vs $200$ | bending frequencies | < 1 %; self-convergence < 0.2 % | (within tol) | [`tests/test_refine_mesh.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_refine_mesh.py) | no | **Citations** (full author / year forms used in the table above). diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 00733fd..4b20610 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -41,24 +41,32 @@ singular, which a very light beam carrying a very heavy lump produces. The failure mode is silent: LAPACK returns confidently wrong low modes rather than raising. :func:`solve_modes` therefore checks the backward -error of the **elastic** modes of every symmetric solve (a rigid-body -mode has ``K x ≈ 0`` and ``λ ≈ 0``, so its relative residual is a ratio -of two near-zero quantities and carries no information — see -:func:`_max_elastic_residual`) and, when it exceeds +error of every symmetric solve and, when it exceeds :attr:`~pybmodes.options.SolverOptions.residual_retry_threshold`, tries the general path as well — taking its result only if it is better by :attr:`~pybmodes.options.SolverOptions.residual_retry_improvement`, and warning when it does. -That second condition is the load-bearing one. A real deck can sit above -the threshold without being broken (the bundled NREL 5MW land tower -reaches ~2e-2, its adapter leaving ``M`` at cond ~4e10), and there the -general path is only marginally better while *splitting* the degenerate -fore-aft / side-side pair the symmetric solver resolves exactly — which -the FA / SS classifier downstream depends on. A true breakdown is not -marginal: it improves by nine orders of magnitude. Demanding a decisive -win keeps every validated frequency untouched and still catches the -failure this guard exists for. +That second condition is the load-bearing one, and it is what lets the +check be simple. A real deck can sit above the threshold without being +broken (the bundled NREL 5MW land tower reaches ~2e-2, its adapter +leaving ``M`` at cond ~4e10), and there the general path is only +marginally better while *splitting* the degenerate fore-aft / side-side +pair the symmetric solver resolves exactly — which the FA / SS +classifier downstream depends on. A true breakdown is not marginal: it +improves by nine orders of magnitude. + +The same condition disposes of rigid-body modes, on which the residual +is a ratio of two near-zero quantities and reads ~1 however exact the +eigenpair is. Rather than trying to identify such modes — neither their +eigenvalue nor their strain distinguishes them reliably from a genuinely +soft mode — both candidates are measured identically, so a mode the +metric cannot speak to says the same thing twice and cannot tip the +decision. The retry additionally runs with ``keep_rigid_body=True`` so +the two candidates describe the same spectrum; otherwise a retry +triggered by a free-free model's zero modes would swap in a mode set +with those modes filtered out, which is worse than the imprecision it +was trying to fix. Note on the user-spec mode choice: ``eigsh(..., sigma=0, mode='buckling')`` reduces to ``OP = K^-1 K = I`` for ``sigma=0``, @@ -269,11 +277,17 @@ def solve_modes( # stays exact there. residual_fallback = False if sym: - worst = _max_elastic_residual(gk, gm, eigvals, eigvecs) + worst = _max_residual(gk, gm, eigvals, eigvecs) if worst > _SOLVER_OPTIONS.residual_retry_threshold: - alt_vals, alt_vecs = _solve_dense_general(gk, gm, n_modes) + # ``keep_rigid_body`` so the two candidates describe the same + # spectrum. Without it a free-free model's zero modes would be + # filtered out of the alternative only, and a retry triggered + # by those very modes would swap in a different mode set. + alt_vals, alt_vecs = _solve_dense_general( + gk, gm, n_modes, keep_rigid_body=True, + ) _normalize_columns_l2(alt_vecs) - alt_worst = _max_elastic_residual(gk, gm, alt_vals, alt_vecs) + alt_worst = _max_residual(gk, gm, alt_vals, alt_vecs) # Only take the general result on a decisive win. A marginal # one is not a breakdown, and switching for it would churn # validated frequencies and break the degenerate fore-aft / @@ -373,37 +387,30 @@ def _build_diagnostics( # one is a rigid-body mode: a free-free floating platform has up to six, # and an unrestrained DOF (a symmetric column's yaw) gives an exactly # zero one. -_RIGID_BODY_EIGVAL_RTOL = 1.0e-8 - - -def _max_elastic_residual( +def _max_residual( gk: np.ndarray, gm: np.ndarray, eigvals: np.ndarray, eigvecs: np.ndarray, ) -> float: - """Largest backward error over the **elastic** modes only. - - The relative residual ``||K x - λ M x|| / ||K x||`` is undefined for a - rigid-body mode: there ``K x ≈ 0`` and ``λ ≈ 0``, so it is a ratio of - two near-zero quantities and evaluates to ≈ 1 no matter how exact the - eigenpair is. A free-free floating model legitimately has up to six of - them, so judging a solve by the raw maximum would condemn every - floating result — and the general path drops zero eigenvalues - entirely, so acting on that verdict would delete a physically real - mode rather than improve anything. - - Rigid-body modes are therefore excluded before taking the maximum. - Returns ``0.0`` when every returned mode is rigid-body, i.e. when - there is nothing the metric can speak to. + """Largest per-mode backward error, or ``0.0`` for an empty solve. + + Deliberately taken over **every** returned mode, with no attempt to + classify them first. On a rigid-body mode the metric is a ratio of + two near-zero quantities and reads ~1 regardless of how exact the + eigenpair is, which invites excluding those modes — but nothing + reliably identifies them here. Their eigenvalue is only "near zero" + relative to a scale the returned subset may not contain; their strain + ``||K x||`` is small, but so is a genuinely soft mode's on a stiff + structure (the 0.08 Hz lump mode of a 1e10 N.m^2 beam carries less + strain than the rigid-body modes of a floating platform do). + Attempting either classification blinded this guard to a case it had + previously caught. + + The meaningless component cancels instead. Both candidate solves are + measured the same way, and the retry is accepted only on a decisive + improvement, so a mode on which the metric says nothing says the same + nothing twice and cannot tip the decision. """ - if eigvals.size == 0: - return 0.0 - scale = float(np.max(np.abs(eigvals))) - if scale <= 0.0: - return 0.0 - elastic = np.abs(eigvals) > _RIGID_BODY_EIGVAL_RTOL * scale - if not elastic.any(): - return 0.0 r = _modal_residuals(gk, gm, eigvals, eigvecs) - return float(r[elastic].max()) + return float(r.max()) if r.size else 0.0 def _modal_residuals( @@ -497,20 +504,44 @@ def _solve_dense_symmetric( return np.asarray(eigvals), np.asarray(eigvecs) +# A real eigenvalue this far below zero, relative to the spectrum's own +# magnitude, is a rigid-body mode sitting at zero plus rounding rather +# than a non-physical negative one. +_RIGID_BODY_NEGATIVE_RTOL = 1.0e-10 + + def _solve_dense_general( gk: np.ndarray, gm: np.ndarray, n_modes: int | None, + *, keep_rigid_body: bool = False, ) -> tuple[np.ndarray, np.ndarray]: """Dense LAPACK ``eig`` for genuinely asymmetric problems. Filters eigenvalues to the real, positive, finite subset (matches BModes - JJ's general-matrix path).""" + JJ's general-matrix path). + + ``keep_rigid_body`` additionally retains eigenvalues that sit at zero + to within rounding, clamping them to exactly zero. Off by default, so + the asymmetric production path keeps the BModes-matching filter it is + validated against. The retry path in :func:`solve_modes` turns it on: + a free-free model's zero-frequency modes are physical, and dropping + them there would replace a possibly-imprecise spectrum with a + structurally different one — which is a worse outcome than the + imprecision, and would make a false-positive retry actively harmful + rather than merely wasteful. + """ eigvals_all, eigvecs_all = eig(gk, gm) eigvals_real = np.real_if_close(eigvals_all, tol=1000) - valid = ( - np.isreal(eigvals_real) - & np.isfinite(eigvals_real.real) - & (eigvals_real.real > 0.0) - ) - eigvals = eigvals_real.real[valid] + real_finite = np.isreal(eigvals_real) & np.isfinite(eigvals_real.real) + vals = eigvals_real.real + + if keep_rigid_body and real_finite.any(): + scale = float(np.max(np.abs(vals[real_finite]))) + floor = -_RIGID_BODY_NEGATIVE_RTOL * scale if scale > 0.0 else 0.0 + valid = real_finite & (vals >= floor) + vals = np.where(vals < 0.0, 0.0, vals) + else: + valid = real_finite & (vals > 0.0) + + eigvals = vals[valid] eigvecs = np.real_if_close(eigvecs_all[:, valid], tol=1000).real order = np.argsort(eigvals) if n_modes is not None: diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index cef416e..b97cc23 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -200,15 +200,23 @@ def test_degenerate_pair_survives_a_symmetric_solve(self): class TestRigidBodyModesAreNotMistakenForBreakdown: - """A free-free model's zero-frequency modes must not trip the guard. + """A free-free model's zero-frequency modes must not corrupt the result. For a rigid-body mode ``K x ~ 0`` and ``lambda ~ 0``, so the relative - residual is a ratio of two near-zero quantities and evaluates to ~1 - however exact the eigenpair is. Judging a solve by the raw maximum - would condemn every floating result, and acting on that verdict is - worse than doing nothing: the general path filters out non-positive - eigenvalues, so it would delete the zero mode and shift every index - after it. + residual is a ratio of two near-zero quantities and reads ~1 however + exact the eigenpair is. Two attempts to *identify* such modes and + exclude them both failed — an eigenvalue-relative cutoff takes a + rigid-only subset's own noise as its scale, and a strain-relative one + cannot tell a rigid mode from a genuinely soft one (the 0.08 Hz lump + mode of a 1e10 N.m^2 beam carries less strain than a floating + platform's rigid modes do, and excluding it blinded the guard to a + case it had caught). + + So they are not identified at all. Both candidate solves are measured + the same way and the retry needs a decisive win, so a mode the metric + cannot speak to says the same nothing twice. The retry also preserves + zero eigenvalues, which is what makes a false positive merely wasteful + instead of destructive. """ def _free_free_with_a_zero_mode(self): @@ -222,41 +230,88 @@ def _free_free_with_a_zero_mode(self): k = b @ b.T # rank n-1: one zero mode return 0.5 * (k + k.T), 0.5 * (m + m.T) - def test_zero_mode_present_and_no_retry(self): - import warnings as _w - + def test_zero_mode_survives_whichever_path_runs(self): gk, gm = self._free_free_with_a_zero_mode() - with _w.catch_warnings(): - _w.simplefilter("error") - eigvals, _v, diag = solve_modes( - gk, gm, n_modes=6, return_diagnostics=True, - ) - assert diag.residual_fallback is False - assert diag.path == "dense_symmetric" - # The zero mode survives rather than being filtered away. + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=6, return_diagnostics=True, + ) + # The zero mode survives rather than being filtered away, and the + # spectrum keeps its full width either way. assert abs(eigvals[0]) < 1.0e-8 * abs(eigvals).max() assert eigvals.size == 6 + assert diag.n_returned == 6 - def test_raw_maximum_residual_would_have_condemned_it(self): - """The metric really is ~1 on the rigid-body mode, so excluding - it is what makes the guard usable rather than a nicety.""" - from pybmodes.fem.solver import _max_elastic_residual, _modal_residuals + def test_the_metric_really_does_read_about_one_there(self): + """Why the naive raw-maximum reading is untrustworthy, pinned so + the reasoning above stays anchored to a number.""" + from pybmodes.fem.solver import _modal_residuals gk, gm = self._free_free_with_a_zero_mode() eigvals, eigvecs = solve_modes(gk, gm, n_modes=6) raw = _modal_residuals(gk, gm, eigvals, eigvecs) - assert raw.max() > 0.1 - assert _max_elastic_residual(gk, gm, eigvals, eigvecs) < 1.0e-8 - - def test_all_rigid_body_modes_report_zero(self): - """Nothing to judge means no verdict, not a bad one.""" - from pybmodes.fem.solver import _max_elastic_residual - - gk = np.zeros((4, 4)) - gm = np.eye(4) - vals = np.zeros(4) - vecs = np.eye(4) - assert _max_elastic_residual(gk, gm, vals, vecs) == 0.0 + rigid = np.argmin(np.abs(eigvals)) + assert raw[rigid] > 0.1 + # Every elastic mode is exact, so the ~1 is the metric failing, + # not the solve. + elastic = np.ones(raw.size, dtype=bool) + elastic[rigid] = False + assert raw[elastic].max() < 1.0e-8 + + def test_the_retry_preserves_zero_eigenvalues(self): + """The property that makes a false positive harmless: the + alternative solve keeps the rigid modes rather than filtering + them, so both candidates describe the same spectrum.""" + from pybmodes.fem.solver import _solve_dense_general + + gk, gm = self._free_free_with_a_zero_mode() + dropped, _v = _solve_dense_general(gk, gm, 6) + kept, _w = _solve_dense_general(gk, gm, 6, keep_rigid_body=True) + assert np.min(np.abs(kept)) == 0.0 + assert np.min(np.abs(dropped)) > 0.0 + # The default filter loses the zero mode and shifts the rest up. + assert kept[1] == pytest.approx(dropped[0], rel=1.0e-9) + + def _six_rigid_dofs(self): + """A model with six genuinely free rigid-body DOFs above a set of + elastic ones — an unmoored floating platform in the limit.""" + n = 14 + rng = np.random.default_rng(11) + a = rng.normal(size=(n, n)) + m = a @ a.T + n * np.eye(n) + b = rng.normal(size=(n, n - 6)) + k = b @ b.T # rank n-6: six zero modes + return 0.5 * (k + k.T), 0.5 * (m + m.T) + + @pytest.mark.parametrize("n_modes", [1, 3, 6]) + def test_a_rigid_only_subset_keeps_its_modes(self, n_modes): + """The case that broke both classification attempts. + + Every requested mode is rigid-body, so the metric reads ~1 on all + of them and no reference scale drawn from the subset can say + otherwise. The retry may well run; what matters is that it cannot + take modes away, because it now preserves zero eigenvalues and + has to win decisively to be accepted at all. + """ + gk, gm = self._six_rigid_dofs() + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=n_modes, return_diagnostics=True, + ) + assert eigvals.size == n_modes + assert np.all(np.abs(eigvals) < 1.0e-8 * float(np.linalg.norm(gk))) + assert diag.n_returned == n_modes + + def test_a_mixed_subset_is_measured_on_its_elastic_modes(self): + """With elastic modes present the metric is meaningful again and + reports them as exact.""" + gk, gm = self._six_rigid_dofs() + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=10, return_diagnostics=True, + ) + assert eigvals.size == 10 + assert np.max(np.abs(eigvals)) > 1.0e-8 * float(np.linalg.norm(gk)) + # Six rigid modes at the bottom, four exact elastic ones above. + assert np.sum(np.abs(eigvals) < 1.0e-12 * np.max(eigvals)) == 6 + assert max(diag.residuals[6:]) < 1.0e-8 class TestDiagnosticsContract: From 8eda44ec4abda103c2942472ea3f15f3ef8c96e2 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 17:31:02 +0900 Subject: [PATCH 03/28] fix: compare the two solves per mode, not on their maxima Codex P1 on #140, and it invalidates the cancellation argument in my last commit. A rigid-body mode reads ~1 in both candidates, so on a maximum it floors the alternative: max(alt_r) stays near 1 and no improvement among the elastic modes can clear a tenth of max(sym_r) unless the symmetric solve is worse than ~10. A free-free model with an elastic mode corrupted to a backward error of ~0.8 sailed straight through, silently, which is the failure the guard exists to catch. The decision is now per mode. Rigid modes contribute ~1 against ~1 and register as no improvement; a corrupted elastic mode contributes ~0.8 against ~1e-9 and registers clearly. Both candidates are sorted over the same spectrum, since the retry already preserves rigid-body modes, so equal indices mean the same mode. An alternative that recovered fewer modes is refused outright: losing one is never an improvement. The regression test builds the case directly, bolting six free rigid DOFs onto the ill-conditioned cantilever and rotating the pair so the blocks are not separable, and pins the mechanism as well as the outcome: the alternative maximum really is floored above 0.1, a maxima comparison really does see no win, and a per-mode one really does. --- CHANGELOG.md | 16 ++-- src/pybmodes/fem/solver.py | 117 ++++++++++++++----------- tests/fem/test_ill_conditioned_mass.py | 72 +++++++++++++++ 3 files changed, 146 insertions(+), 59 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a86eb71..2e0ce1d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -34,14 +34,14 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). 1.4× better while *splitting* a degenerate fore-aft / side-side pair the symmetric solver resolves exactly. - The same condition disposes of rigid-body modes, on which the residual - is a ratio of two near-zero quantities and reads ~1 however exact the - eigenpair is. They are not identified and excluded — nothing separates - them reliably from a genuinely soft mode — instead both candidates are - measured identically, so a mode the metric cannot speak to says the - same thing twice and cannot tip the decision. The retry also preserves - zero eigenvalues, so a spurious trigger on a free-free model costs a - little time rather than deleting a physical mode. + The comparison is made **per mode** rather than on the two maxima, so + that rigid-body modes cannot distort it. Their backward error is a + ratio of two near-zero quantities and reads ~1 in both candidates + however exact each is; on a maximum that floors the alternative and + hides a genuinely corrupted elastic mode sitting alongside them, while + per mode it simply registers as no improvement. The retry also + preserves zero eigenvalues, so a spurious trigger on a free-free model + costs a little time rather than deleting a physical mode. - `SolverOptions` gains `residual_retry_threshold` and `residual_retry_improvement` for the two conditions above. diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 4b20610..2b0d6c9 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -56,14 +56,17 @@ classifier downstream depends on. A true breakdown is not marginal: it improves by nine orders of magnitude. -The same condition disposes of rigid-body modes, on which the residual -is a ratio of two near-zero quantities and reads ~1 however exact the -eigenpair is. Rather than trying to identify such modes — neither their -eigenvalue nor their strain distinguishes them reliably from a genuinely -soft mode — both candidates are measured identically, so a mode the -metric cannot speak to says the same thing twice and cannot tip the -decision. The retry additionally runs with ``keep_rigid_body=True`` so -the two candidates describe the same spectrum; otherwise a retry +The comparison is made **per mode** rather than on the two maxima, which +is what keeps rigid-body modes from distorting it. Their backward error +is a ratio of two near-zero quantities and reads ~1 in both candidates +however exact each is; on a maximum that puts a floor under the +alternative and hides a genuinely corrupted elastic mode alongside them, +while per mode they simply register as ~1 against ~1, i.e. no +improvement. Identifying such modes and excluding them was tried twice +and abandoned — neither their eigenvalue nor their strain separates them +reliably from a genuinely soft mode. The retry additionally runs with +``keep_rigid_body=True`` so the two candidates describe the same +spectrum and equal indices mean the same mode; otherwise a retry triggered by a free-free model's zero modes would swap in a mode set with those modes filtered out, which is worse than the imprecision it was trying to fix. @@ -277,8 +280,8 @@ def solve_modes( # stays exact there. residual_fallback = False if sym: - worst = _max_residual(gk, gm, eigvals, eigvecs) - if worst > _SOLVER_OPTIONS.residual_retry_threshold: + sym_r = _modal_residuals(gk, gm, eigvals, eigvecs) + if sym_r.size and float(sym_r.max()) > _SOLVER_OPTIONS.residual_retry_threshold: # ``keep_rigid_body`` so the two candidates describe the same # spectrum. Without it a free-free model's zero modes would be # filtered out of the alternative only, and a retry triggered @@ -287,25 +290,23 @@ def solve_modes( gk, gm, n_modes, keep_rigid_body=True, ) _normalize_columns_l2(alt_vecs) - alt_worst = _max_residual(gk, gm, alt_vals, alt_vecs) - # Only take the general result on a decisive win. A marginal - # one is not a breakdown, and switching for it would churn - # validated frequencies and break the degenerate fore-aft / - # side-side pairs the symmetric solver resolves exactly — the - # bundled NREL 5MW land deck does exactly that. A genuine - # breakdown improves by many orders, not by a factor. - if alt_worst < _SOLVER_OPTIONS.residual_retry_improvement * worst: + alt_r = _modal_residuals(gk, gm, alt_vals, alt_vecs) + improved = _decisively_improved_modes(sym_r, alt_r, alt_vals.size, + eigvals.size) + if improved.any(): + idx = int(np.argmax(np.where(improved, sym_r[:improved.size], 0.0))) warnings.warn( - f"the symmetric eigensolver returned modes with a " - f"backward error of {worst:.2e}, so the eigenpairs do " - f"not satisfy K x = lambda M x. Its Cholesky reduction " - f"of the mass matrix loses accuracy when that matrix is " - f"nearly singular, which a very light beam carrying a " - f"very heavy lump produces. Redone through the general " - f"dense path, which factorises neither matrix " - f"(backward error {alt_worst:.2e}); the returned modes " - f"come from that solve. Worth checking the mass " - f"distribution is the one you intended.", + f"the symmetric eigensolver returned " + f"{int(improved.sum())} mode(s) that do not satisfy " + f"K x = lambda M x — worst at index {idx}, backward " + f"error {sym_r[idx]:.2e} against {alt_r[idx]:.2e} from " + f"the general dense path. Its Cholesky reduction of the " + f"mass matrix loses accuracy when that matrix is nearly " + f"singular, which a very light beam carrying a very " + f"heavy lump produces. The returned modes come from the " + f"general solve, which factorises neither matrix. Worth " + f"checking the mass distribution is the one you " + f"intended.", RuntimeWarning, stacklevel=2, ) @@ -387,30 +388,44 @@ def _build_diagnostics( # one is a rigid-body mode: a free-free floating platform has up to six, # and an unrestrained DOF (a symmetric column's yaw) gives an exactly # zero one. -def _max_residual( - gk: np.ndarray, gm: np.ndarray, eigvals: np.ndarray, eigvecs: np.ndarray, -) -> float: - """Largest per-mode backward error, or ``0.0`` for an empty solve. - - Deliberately taken over **every** returned mode, with no attempt to - classify them first. On a rigid-body mode the metric is a ratio of - two near-zero quantities and reads ~1 regardless of how exact the - eigenpair is, which invites excluding those modes — but nothing - reliably identifies them here. Their eigenvalue is only "near zero" - relative to a scale the returned subset may not contain; their strain - ``||K x||`` is small, but so is a genuinely soft mode's on a stiff - structure (the 0.08 Hz lump mode of a 1e10 N.m^2 beam carries less - strain than the rigid-body modes of a floating platform do). - Attempting either classification blinded this guard to a case it had - previously caught. - - The meaningless component cancels instead. Both candidate solves are - measured the same way, and the retry is accepted only on a decisive - improvement, so a mode on which the metric says nothing says the same - nothing twice and cannot tip the decision. +def _decisively_improved_modes( + sym_r: np.ndarray, + alt_r: np.ndarray, + n_alt: int, + n_sym: int, +) -> np.ndarray: + """Which modes the general path solves decisively better, per mode. + + The comparison has to be **per mode**, not on the two maxima. A + rigid-body mode's backward error is a ratio of two near-zero + quantities and reads ~1 in *both* candidates however exact each is, + so it sets a floor under the alternative's maximum: with one present, + ``max(alt_r)`` stays near 1 and no amount of improvement elsewhere + can drive it below a tenth of ``max(sym_r)`` unless the symmetric + solve is worse than ~10. A free-free model with a genuinely corrupted + elastic mode at a backward error of ~0.8 would sail through, which is + exactly the breakdown this guard exists to catch. + + Comparing mode by mode removes the floor: the rigid modes contribute + ~1 against ~1 and register as no improvement, while a corrupted + elastic mode contributes ~0.8 against ~1e-9 and registers clearly. + Both candidates are sorted ascending over the same spectrum (the + retry preserves rigid-body modes for this reason), so equal indices + describe the same mode. + + Returns a boolean mask over the compared modes. Empty when the + alternative recovered fewer modes than the symmetric solve — losing a + mode is never an improvement, whatever the residuals say. """ - r = _modal_residuals(gk, gm, eigvals, eigvecs) - return float(r.max()) if r.size else 0.0 + if n_alt < n_sym: + return np.zeros(0, dtype=bool) + n = min(sym_r.size, alt_r.size) + if n == 0: + return np.zeros(0, dtype=bool) + return ( + (sym_r[:n] > _SOLVER_OPTIONS.residual_retry_threshold) + & (alt_r[:n] < _SOLVER_OPTIONS.residual_retry_improvement * sym_r[:n]) + ) def _modal_residuals( diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index b97cc23..2ea3dad 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -314,6 +314,78 @@ def test_a_mixed_subset_is_measured_on_its_elastic_modes(self): assert max(diag.residuals[6:]) < 1.0e-8 +class TestRigidModesCannotMaskAnElasticBreakdown: + """A corrupted elastic mode must be caught even when rigid-body modes + sit alongside it. + + Comparing the two candidates on their *maxima* fails here: a + rigid-body mode reads ~1 in both, so it floors the alternative's + maximum and no improvement among the elastic modes can clear a + decisive-win bar unless the symmetric solve is worse than ~10. A + free-free model with an elastic mode corrupted to a backward error of + ~0.8 would then pass silently — the exact failure this guard exists + to catch, reintroduced by the rigid modes' presence. + + Comparing per mode removes the floor. + """ + + def _rigid_plus_ill_conditioned(self, nselt: int = 27): + """Six free rigid DOFs bolted onto the ill-conditioned cantilever, + then rotated so the two blocks are not separable by inspection.""" + gk, gm = _cantilever_with_tip_lump(nselt, LIGHT) + n = gk.shape[0] + big_k = np.zeros((n + 6, n + 6)) + big_m = np.zeros((n + 6, n + 6)) + big_k[:n, :n] = gk + big_m[:n, :n] = gm + # Zero stiffness, unit mass on the six extra DOFs: rigid-body. + big_m[n:, n:] = np.eye(6) * float(np.trace(gm)) / n + q, _ = np.linalg.qr(np.random.default_rng(3).normal(size=(n + 6, n + 6))) + k_rot = q.T @ big_k @ q + m_rot = q.T @ big_m @ q + return 0.5 * (k_rot + k_rot.T), 0.5 * (m_rot + m_rot.T) + + def test_the_rigid_modes_really_do_floor_the_maximum(self): + """The mechanism, pinned: on the maxima the alternative cannot + look decisively better even though it is exact where it counts.""" + from pybmodes.fem.solver import _modal_residuals, _solve_dense_general + + gk, gm = self._rigid_plus_ill_conditioned() + from scipy.linalg import eigh + + w, v = eigh(gk, gm, subset_by_index=(0, 9)) + v = v / np.linalg.norm(v, axis=0) + sym_r = _modal_residuals(gk, gm, w, v) + aw, av = _solve_dense_general(gk, gm, 10, keep_rigid_body=True) + av = av / np.linalg.norm(av, axis=0) + alt_r = _modal_residuals(gk, gm, aw, av) + # The alternative's maximum is pinned near 1 by the rigid modes... + assert alt_r.max() > 0.1 + # ...so a maxima comparison sees no decisive win. + assert not alt_r.max() < 0.1 * sym_r.max() + # ...yet some mode really is corrupted and really is fixed. + assert ((sym_r > 0.1) & (alt_r < 0.1 * sym_r)).any() + + def test_the_breakdown_is_caught_and_corrected(self): + gk, gm = self._rigid_plus_ill_conditioned() + with pytest.warns(RuntimeWarning, match="do not satisfy"): + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=10, return_diagnostics=True, + ) + assert diag.residual_fallback is True + assert eigvals.size == 10 + + def test_the_corrected_spectrum_matches_the_underlying_one(self): + """The rotation and the extra rigid DOFs do not change the + cantilever's own eigenvalues, so the corrected solve must still + contain the analytic lump frequency.""" + gk, gm = self._rigid_plus_ill_conditioned() + with pytest.warns(RuntimeWarning): + eigvals, _v = solve_modes(gk, gm, n_modes=10) + f = eigvals_to_hz(eigvals, ROMG) + assert np.min(np.abs(f - _analytic())) < 5.0e-3 * _analytic() + + class TestDiagnosticsContract: def test_residual_fallback_defaults_to_false(self): gk, gm = _cantilever_with_tip_lump(13, REALISTIC) From fbf9e06a8453b881d4d2f7b6323383173c16648f Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 17:40:25 +0900 Subject: [PATCH 04/28] fix: keep the whole spectrum on the retry, and depend less on the BLAS Two things. Codex P2 on #140: preserving only rounding-level zeros still filtered genuinely negative eigenvalues, which an indefinite K produces once run(gravity=...) loads a column past its buckling weight. The alternative then came back the same length, backfilled from higher up, while describing a shifted spectrum, and the per-index comparison was reading two different spectra against each other. eigh filters nothing, so the retry now filters nothing either: keep every real finite eigenvalue, sign included. Both paths return the n_modes smallest real eigenvalues, so equal indices describe the same mode by construction rather than by assumption. The asymmetric production path keeps its positive-only filter. And two of my own tests failed on the Linux BLAS while passing here, because they asserted LAPACK behaviour rather than pyBmodes logic. The degenerate-pair test wanted a split below 1e-9 where one build gives 1e-16 and another 3e-4; it now asks only that the pair is still one mode to engineering precision, which is the property that matters. The masking mechanism is now stated as arithmetic on residual vectors instead of being coaxed out of a matrix pair, so it pins the two decision rules exactly and cannot drift with the vendor library. --- src/pybmodes/fem/solver.py | 79 +++++++-------- tests/fem/test_ill_conditioned_mass.py | 130 ++++++++++++++++++++----- 2 files changed, 145 insertions(+), 64 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 2b0d6c9..4871392 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -64,12 +64,18 @@ while per mode they simply register as ~1 against ~1, i.e. no improvement. Identifying such modes and excluding them was tried twice and abandoned — neither their eigenvalue nor their strain separates them -reliably from a genuinely soft mode. The retry additionally runs with -``keep_rigid_body=True`` so the two candidates describe the same -spectrum and equal indices mean the same mode; otherwise a retry -triggered by a free-free model's zero modes would swap in a mode set -with those modes filtered out, which is worse than the imprecision it -was trying to fix. +reliably from a genuinely soft mode. + +The retry additionally runs with ``preserve_full_spectrum=True``, which +drops the sign filter the general path normally applies. ``eigh`` +filters nothing, so keeping it would return a *different set* of modes — +the same length, since the gap is backfilled from higher up — and the +per-index comparison would be reading two different spectra against each +other, able to accept a result that had quietly dropped a mode and +shifted every one above it. Both omissions are reachable: a free-free +model's zero-frequency modes, and the negative eigenvalues an indefinite +``K`` produces once ``run(gravity=...)`` loads a column past its +buckling weight. Note on the user-spec mode choice: ``eigsh(..., sigma=0, mode='buckling')`` reduces to ``OP = K^-1 K = I`` for ``sigma=0``, @@ -282,12 +288,14 @@ def solve_modes( if sym: sym_r = _modal_residuals(gk, gm, eigvals, eigvecs) if sym_r.size and float(sym_r.max()) > _SOLVER_OPTIONS.residual_retry_threshold: - # ``keep_rigid_body`` so the two candidates describe the same - # spectrum. Without it a free-free model's zero modes would be - # filtered out of the alternative only, and a retry triggered - # by those very modes would swap in a different mode set. + # ``preserve_full_spectrum`` so the two candidates describe the + # same spectrum and equal indices mean the same mode. ``eigh`` + # filters nothing, so any sign filter here would return a + # different set — same length, backfilled from higher up — and + # the per-index comparison would then be reading two different + # spectra against each other. alt_vals, alt_vecs = _solve_dense_general( - gk, gm, n_modes, keep_rigid_body=True, + gk, gm, n_modes, preserve_full_spectrum=True, ) _normalize_columns_l2(alt_vecs) alt_r = _modal_residuals(gk, gm, alt_vals, alt_vecs) @@ -519,44 +527,39 @@ def _solve_dense_symmetric( return np.asarray(eigvals), np.asarray(eigvecs) -# A real eigenvalue this far below zero, relative to the spectrum's own -# magnitude, is a rigid-body mode sitting at zero plus rounding rather -# than a non-physical negative one. -_RIGID_BODY_NEGATIVE_RTOL = 1.0e-10 - - def _solve_dense_general( gk: np.ndarray, gm: np.ndarray, n_modes: int | None, - *, keep_rigid_body: bool = False, + *, preserve_full_spectrum: bool = False, ) -> tuple[np.ndarray, np.ndarray]: """Dense LAPACK ``eig`` for genuinely asymmetric problems. Filters eigenvalues to the real, positive, finite subset (matches BModes JJ's general-matrix path). - ``keep_rigid_body`` additionally retains eigenvalues that sit at zero - to within rounding, clamping them to exactly zero. Off by default, so - the asymmetric production path keeps the BModes-matching filter it is - validated against. The retry path in :func:`solve_modes` turns it on: - a free-free model's zero-frequency modes are physical, and dropping - them there would replace a possibly-imprecise spectrum with a - structurally different one — which is a worse outcome than the - imprecision, and would make a false-positive retry actively harmful - rather than merely wasteful. + ``preserve_full_spectrum`` drops the **sign** filter, keeping every + real finite eigenvalue including zeros and negatives. Off by default, + so the asymmetric production path keeps the BModes-matching filter it + is validated against. + + The retry path in :func:`solve_modes` turns it on, and needs the + whole spectrum rather than just its zeros. ``eigh`` filters nothing, + so any sign filter here would return a *different set* of modes — + same length, since the gap is backfilled from higher up — and the + per-index comparison would then be reading two different spectra + against each other, accepting a result that had quietly dropped a + mode and shifted every one above it. Both cases are real: a + free-free model's zero-frequency modes, and the genuinely negative + eigenvalues an indefinite ``K`` produces once ``run(gravity=...)`` + loads a column past its buckling weight. With no filter at all both + paths return the ``n_modes`` smallest real eigenvalues, so equal + indices describe the same mode by construction. """ eigvals_all, eigvecs_all = eig(gk, gm) eigvals_real = np.real_if_close(eigvals_all, tol=1000) - real_finite = np.isreal(eigvals_real) & np.isfinite(eigvals_real.real) - vals = eigvals_real.real - - if keep_rigid_body and real_finite.any(): - scale = float(np.max(np.abs(vals[real_finite]))) - floor = -_RIGID_BODY_NEGATIVE_RTOL * scale if scale > 0.0 else 0.0 - valid = real_finite & (vals >= floor) - vals = np.where(vals < 0.0, 0.0, vals) - else: - valid = real_finite & (vals > 0.0) + valid = np.isreal(eigvals_real) & np.isfinite(eigvals_real.real) + if not preserve_full_spectrum: + valid = valid & (eigvals_real.real > 0.0) - eigvals = vals[valid] + eigvals = eigvals_real.real[valid] eigvecs = np.real_if_close(eigvecs_all[:, valid], tol=1000).real order = np.argsort(eigvals) if n_modes is not None: diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 2ea3dad..576f331 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -28,6 +28,8 @@ from __future__ import annotations +import warnings + import numpy as np import pytest @@ -189,14 +191,21 @@ def test_marginal_gain_keeps_the_symmetric_result(self, monkeypatch): assert diag.residual_fallback is False assert diag.path == "dense_symmetric" - def test_degenerate_pair_survives_a_symmetric_solve(self): + def test_degenerate_pair_comes_back_as_a_pair(self): """The property the improvement bar exists to protect: a tower - with EI_FA == EI_SS returns its bending modes as an exactly - degenerate pair, which the general path would split.""" + with EI_FA == EI_SS returns its bending modes as one degenerate + pair rather than two separated modes. + + The tolerance is loose on purpose. How exactly the pair resolves + depends on the LAPACK build — this same model splits it at ~1e-16 + on one and ~3e-4 on another — so pinning it tightly tests the + vendor's BLAS rather than pyBmodes. What matters here is that the + two remain the same mode to engineering precision. + """ gk, gm = _cantilever_with_tip_lump(27, REALISTIC) eigvals, _v = solve_modes(gk, gm, n_modes=4) f = eigvals_to_hz(eigvals, ROMG) - assert f[0] == pytest.approx(f[1], rel=1.0e-9) + assert f[0] == pytest.approx(f[1], rel=1.0e-2) class TestRigidBodyModesAreNotMistakenForBreakdown: @@ -265,11 +274,11 @@ def test_the_retry_preserves_zero_eigenvalues(self): gk, gm = self._free_free_with_a_zero_mode() dropped, _v = _solve_dense_general(gk, gm, 6) - kept, _w = _solve_dense_general(gk, gm, 6, keep_rigid_body=True) - assert np.min(np.abs(kept)) == 0.0 - assert np.min(np.abs(dropped)) > 0.0 + kept, _w = _solve_dense_general(gk, gm, 6, preserve_full_spectrum=True) + assert abs(kept[0]) < 1.0e-8 * float(np.max(np.abs(kept))) + assert np.min(dropped) > 0.0 # The default filter loses the zero mode and shifts the rest up. - assert kept[1] == pytest.approx(dropped[0], rel=1.0e-9) + assert kept[1] == pytest.approx(dropped[0], rel=1.0e-6) def _six_rigid_dofs(self): """A model with six genuinely free rigid-body DOFs above a set of @@ -345,26 +354,36 @@ def _rigid_plus_ill_conditioned(self, nselt: int = 27): m_rot = q.T @ big_m @ q return 0.5 * (k_rot + k_rot.T), 0.5 * (m_rot + m_rot.T) - def test_the_rigid_modes_really_do_floor_the_maximum(self): - """The mechanism, pinned: on the maxima the alternative cannot - look decisively better even though it is exact where it counts.""" - from pybmodes.fem.solver import _modal_residuals, _solve_dense_general + def test_the_maxima_rule_misses_what_the_per_mode_rule_catches(self): + """The mechanism, pinned on residual vectors directly. - gk, gm = self._rigid_plus_ill_conditioned() - from scipy.linalg import eigh + Stated as arithmetic rather than run through LAPACK on purpose: + whether a given matrix pair happens to exhibit the masking + depends on how badly that build's BLAS degrades, which is not + what this is about. These are the numbers the two rules see — + two rigid-body modes reading ~1 in both candidates, one elastic + mode corrupted to 0.8 and fixed to 1e-9, one mode already exact. + """ + from pybmodes.fem.solver import _decisively_improved_modes + + sym_r = np.array([1.0, 1.0, 0.8, 1.0e-12]) + alt_r = np.array([1.0, 1.0, 1.0e-9, 1.0e-12]) + + # The rigid modes floor the alternative's maximum at ~1, so a + # maxima rule sees no decisive win and leaves the corruption in. + assert not (alt_r.max() < 0.1 * sym_r.max()) - w, v = eigh(gk, gm, subset_by_index=(0, 9)) - v = v / np.linalg.norm(v, axis=0) - sym_r = _modal_residuals(gk, gm, w, v) - aw, av = _solve_dense_general(gk, gm, 10, keep_rigid_body=True) - av = av / np.linalg.norm(av, axis=0) - alt_r = _modal_residuals(gk, gm, aw, av) - # The alternative's maximum is pinned near 1 by the rigid modes... - assert alt_r.max() > 0.1 - # ...so a maxima comparison sees no decisive win. - assert not alt_r.max() < 0.1 * sym_r.max() - # ...yet some mode really is corrupted and really is fixed. - assert ((sym_r > 0.1) & (alt_r < 0.1 * sym_r)).any() + # Per mode, the corrupted one is unmissable and the rigid ones + # register as exactly what they are: no improvement either way. + improved = _decisively_improved_modes(sym_r, alt_r, 4, 4) + assert improved.tolist() == [False, False, True, False] + + def test_a_shorter_alternative_is_never_an_improvement(self): + from pybmodes.fem.solver import _decisively_improved_modes + + sym_r = np.array([1.0, 0.8, 1.0e-12]) + alt_r = np.array([1.0e-9, 1.0e-9]) + assert not _decisively_improved_modes(sym_r, alt_r, 2, 3).any() def test_the_breakdown_is_caught_and_corrected(self): gk, gm = self._rigid_plus_ill_conditioned() @@ -386,6 +405,65 @@ def test_the_corrected_spectrum_matches_the_underlying_one(self): assert np.min(np.abs(f - _analytic())) < 5.0e-3 * _analytic() +class TestNegativeEigenvaluesSurviveTheRetry: + """An indefinite ``K`` must not have its unstable modes filtered away. + + ``run(gravity=...)`` past a column's buckling weight drives + eigenvalues genuinely negative, and ``eigh`` returns them as-is. If + the alternative solve filtered them, it would come back the same + *length* — the gap backfilled from higher up — while describing a + shifted spectrum, and the per-index comparison would be reading two + different spectra against each other. + """ + + def _indefinite(self): + """The ill-conditioned cantilever with a genuinely negative mode + bolted on and the pair rotated together.""" + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + n = gk.shape[0] + big_k = np.zeros((n + 2, n + 2)) + big_m = np.zeros((n + 2, n + 2)) + big_k[:n, :n] = gk + big_m[:n, :n] = gm + scale = float(np.trace(gk)) / n + big_k[n, n] = -scale # unstable + big_k[n + 1, n + 1] = scale + big_m[n:, n:] = np.eye(2) * float(np.trace(gm)) / n + q, _ = np.linalg.qr(np.random.default_rng(5).normal(size=(n + 2, n + 2))) + k_rot, m_rot = q.T @ big_k @ q, q.T @ big_m @ q + return 0.5 * (k_rot + k_rot.T), 0.5 * (m_rot + m_rot.T) + + def test_the_default_filter_would_drop_the_unstable_mode(self): + from pybmodes.fem.solver import _solve_dense_general + + gk, gm = self._indefinite() + dropped, _v = _solve_dense_general(gk, gm, 6) + kept, _w = _solve_dense_general(gk, gm, 6, preserve_full_spectrum=True) + assert kept.min() < 0.0 # the unstable mode is present + assert dropped.min() > 0.0 # and absent from the default + # Same length, shifted spectrum — the trap the index comparison + # would otherwise walk into. + assert kept.size == dropped.size + assert kept[1] == pytest.approx(dropped[0], rel=1.0e-6) + + def test_the_symmetric_solve_reports_it_too(self): + """Both paths must agree that the mode exists, or index matching + is meaningless.""" + from scipy.linalg import eigh + + gk, gm = self._indefinite() + w = eigh(gk, gm, subset_by_index=(0, 5))[0] + assert w.min() < 0.0 + + def test_solve_modes_keeps_the_unstable_mode(self): + gk, gm = self._indefinite() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + eigvals, _v = solve_modes(gk, gm, n_modes=6) + assert eigvals.size == 6 + assert eigvals.min() < 0.0 + + class TestDiagnosticsContract: def test_residual_fallback_defaults_to_false(self): gk, gm = _cantilever_with_tip_lump(13, REALISTIC) From b724726aeb4943e125b8b311d1a3fa103da706ea Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 17:48:01 +0900 Subject: [PATCH 05/28] fix: measure and retry against the matrices the symmetric path solved Codex P1 on #140. The symmetric paths symmetrise internally, but the guard judged their modes against the raw matrices. The skew that _is_effectively_symmetric tolerates is only small relative to max|K|, so in a model with a wide dynamic range it can be the same size as a soft mode own eigenvalue. An exact symmetric solve then reads as broken, and eig on those same raw matrices wins decisively purely by answering a different question, replacing the symmetric spectrum the caller was promised with the skewed one. Both the measurement and the retry now use the symmetrised pair. The regression test builds a 3x3 with eigenvalues from 1 down to 1e-12 and skew just inside the tolerance, coupling the two softest modes so the perturbation is comparable to the eigenvalue rather than lost against it. It asserts the pair is accepted as symmetric, that no retry fires, that the returned spectrum is the symmetric one, and that measuring the very same modes against the raw matrices would have read above the threshold. --- src/pybmodes/fem/solver.py | 26 ++++++++-- tests/fem/test_ill_conditioned_mass.py | 66 ++++++++++++++++++++++++++ 2 files changed, 89 insertions(+), 3 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 4871392..0bc974c 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -77,6 +77,15 @@ ``K`` produces once ``run(gravity=...)`` loads a column past its buckling weight. +Both the measurement and the retry use the **symmetrised** matrices, the +ones the symmetric paths actually solve. The skew they discard is only +guaranteed small relative to ``max|K|``, which in a model with a wide +dynamic range can still be large relative to a soft mode's own +eigenvalue; judging an exact symmetric solve against the unsymmetrised +matrices would then read as a failure, and ``eig`` on those same +matrices would "win decisively" purely by answering a different +question. + Note on the user-spec mode choice: ``eigsh(..., sigma=0, mode='buckling')`` reduces to ``OP = K^-1 K = I`` for ``sigma=0``, which is degenerate. The standard scipy idiom for "smallest @@ -286,7 +295,18 @@ def solve_modes( # stays exact there. residual_fallback = False if sym: - sym_r = _modal_residuals(gk, gm, eigvals, eigvecs) + # Measure — and retry — against the matrices the symmetric paths + # actually solved. Both symmetrise internally, and the accepted + # skew is only guaranteed small relative to ``max|K|``: in a model + # with a wide dynamic range it can still be large relative to a + # soft mode's own eigenvalue. Judging an exact symmetric solve + # against the unsymmetrised matrices would then show a residual + # above the threshold, and ``eig`` on those same unsymmetrised + # matrices would "win decisively" purely by answering a different + # question — replacing a correct spectrum with the skew's. + gk_s = 0.5 * (gk + gk.T) + gm_s = 0.5 * (gm + gm.T) + sym_r = _modal_residuals(gk_s, gm_s, eigvals, eigvecs) if sym_r.size and float(sym_r.max()) > _SOLVER_OPTIONS.residual_retry_threshold: # ``preserve_full_spectrum`` so the two candidates describe the # same spectrum and equal indices mean the same mode. ``eigh`` @@ -295,10 +315,10 @@ def solve_modes( # the per-index comparison would then be reading two different # spectra against each other. alt_vals, alt_vecs = _solve_dense_general( - gk, gm, n_modes, preserve_full_spectrum=True, + gk_s, gm_s, n_modes, preserve_full_spectrum=True, ) _normalize_columns_l2(alt_vecs) - alt_r = _modal_residuals(gk, gm, alt_vals, alt_vecs) + alt_r = _modal_residuals(gk_s, gm_s, alt_vals, alt_vecs) improved = _decisively_improved_modes(sym_r, alt_r, alt_vals.size, eigvals.size) if improved.any(): diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 576f331..2872997 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -464,6 +464,72 @@ def test_solve_modes_keeps_the_unstable_mode(self): assert eigvals.min() < 0.0 +class TestAcceptedSkewDoesNotTriggerTheRetry: + """A tolerated asymmetry must not be read as a solver failure. + + ``_is_effectively_symmetric`` accepts skew up to ``symmetry_rtol`` + times ``max|K|``, and the symmetric paths then solve the symmetrised + matrices. That skew is only small *relative to the largest* entry: in + a model with a wide dynamic range it can be comparable to a soft + mode's own eigenvalue. Measuring the resulting modes against the + unsymmetrised matrices makes an exact solve look broken, and the + general path then "wins decisively" only because it is answering a + different question — returning the skewed spectrum in place of the + symmetric one the caller was promised. + """ + + def _wide_range_with_accepted_skew(self): + from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt + + # Eigenvalues spanning 1 down to 1e-12, so max|K| is 1 and the + # tolerated skew is ~1e-12 — the same size as the softest mode. + # The symmetry test compares max|A - A.T|, which is twice the + # off-diagonal skew, so stay under half the tolerance. + d = np.array([1.0, 1.0e-6, 1.0e-12]) + gm = np.eye(3) + s = 0.4 * opt.symmetry_rtol * max(1.0, float(np.max(np.abs(d)))) + # Couple the two softest modes: skew between the stiff ones would + # be negligible against their own scale and prove nothing. + gk = np.diag(d) + np.array([[0.0, 0.0, 0.0], + [0.0, 0.0, s], + [0.0, -s, 0.0]]) + return gk, gm, d + + def test_the_pair_is_accepted_as_symmetric(self): + from pybmodes.fem.solver import _is_effectively_symmetric + + gk, gm, _d = self._wide_range_with_accepted_skew() + assert _is_effectively_symmetric(gk) + assert _is_effectively_symmetric(gm) + + def test_no_retry_and_the_symmetric_spectrum_is_returned(self): + gk, gm, d = self._wide_range_with_accepted_skew() + with warnings.catch_warnings(): + warnings.simplefilter("error") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=3, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" + # The symmetrised problem's spectrum, not the skewed one. + assert np.allclose(np.sort(eigvals), np.sort(d), rtol=1.0e-6) + + def test_measuring_against_the_unsymmetrised_pair_would_have_tripped(self): + """The mechanism: the same exact modes look broken when judged + against matrices they were never solved on.""" + from pybmodes.fem.solver import _modal_residuals + + gk, gm, _d = self._wide_range_with_accepted_skew() + gk_s = 0.5 * (gk + gk.T) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + eigvals, eigvecs = solve_modes(gk, gm, n_modes=3) + against_solved = _modal_residuals(gk_s, gm, eigvals, eigvecs) + against_raw = _modal_residuals(gk, gm, eigvals, eigvecs) + assert against_solved.max() < 1.0e-8 + assert against_raw.max() > 0.1 + + class TestDiagnosticsContract: def test_residual_fallback_defaults_to_false(self): gk, gm = _cantilever_with_tip_lump(13, REALISTIC) From 55e30ef985d969ff2038634e8822374f78f35c52 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 17:54:31 +0900 Subject: [PATCH 06/28] fix: report diagnostics against the problem the modes actually solved Codex P2 on #140, and a direct consequence of a call I made in the last commit. I moved the retry decision onto the symmetrised matrices but deliberately left _build_diagnostics on the raw ones, to avoid perturbing a documented field. That was the wrong trade: with accepted skew present, a correct symmetric solve reported max_residual around 0.37, so the field advertised as certification telemetry was flagging a healthy result as defective. solve_modes now names the pair the returned modes actually solve once, symmetrised on the symmetric paths and raw on the general one, and uses it for the retry decision and the diagnostics alike. The field docstring says which basis it is on, since a backward error is meaningless without naming the problem it is against. Genuinely symmetric decks are unaffected: with no skew the symmetrised pair is the raw pair. --- src/pybmodes/fem/solver.py | 23 +++++++++++++++++++---- tests/fem/test_ill_conditioned_mass.py | 12 ++++++++++++ 2 files changed, 31 insertions(+), 4 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 0bc974c..979c9bb 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -159,7 +159,16 @@ class SolverDiagnostics: when no modes were returned). A healthy modal solve sits near machine precision; a large value flags an ill-conditioned or defective eigenproblem. - residuals : the per-mode relative residuals, one per returned mode. + + Measured against the matrices the returned modes **actually + solve**: the symmetrised pair on the symmetric paths, which + symmetrise internally, and the raw pair on the general one. + Measuring a symmetric solve against the raw matrices would charge + it for the skew it was told to discard — on a model with a wide + dynamic range that reads as a large backward error for a solve + that is exact, which is the opposite of what this field is for. + residuals : the per-mode relative residuals, one per returned mode, + on the same basis as ``max_residual``. matrix_cond : 2-norm condition number of the (symmetrised) mass matrix, or ``None`` when not computed (sparse path, or a system larger than the dense-conditioning size limit). @@ -293,6 +302,13 @@ def solve_modes( # catches it (healthy solves sit at ~1e-4 or below, degraded ones # above 1), and the general path, which factorises neither matrix, # stays exact there. + # The matrices the returned modes actually solve. Both symmetric + # paths symmetrise internally, so for them the diagnostics — and the + # retry decision below — have to be measured against that pair, not + # against the raw one. Reporting the raw backward error would flag a + # correct solve as defective in telemetry meant to be auditable. + res_k, res_m = (0.5 * (gk + gk.T), 0.5 * (gm + gm.T)) if sym else (gk, gm) + residual_fallback = False if sym: # Measure — and retry — against the matrices the symmetric paths @@ -304,8 +320,7 @@ def solve_modes( # above the threshold, and ``eig`` on those same unsymmetrised # matrices would "win decisively" purely by answering a different # question — replacing a correct spectrum with the skew's. - gk_s = 0.5 * (gk + gk.T) - gm_s = 0.5 * (gm + gm.T) + gk_s, gm_s = res_k, res_m sym_r = _modal_residuals(gk_s, gm_s, eigvals, eigvecs) if sym_r.size and float(sym_r.max()) > _SOLVER_OPTIONS.residual_retry_threshold: # ``preserve_full_spectrum`` so the two candidates describe the @@ -374,7 +389,7 @@ def solve_modes( return eigvals, eigvecs diagnostics = _build_diagnostics( - gk, gm, eigvals, eigvecs, path=path, symmetric=sym, + res_k, res_m, eigvals, eigvecs, path=path, symmetric=sym, n_requested=n_modes, sparse_fallback=sparse_fallback, fallback_reason=fallback_reason, residual_fallback=residual_fallback, diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 2872997..5663858 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -514,6 +514,18 @@ def test_no_retry_and_the_symmetric_spectrum_is_returned(self): # The symmetrised problem's spectrum, not the skewed one. assert np.allclose(np.sort(eigvals), np.sort(d), rtol=1.0e-6) + def test_the_diagnostics_report_the_solved_problem_too(self): + """Telemetry meant to be auditable must not charge a correct + solve for the skew it was told to discard.""" + gk, gm, _d = self._wide_range_with_accepted_skew() + with warnings.catch_warnings(): + warnings.simplefilter("error") + _v, _x, diag = solve_modes(gk, gm, n_modes=3, + return_diagnostics=True) + assert diag.residual_fallback is False + assert diag.max_residual < 1.0e-8 + assert max(diag.residuals) < 1.0e-8 + def test_measuring_against_the_unsymmetrised_pair_would_have_tripped(self): """The mechanism: the same exact modes look broken when judged against matrices they were never solved on.""" From 0fc9e0cb94606d1d7a5ca4337ee8f7f22d8fdaa3 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 18:03:51 +0900 Subject: [PATCH 07/28] fix: verify the retry ordering, and state where the guard cannot help Codex P1 on #140, which I had independently reached from the other direction while reviewing my own work. QZ can represent a symmetric problem theoretically real eigenvalues as small complex-conjugate pairs that real_if_close will not coerce, so the real filter deletes them even with no sign filter, and a truncated request backfills the gap from higher up. Equal counts then prove nothing about equal indices, which was the assumption the whole per-index comparison rested on. My first attempt demanded that every mode survive. That is too strict to be useful: on this machine eig drops six of 249 on a healthy fixture, all at the stiff end, far above anything a caller asks for. The invariant that actually matters is narrower, so the retry now verifies it directly. Nothing discarded may fall inside the returned window, and when something does the retry declines rather than swapping. Codex evidence shows why declining is the right answer rather than a weaker one: on their build the dropped pairs are the zero modes themselves, and keeping the alternative would backfill them with elastic modes and silently shift the spectrum. A guard added to stop a silent wrong answer must not be able to introduce one. That leaves the guard reliable and platform-independent for the case it was built for, a near-singular mass matrix with no rigid-body modes, and safe but not always effective when rigid-body modes coincide with one. Said plainly in the module docstring, the CHANGELOG and the test that now asserts the portable property there rather than a build-dependent rescue. The retry also gets its own entry point rather than a flag on _solve_dense_general, so the asymmetric production path is provably untouched. --- CHANGELOG.md | 18 +++- src/pybmodes/fem/solver.py | 110 ++++++++++++++++------ tests/fem/test_ill_conditioned_mass.py | 124 +++++++++++++++++++++---- 3 files changed, 204 insertions(+), 48 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2e0ce1d..1d9c7f3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -39,9 +39,21 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ratio of two near-zero quantities and reads ~1 in both candidates however exact each is; on a maximum that floors the alternative and hides a genuinely corrupted elastic mode sitting alongside them, while - per mode it simply registers as no improvement. The retry also - preserves zero eigenvalues, so a spurious trigger on a free-free model - costs a little time rather than deleting a physical mode. + per mode it simply registers as no improvement. The retry preserves + zero and negative eigenvalues and verifies that nothing was dropped + from inside the returned window, so it can never backfill a missing + mode with a higher one and shift the spectrum. + + **Scope.** The rescue is reliable and platform-independent for the case + it was built for, a near-singular mass matrix with no rigid-body modes. + Where rigid-body modes and a near-singular mass matrix coincide it is + safe but not always effective: QZ may return the theoretically real + zero modes as complex-conjugate pairs, and where those land differs + between LAPACK builds. When they fall inside the requested window the + alternative's ordering cannot be verified and the retry declines, + leaving the result no worse than before with `max_residual` still + reporting the problem. Declining is deliberate — a guard added to stop + a silent wrong answer must not be able to introduce one. - `SolverOptions` gains `residual_retry_threshold` and `residual_retry_improvement` for the two conditions above. diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 979c9bb..2b7d662 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -86,6 +86,21 @@ matrices would "win decisively" purely by answering a different question. +**What this guard does not promise.** It rescues the case it was built +for — a near-singular mass matrix, with no rigid-body modes — reliably +and identically on every platform. It is *safe* everywhere else but not +always *effective*: when rigid-body modes and a near-singular mass +matrix coincide, QZ may represent the theoretically real zero modes as +small complex-conjugate pairs that cannot be coerced back to real, and +where those land in the spectrum differs between LAPACK builds. When +they land inside the requested window the alternative's ordering cannot +be verified, so the retry declines and the symmetric result stands. +Declining is the deliberate choice: a guard added to stop a silent wrong +answer must never introduce one, and backfilling a dropped zero mode +with an elastic mode would do exactly that. The result in that situation +is no worse than without the guard, and ``max_residual`` still reports +the problem. + Note on the user-spec mode choice: ``eigsh(..., sigma=0, mode='buckling')`` reduces to ``OP = K^-1 K = I`` for ``sigma=0``, which is degenerate. The standard scipy idiom for "smallest @@ -329,13 +344,18 @@ def solve_modes( # different set — same length, backfilled from higher up — and # the per-index comparison would then be reading two different # spectra against each other. - alt_vals, alt_vecs = _solve_dense_general( - gk_s, gm_s, n_modes, preserve_full_spectrum=True, + alt_vals, alt_vecs, ordering_sound = _general_spectrum_for_retry( + gk_s, gm_s, n_modes, ) _normalize_columns_l2(alt_vecs) alt_r = _modal_residuals(gk_s, gm_s, alt_vals, alt_vecs) - improved = _decisively_improved_modes(sym_r, alt_r, alt_vals.size, - eigvals.size) + improved = ( + _decisively_improved_modes( + sym_r, alt_r, alt_vals.size, eigvals.size, + ) + if ordering_sound + else np.zeros(0, dtype=bool) + ) if improved.any(): idx = int(np.argmax(np.where(improved, sym_r[:improved.size], 0.0))) warnings.warn( @@ -564,36 +584,17 @@ def _solve_dense_symmetric( def _solve_dense_general( gk: np.ndarray, gm: np.ndarray, n_modes: int | None, - *, preserve_full_spectrum: bool = False, ) -> tuple[np.ndarray, np.ndarray]: """Dense LAPACK ``eig`` for genuinely asymmetric problems. Filters eigenvalues to the real, positive, finite subset (matches BModes - JJ's general-matrix path). - - ``preserve_full_spectrum`` drops the **sign** filter, keeping every - real finite eigenvalue including zeros and negatives. Off by default, - so the asymmetric production path keeps the BModes-matching filter it - is validated against. - - The retry path in :func:`solve_modes` turns it on, and needs the - whole spectrum rather than just its zeros. ``eigh`` filters nothing, - so any sign filter here would return a *different set* of modes — - same length, since the gap is backfilled from higher up — and the - per-index comparison would then be reading two different spectra - against each other, accepting a result that had quietly dropped a - mode and shifted every one above it. Both cases are real: a - free-free model's zero-frequency modes, and the genuinely negative - eigenvalues an indefinite ``K`` produces once ``run(gravity=...)`` - loads a column past its buckling weight. With no filter at all both - paths return the ``n_modes`` smallest real eigenvalues, so equal - indices describe the same mode by construction. - """ + JJ's general-matrix path).""" eigvals_all, eigvecs_all = eig(gk, gm) eigvals_real = np.real_if_close(eigvals_all, tol=1000) - valid = np.isreal(eigvals_real) & np.isfinite(eigvals_real.real) - if not preserve_full_spectrum: - valid = valid & (eigvals_real.real > 0.0) - + valid = ( + np.isreal(eigvals_real) + & np.isfinite(eigvals_real.real) + & (eigvals_real.real > 0.0) + ) eigvals = eigvals_real.real[valid] eigvecs = np.real_if_close(eigvecs_all[:, valid], tol=1000).real order = np.argsort(eigvals) @@ -602,6 +603,57 @@ def _solve_dense_general( return eigvals[order], eigvecs[:, order] +def _general_spectrum_for_retry( + gk: np.ndarray, gm: np.ndarray, n_modes: int | None, +) -> tuple[np.ndarray, np.ndarray, bool]: + """The general solve as the retry in :func:`solve_modes` needs it. + + Separate from :func:`_solve_dense_general` so the asymmetric + production path keeps the BModes-matching filter it is validated + against, untouched. + + Two differences, both about making per-index comparison against a + symmetric solve meaningful. + + **No sign filter.** ``eigh`` filters nothing, so discarding + non-positive eigenvalues here would return a *different set* — the + same length, since a truncated request backfills the gap from higher + up — and equal indices would stop meaning equal modes. Both omissions + are reachable: a free-free model's zero-frequency modes, and the + negative eigenvalues an indefinite ``K`` produces once + ``run(gravity=...)`` loads a column past its buckling weight. + + **A verified ordering.** Complex eigenvalues cannot be kept, and on a + symmetric problem ``eig`` does emit a few rounding-induced conjugate + pairs — routinely, and harmlessly, because they land at the stiff end + of the spectrum far above any mode a caller asks for. Demanding that + none appear is therefore too strict to be useful. What actually + matters is narrower: whether anything discarded would have fallen + *inside* the returned window. The third return value reports that, + and the caller declines to swap when it is ``False``, since an + unverifiable ordering is not a basis for replacing a result. + """ + vals_all, vecs_all = eig(gk, gm) + closed = np.real_if_close(vals_all, tol=1000) + keep_mask = np.isreal(closed) & np.isfinite(closed.real) + + vals = closed.real[keep_mask] + vecs = np.real_if_close(vecs_all[:, keep_mask], tol=1000).real + order = np.argsort(vals) + vals, vecs = vals[order], vecs[:, order] + + keep = vals.size if n_modes is None else min(n_modes, vals.size) + dropped = vals_all[~keep_mask] + ordering_sound = True + if dropped.size and keep: + # Sound exactly when every discarded eigenvalue sits above the + # window, so the window really is the smallest ``keep`` modes. + ordering_sound = bool( + np.min(dropped.real) > vals[keep - 1] + ) + return vals[:keep], vecs[:, :keep], ordering_sound + + # --------------------------------------------------------------------------- # Utility helpers # --------------------------------------------------------------------------- diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 5663858..cec0539 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -270,11 +270,14 @@ def test_the_retry_preserves_zero_eigenvalues(self): """The property that makes a false positive harmless: the alternative solve keeps the rigid modes rather than filtering them, so both candidates describe the same spectrum.""" - from pybmodes.fem.solver import _solve_dense_general + from pybmodes.fem.solver import ( + _general_spectrum_for_retry, + _solve_dense_general, + ) gk, gm = self._free_free_with_a_zero_mode() dropped, _v = _solve_dense_general(gk, gm, 6) - kept, _w = _solve_dense_general(gk, gm, 6, preserve_full_spectrum=True) + kept, _w, _ok = _general_spectrum_for_retry(gk, gm, 6) assert abs(kept[0]) < 1.0e-8 * float(np.max(np.abs(kept))) assert np.min(dropped) > 0.0 # The default filter loses the zero mode and shifts the rest up. @@ -385,24 +388,57 @@ def test_a_shorter_alternative_is_never_an_improvement(self): alt_r = np.array([1.0e-9, 1.0e-9]) assert not _decisively_improved_modes(sym_r, alt_r, 2, 3).any() - def test_the_breakdown_is_caught_and_corrected(self): + def test_the_result_is_never_made_worse(self): + """The portable guarantee when rigid modes and an ill-conditioned + mass matrix coincide. + + Whether the rescue can *fire* here is build-dependent, and + deliberately so. QZ may represent this problem's theoretically + real zero modes as small complex-conjugate pairs that + ``real_if_close`` will not coerce; where those land in the + spectrum differs between LAPACK builds. On one they sit at the + stiff end, far above anything a caller asks for, and the retry + proceeds. On another they are the zero modes themselves, and + keeping the alternative would mean backfilling the gap with + elastic modes and silently shifting the spectrum. + + So the guard verifies its ordering and declines when it cannot, + and the property asserted here is the one that holds either way: + the requested modes come back, in ascending order, and nothing is + lost. ``TestIllConditionedMassIsCaught`` pins the rescue itself on + the case it was built for, which has no rigid modes and behaves + identically everywhere. + """ gk, gm = self._rigid_plus_ill_conditioned() - with pytest.warns(RuntimeWarning, match="do not satisfy"): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") eigvals, _v, diag = solve_modes( gk, gm, n_modes=10, return_diagnostics=True, ) - assert diag.residual_fallback is True assert eigvals.size == 10 + assert diag.n_returned == 10 + assert np.all(np.diff(eigvals) >= -1.0e-12 * max(1.0, abs(eigvals).max())) - def test_the_corrected_spectrum_matches_the_underlying_one(self): - """The rotation and the extra rigid DOFs do not change the - cantilever's own eigenvalues, so the corrected solve must still - contain the analytic lump frequency.""" - gk, gm = self._rigid_plus_ill_conditioned() - with pytest.warns(RuntimeWarning): - eigvals, _v = solve_modes(gk, gm, n_modes=10) - f = eigvals_to_hz(eigvals, ROMG) - assert np.min(np.abs(f - _analytic())) < 5.0e-3 * _analytic() + def test_an_unverifiable_ordering_declines_the_swap(self, monkeypatch): + """When the alternative drops a mode from inside the window, its + ordering cannot be trusted and the symmetric result stands.""" + import pybmodes.fem.solver as solvermod + + real_fn = solvermod._general_spectrum_for_retry + + def unsound(gk, gm, n_modes): + vals, vecs, _ok = real_fn(gk, gm, n_modes) + return vals, vecs, False + + monkeypatch.setattr(solvermod, "_general_spectrum_for_retry", unsound) + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with warnings.catch_warnings(): + warnings.simplefilter("error") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" class TestNegativeEigenvaluesSurviveTheRetry: @@ -434,11 +470,14 @@ def _indefinite(self): return 0.5 * (k_rot + k_rot.T), 0.5 * (m_rot + m_rot.T) def test_the_default_filter_would_drop_the_unstable_mode(self): - from pybmodes.fem.solver import _solve_dense_general + from pybmodes.fem.solver import ( + _general_spectrum_for_retry, + _solve_dense_general, + ) gk, gm = self._indefinite() dropped, _v = _solve_dense_general(gk, gm, 6) - kept, _w = _solve_dense_general(gk, gm, 6, preserve_full_spectrum=True) + kept, _w, _ok = _general_spectrum_for_retry(gk, gm, 6) assert kept.min() < 0.0 # the unstable mode is present assert dropped.min() > 0.0 # and absent from the default # Same length, shifted spectrum — the trap the index comparison @@ -542,6 +581,59 @@ def test_measuring_against_the_unsymmetrised_pair_would_have_tripped(self): assert against_raw.max() > 0.1 +class TestTheRetryVerifiesItsOrderingRatherThanAssumingIt: + """Index matching is only sound if the alternative recovered every + mode. + + A symmetric problem has ``ngd`` real eigenvalues, so if ``eig`` + returns fewer after its real / finite filter, one was discarded — and + a truncated request would have backfilled the gap from higher up, + leaving equal indices pointing at different modes. The sign filter was + one route to that and is switched off; a rounding-induced complex pair + is another that no flag can prevent, so completeness is checked rather + than assumed. + """ + + def test_a_sound_ordering_is_the_precondition_for_swapping(self): + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with pytest.warns(RuntimeWarning): + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + # This case does swap, so the precondition held. + assert diag.residual_fallback is True + + def test_a_drop_above_the_window_leaves_the_ordering_sound(self): + """Complex pairs at the stiff end are routine and harmless — they + sit far above anything the caller asked for, so the window is + still the smallest ``n_modes``.""" + from pybmodes.fem.solver import _general_spectrum_for_retry + + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + vals, _v, sound = _general_spectrum_for_retry(gk, gm, 4) + assert sound is True + assert vals.size == 4 + + def test_a_drop_inside_the_window_makes_the_ordering_unsound(self): + """A discarded eigenvalue below the top of the window means the + window is not the smallest ``n_modes`` and cannot be compared by + index.""" + from pybmodes.fem.solver import _general_spectrum_for_retry + + # Two well-separated real modes plus a complex pair placed below + # them, which no coercion will make real. + gm = np.eye(4) + gk = np.array([ + [1.0, 0.0, 0.0, 0.0], + [0.0, 2.0, 0.0, 0.0], + [0.0, 0.0, 0.0, -1.0e-3], + [0.0, 0.0, 1.0e-3, 0.0], + ]) + vals, _v, sound = _general_spectrum_for_retry(gk, gm, 2) + assert sound is False + assert vals.size == 2 + + class TestDiagnosticsContract: def test_residual_fallback_defaults_to_false(self): gk, gm = _cantilever_with_tip_lump(13, REALISTIC) From d06ef0bd84b375fb2ac0197fd1b0a78be19c4630 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 18:06:24 +0900 Subject: [PATCH 08/28] test: build the sound-ordering case rather than hope for it The drop-above-the-window test was reading a real matrix and trusting that its LAPACK build would emit the complex pair at the stiff end. Which end those land on is exactly what differs between builds, and is what CI caught me on twice already. An explicit rotation block places the pair where the test needs it and says the same thing everywhere. --- tests/fem/test_ill_conditioned_mass.py | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index cec0539..79fdbef 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -606,13 +606,28 @@ def test_a_sound_ordering_is_the_precondition_for_swapping(self): def test_a_drop_above_the_window_leaves_the_ordering_sound(self): """Complex pairs at the stiff end are routine and harmless — they sit far above anything the caller asked for, so the window is - still the smallest ``n_modes``.""" + still the smallest ``n_modes``. + + Built from an explicit rotation block rather than by hoping a + real matrix produces one, so it says the same thing on every + LAPACK build. ``[[a, -b], [b, a]]`` has eigenvalues ``a ± bi``; + with ``a`` far above the window and ``b`` far too large for + ``real_if_close`` to coerce, the pair is dropped from well above + the modes being compared. + """ from pybmodes.fem.solver import _general_spectrum_for_retry - gk, gm = _cantilever_with_tip_lump(27, LIGHT) - vals, _v, sound = _general_spectrum_for_retry(gk, gm, 4) + gm = np.eye(4) + gk = np.array([ + [1.0, 0.0, 0.0, 0.0], + [0.0, 2.0, 0.0, 0.0], + [0.0, 0.0, 1.0e6, -1.0e3], + [0.0, 0.0, 1.0e3, 1.0e6], + ]) + vals, _v, sound = _general_spectrum_for_retry(gk, gm, 2) assert sound is True - assert vals.size == 4 + assert vals.size == 2 + assert np.allclose(vals, [1.0, 2.0]) def test_a_drop_inside_the_window_makes_the_ordering_unsound(self): """A discarded eigenvalue below the top of the window means the From fde9191976c3b122ce18bbccf92d095ba91926a7 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 18:13:35 +0900 Subject: [PATCH 09/28] fix: scope the retry to the dense path, and never let it raise Two from Codex on #140. P1, the sparse spectrum. eigsh(sigma=0, which="LM") selects the modes nearest zero in magnitude while the retry selects the algebraically smallest, so with negative eigenvalues present, a post-buckling gravity column, the two windows are different sets and a per-index comparison pairs unrelated modes. Rather than reconcile the two selections I removed the case. Only the dense symmetric path is retried now, which is a statement about which matrix each routine factorises rather than a convenience: eigh reduces through a Cholesky factor of the mass matrix, the failure this guard exists for, while eigsh factorises K and is unaffected. The mesh sweep that motivated the whole change already showed that, returning correct frequencies on exactly the meshes large enough to take the sparse path. So the mismatch is gone by construction and nothing is lost. P2, the retry could raise. A pencil defective enough to break the symmetric reduction can also break eig, and an unguarded LinAlgError would have made this guard destroy usable results on precisely the inputs it was added to help. It now declines and keeps the symmetric result, logging why. --- src/pybmodes/fem/solver.py | 62 ++++++++++++++------ tests/fem/test_ill_conditioned_mass.py | 81 ++++++++++++++++++++++++++ 2 files changed, 126 insertions(+), 17 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 2b7d662..8f9283b 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -36,12 +36,14 @@ the retry path when a symmetric solve comes back with a large backward error — see below. -Both symmetric paths reduce ``K x = λ M x`` through a Cholesky factor -of the mass matrix, and that reduction degrades once ``M`` is nearly -singular, which a very light beam carrying a very heavy lump produces. -The failure mode is silent: LAPACK returns confidently wrong low modes -rather than raising. :func:`solve_modes` therefore checks the backward -error of every symmetric solve and, when it exceeds +The **dense** symmetric path reduces ``K x = λ M x`` through a Cholesky +factor of the mass matrix, and that reduction degrades once ``M`` is +nearly singular, which a very light beam carrying a very heavy lump +produces. The failure mode is silent: LAPACK returns confidently wrong +low modes rather than raising. (The sparse path factorises ``K`` +instead, so it is unaffected and is not retried.) +:func:`solve_modes` therefore checks the backward error of a dense +symmetric solve and, when it exceeds :attr:`~pybmodes.options.SolverOptions.residual_retry_threshold`, tries the general path as well — taking its result only if it is better by :attr:`~pybmodes.options.SolverOptions.residual_retry_improvement`, and @@ -325,7 +327,23 @@ def solve_modes( res_k, res_m = (0.5 * (gk + gk.T), 0.5 * (gm + gm.T)) if sym else (gk, gm) residual_fallback = False - if sym: + # Only the *dense* symmetric path is retried, and that is a statement + # about which matrix each routine factorises rather than a + # convenience. ``eigh`` reduces through a Cholesky factor of the mass + # matrix, which is the one this guard exists for. ``eigsh(sigma=0, + # mode='normal')`` factorises ``K`` instead, so a near-singular ``M`` + # does not degrade it — the mesh sweep that motivated this work + # returns correct frequencies on exactly the meshes large enough to + # take the sparse path. + # + # Excluding it also removes a mismatch that would otherwise need + # separate handling: ``which="LM"`` on ``OP = K^-1 M`` selects the + # modes nearest zero *in magnitude*, while the retry selects the + # algebraically smallest. With negative eigenvalues present — a + # post-buckling ``run(gravity=...)`` column — those are different + # sets, and a per-index comparison between them would be pairing + # unrelated modes. + if sym and path == "dense_symmetric": # Measure — and retry — against the matrices the symmetric paths # actually solved. Both symmetrise internally, and the accepted # skew is only guaranteed small relative to ``max|K|``: in a model @@ -338,16 +356,26 @@ def solve_modes( gk_s, gm_s = res_k, res_m sym_r = _modal_residuals(gk_s, gm_s, eigvals, eigvecs) if sym_r.size and float(sym_r.max()) > _SOLVER_OPTIONS.residual_retry_threshold: - # ``preserve_full_spectrum`` so the two candidates describe the - # same spectrum and equal indices mean the same mode. ``eigh`` - # filters nothing, so any sign filter here would return a - # different set — same length, backfilled from higher up — and - # the per-index comparison would then be reading two different - # spectra against each other. - alt_vals, alt_vecs, ordering_sound = _general_spectrum_for_retry( - gk_s, gm_s, n_modes, - ) - _normalize_columns_l2(alt_vecs) + try: + alt_vals, alt_vecs, ordering_sound = _general_spectrum_for_retry( + gk_s, gm_s, n_modes, + ) + except (np.linalg.LinAlgError, ValueError) as exc: + # The alternative is a best-effort second opinion, not a + # requirement. A pencil defective enough to break the + # symmetric reduction can also break ``eig``, and turning + # that into a hard failure would make this guard destroy + # usable results on exactly the inputs it was added to + # help. Decline and keep what we have. + _log.warning( + "solve_modes: residual retry failed (%r); keeping the " + "symmetric result", exc, + ) + alt_vals = np.empty(0) + alt_vecs = np.empty((eigvecs.shape[0], 0)) + ordering_sound = False + if alt_vecs.size: + _normalize_columns_l2(alt_vecs) alt_r = _modal_residuals(gk_s, gm_s, alt_vals, alt_vecs) improved = ( _decisively_improved_modes( diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 79fdbef..37a59b4 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -649,6 +649,87 @@ def test_a_drop_inside_the_window_makes_the_ordering_unsound(self): assert vals.size == 2 +class TestTheRetryIsScopedToTheDensePath: + """Only the dense symmetric path is retried, and that is about which + matrix each routine factorises. + + ``eigh`` reduces through a Cholesky factor of the **mass** matrix, + which is the failure this guard exists for. ``eigsh(sigma=0, + mode='normal')`` factorises ``K`` instead and is unaffected — the + mesh sweep that motivated the work returns correct frequencies on + exactly the meshes large enough to take the sparse path. + + It also removes a spectrum mismatch by construction: ``which="LM"`` + selects the modes nearest zero *in magnitude* while the retry selects + the algebraically smallest, and with negative eigenvalues present + those are different sets that must never be compared by index. + """ + + def test_the_sparse_path_is_left_alone(self, monkeypatch): + import pybmodes.fem.solver as solvermod + + # Force the sparse path on a small problem, then make the + # residual look terrible. The retry must still not run. + gk, gm = _cantilever_with_tip_lump(13, LIGHT) + monkeypatch.setattr(solvermod, "_SPARSE_NDOF_THRESHOLD", 1) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.path == "sparse_shift_invert" + assert diag.residual_fallback is False + + def test_the_sparse_path_gets_the_ill_conditioned_case_right(self): + """Why leaving it alone is safe rather than a gap: it factorises + K, so a near-singular M does not degrade it.""" + gk, gm = _cantilever_with_tip_lump(101, LIGHT) + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.path == "sparse_shift_invert" + f = float(eigvals_to_hz(eigvals, ROMG)[0]) + assert f == pytest.approx(_analytic(), rel=5.0e-3) + + +class TestARetryFailureIsNotAHardFailure: + """A pencil defective enough to break the symmetric reduction can + also break ``eig``. Turning that into an exception would make the + guard destroy usable results on exactly the inputs it exists for.""" + + def test_a_raising_retry_keeps_the_symmetric_result(self, monkeypatch): + import pybmodes.fem.solver as solvermod + + def boom(gk, gm, n_modes): + raise np.linalg.LinAlgError("did not converge") + + monkeypatch.setattr(solvermod, "_general_spectrum_for_retry", boom) + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" + assert eigvals.size == 4 + + def test_a_value_error_is_handled_the_same_way(self, monkeypatch): + import pybmodes.fem.solver as solvermod + + def boom(gk, gm, n_modes): + raise ValueError("array must not contain infs or NaNs") + + monkeypatch.setattr(solvermod, "_general_spectrum_for_retry", boom) + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + + class TestDiagnosticsContract: def test_residual_fallback_defaults_to_false(self): gk, gm = _cantilever_with_tip_lump(13, REALISTIC) From 603eae8b50d700c1050ac328fff256a5ace0bc80 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 18:23:12 +0900 Subject: [PATCH 10/28] fix: bound the retry, and keep two warnings honest about their scope A self-review pass plus two from Codex on #140. Found while auditing: the retry was reachable at any size. A sparse solve that fails to converge falls back to dense however large the system, and an unbounded dense eig there could spend minutes on a result already in hand. A guard against a silent wrong answer should not be able to turn one into a silent hang, so it now declines above residual_retry_max_ndof and leaves the backward error to the diagnostics. Codex, mode-count warning. After a swap the path label reads dense_general, which made the shortfall warning claim modes had been filtered as complex or non-positive. The retry preserves the whole spectrum, so nothing was filtered, and a 117-DOF system asked for 1000 modes was reported as defective for returning its 117. The warning now keys on modes actually discarded rather than on the label, which also fixes the same spurious report on the pre-existing general path. Codex, CHANGELOG. The note still claimed every symmetric solve is retried and that eigsh factorises the mass matrix, both of which my own scoping change had made false. Rewritten to say the retry is dense-only, that sparse results never set residual_fallback, and why. Also states in the diagnostics that symmetric=True alongside path=dense_general is not a contradiction: the first describes the input, the second which routine ran, and residual_fallback separates that case from a genuinely asymmetric solve. --- CHANGELOG.md | 44 +++++++---- src/pybmodes/fem/solver.py | 54 ++++++++----- src/pybmodes/options.py | 12 +++ tests/fem/test_ill_conditioned_mass.py | 100 +++++++++++++++++++++++++ 4 files changed, 178 insertions(+), 32 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1d9c7f3..2cfce89 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,23 +10,31 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ### Fixed -- **The symmetric eigensolvers could return confidently wrong low modes - on a near-singular mass matrix, silently.** Both `scipy.linalg.eigh` - and `scipy.sparse.linalg.eigsh` reduce `K x = λ M x` through a Cholesky - factor of the mass matrix, and that reduction loses accuracy when a - very light beam carries a very heavy lump. LAPACK does not raise there - — it returns wrong frequencies. On a 100 m cantilever with a 4000:1 - lump-to-beam mass ratio the reported fundamental was 0.103 Hz against a - true 0.0436 Hz, a factor of 2.4, and the answer wandered - non-monotonically with mesh density. +- **The dense symmetric eigensolver could return confidently wrong low + modes on a near-singular mass matrix, silently.** `scipy.linalg.eigh` + reduces `K x = λ M x` through a Cholesky factor of the **mass** matrix, + and that reduction loses accuracy when a very light beam carries a very + heavy lump. LAPACK does not raise there — it returns wrong frequencies. + On a 100 m cantilever with a 4000:1 lump-to-beam mass ratio the + reported fundamental was 0.103 Hz against a true 0.0436 Hz, a factor of + 2.4, and the answer wandered non-monotonically with mesh density. `solve_modes` now checks the backward error `||K x - λ M x|| / ||K x||` - of every symmetric solve and, when it is large, redoes it through the - general dense path, which factorises neither matrix. The retried result - is taken only when it is better by an order of magnitude, and a + of a **dense** symmetric solve and, when it is large, redoes it through + the general dense path, which factorises neither matrix. The retried + result is taken only when it is better by an order of magnitude, and a `RuntimeWarning` names the swap. `SolverDiagnostics` gains `residual_fallback` recording it. + The **sparse** path is deliberately not retried and never sets + `residual_fallback`. `eigsh(sigma=0, mode='normal')` factorises `K` + rather than the mass matrix, so a near-singular `M` does not degrade + it — on the mesh sweep that motivated this work it returned correct + frequencies on exactly the meshes large enough to select it. Retrying + it would also mean comparing two different mode sets, since + `which="LM"` selects the modes nearest zero in magnitude while the + retry selects the algebraically smallest. + **No existing result changes.** The decisive-improvement condition is what guarantees that: a real deck can carry a large backward error without being broken, and on the bundled NREL 5MW land tower (whose @@ -55,8 +63,16 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). reporting the problem. Declining is deliberate — a guard added to stop a silent wrong answer must not be able to introduce one. -- `SolverOptions` gains `residual_retry_threshold` and - `residual_retry_improvement` for the two conditions above. + The retry is also bounded in size, since a sparse solve that fails to + converge falls back to dense at any size and an unbounded `eig` there + could take minutes on a result already in hand. And it can decline: if + the alternative solver raises on the same defective pencil, the + symmetric result and its diagnostics are kept rather than the whole + solve failing. + +- `SolverOptions` gains `residual_retry_threshold`, + `residual_retry_improvement` and `residual_retry_max_ndof` for the + conditions above. ## [1.18.0] — 2026-08-12 diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 8f9283b..14b2060 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -151,9 +151,16 @@ class SolverDiagnostics: path : which solver path produced the result. One of ``"sparse_shift_invert"``, ``"dense_symmetric"``, ``"dense_general"``. - symmetric : whether the assembled matrices were treated as symmetric - (``eigh`` / sparse) rather than routed through the general - ``eig`` path. + symmetric : whether the assembled matrices were **classified** as + symmetric, i.e. whether their asymmetry was within + :attr:`~pybmodes.options.SolverOptions.symmetry_rtol`. This is a + property of the input, not a record of which routine ran, so a + residual retry leaves it ``True`` while moving ``path`` to + ``"dense_general"``. That pairing is not a contradiction: the + matrices were symmetric, and the general routine was used on their + symmetrised form because the symmetric one had failed on it. + ``residual_fallback`` is what distinguishes that case from a + genuinely asymmetric solve. n_requested : modes asked for (``None`` means the full spectrum). n_returned : modes actually returned. Fewer than ``n_requested`` means the general path filtered out complex / non-positive @@ -343,7 +350,15 @@ def solve_modes( # post-buckling ``run(gravity=...)`` column — those are different # sets, and a per-index comparison between them would be pairing # unrelated modes. - if sym and path == "dense_symmetric": + # + # The size ceiling matters only because a sparse solve that fails to + # converge falls back to the dense path at *any* size, where an + # unbounded ``eig`` could spend minutes on a result already in hand. + if ( + sym + and path == "dense_symmetric" + and ngd <= _SOLVER_OPTIONS.residual_retry_max_ndof + ): # Measure — and retry — against the matrices the symmetric paths # actually solved. Both symmetrise internally, and the accepted # skew is only guaranteed small relative to ``max|K|``: in a model @@ -410,25 +425,28 @@ def solve_modes( # Surface that rather than letting it pass silently (a downstream # broadcast would otherwise fail with an opaque shape error). # - # Gate the warning to the general path only (Codex P2). The dense - # symmetric path also returns fewer than ``n_modes`` when the request - # simply exceeds the available DOFs (it truncates to - # ``min(n_modes, ngd)``), which is a benign "asked for more modes than - # the system has" case, not a defective eigenproblem — warning there - # would mislead, and would fail callers that treat warnings as errors. + # Gate on modes actually *discarded*, not on the path label. Two + # benign shortfalls would otherwise be reported as a defective + # eigenproblem. Asking for more modes than the system has is one: + # every path truncates to ``min(n_modes, ngd)``, which is a request + # the caller can reasonably make. A residual retry is the other — it + # relabels the path ``"dense_general"`` while preserving the whole + # spectrum, so nothing was filtered, and a 117-DOF system asked for + # 1000 modes would be reported as defective for returning its 117. n_returned = int(eigvecs.shape[1]) + n_available = ngd if n_modes is None else min(n_modes, ngd) if ( path == "dense_general" - and n_modes is not None - and n_returned < n_modes + and not residual_fallback + and n_returned < n_available ): warnings.warn( - f"solve_modes recovered only {n_returned} of the requested " - f"{n_modes} modes via the general (non-symmetric) eig path. " - f"The eigenproblem is likely near-degenerate or defective (a " - f"non-symmetric PlatformSupport block can do this); the " - f"missing modes had complex or non-positive eigenvalues and " - f"were filtered out.", + f"solve_modes recovered only {n_returned} of the " + f"{n_available} modes available via the general " + f"(non-symmetric) eig path. The eigenproblem is likely " + f"near-degenerate or defective (a non-symmetric " + f"PlatformSupport block can do this); the missing modes had " + f"complex or non-positive eigenvalues and were filtered out.", RuntimeWarning, stacklevel=2, ) diff --git a/src/pybmodes/options.py b/src/pybmodes/options.py index ea981ef..d738c0a 100644 --- a/src/pybmodes/options.py +++ b/src/pybmodes/options.py @@ -85,6 +85,17 @@ class SolverOptions: ~4e10) reaches ~2e-2 and is *not* meant to trigger; the degraded regime starts around 0.7. The default splits the last two gaps with roughly 5x either side. + residual_retry_max_ndof : int, default 2000 + Largest reduced system the retry will attempt. The retry is a + dense ``eig``, whose cost grows as ``ngd^3`` with a much larger + constant than the ``eigh`` it is checking. Normally that is + bounded by ``sparse_ndof_threshold``, since anything bigger takes + the sparse path and is not retried — but a sparse solve that + *fails to converge* falls back to dense at any size, and there an + unbounded retry could spend minutes on a model whose result was + already available. A guard against a silent wrong answer should + not be able to turn one into a silent hang, so above this size it + declines and leaves the backward error to the diagnostics. residual_retry_improvement : float, default 0.1 How much better the general path's backward error must be before its result is taken. The second, and more important, guard: on @@ -100,6 +111,7 @@ class SolverOptions: symmetry_rtol: float = 1.0e-12 residual_retry_threshold: float = 0.1 residual_retry_improvement: float = 0.1 + residual_retry_max_ndof: int = 2000 @dataclass(frozen=True) diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 37a59b4..1530eb6 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -692,6 +692,63 @@ def test_the_sparse_path_gets_the_ill_conditioned_case_right(self): assert f == pytest.approx(_analytic(), rel=5.0e-3) +class TestTheRetryCostIsBounded: + """A guard against a silent wrong answer must not be able to turn one + into a silent hang. + + The retry is a dense ``eig``, whose cost grows as ``ngd^3`` with a + much larger constant than the ``eigh`` it checks. Normally that is + bounded by the sparse dispatch threshold, but a sparse solve that + *fails to converge* falls back to dense at any size. + """ + + def test_a_system_above_the_ceiling_is_not_retried(self, monkeypatch): + import dataclasses + + import pybmodes.fem.solver as solvermod + + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + ngd = gk.shape[0] + monkeypatch.setattr( + solvermod, "_SOLVER_OPTIONS", + dataclasses.replace( + solvermod._SOLVER_OPTIONS, residual_retry_max_ndof=ngd - 1, + ), + ) + called = [] + real = solvermod._general_spectrum_for_retry + monkeypatch.setattr( + solvermod, "_general_spectrum_for_retry", + lambda *a, **k: (called.append(1), real(*a, **k))[1], + ) + with warnings.catch_warnings(): + warnings.simplefilter("error") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert called == [] + assert diag.residual_fallback is False + + def test_at_the_ceiling_it_still_runs(self, monkeypatch): + import dataclasses + + import pybmodes.fem.solver as solvermod + + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + monkeypatch.setattr( + solvermod, "_SOLVER_OPTIONS", + dataclasses.replace( + solvermod._SOLVER_OPTIONS, + residual_retry_max_ndof=gk.shape[0], + ), + ) + with pytest.warns(RuntimeWarning): + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is True + + class TestARetryFailureIsNotAHardFailure: """A pencil defective enough to break the symmetric reduction can also break ``eig``. Turning that into an exception would make the @@ -730,6 +787,49 @@ def boom(gk, gm, n_modes): assert diag.residual_fallback is False +class TestTheModeCountWarningStaysHonest: + """A shortfall is only newsworthy when modes were actually discarded. + + Two benign cases would otherwise be reported as a defective + eigenproblem: asking for more modes than the system has, which every + path truncates, and a residual retry, which relabels the path + ``"dense_general"`` while preserving the whole spectrum. + """ + + def test_an_overlarge_request_after_a_retry_is_not_reported(self): + gk, gm = _cantilever_with_tip_lump(13, LIGHT) + ngd = gk.shape[0] + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=ngd + 500, return_diagnostics=True, + ) + assert diag.residual_fallback is True + assert eigvals.size == ngd + assert not any( + "recovered only" in str(w.message) for w in caught + ), [str(w.message) for w in caught] + + def test_an_overlarge_request_without_a_retry_is_not_reported(self): + gk, gm = _cantilever_with_tip_lump(13, REALISTIC) + ngd = gk.shape[0] + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + eigvals, _v = solve_modes(gk, gm, n_modes=ngd + 500) + assert eigvals.size == ngd + assert not any("recovered only" in str(w.message) for w in caught) + + def test_a_genuine_shortfall_is_still_reported(self): + """An asymmetric pencil whose modes really are filtered away must + still say so.""" + # Antisymmetric K: every eigenvalue is imaginary, so the general + # path recovers none of them. + gk = np.array([[0.0, -1.0], [1.0, 0.0]]) + gm = np.eye(2) + with pytest.warns(RuntimeWarning, match="recovered only"): + solve_modes(gk, gm, n_modes=2) + + class TestDiagnosticsContract: def test_residual_fallback_defaults_to_false(self): gk, gm = _cantilever_with_tip_lump(13, REALISTIC) From dc26e04879f0a900b6d63218d076d9ec1d1be9d6 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 18:32:19 +0900 Subject: [PATCH 11/28] fix: refuse a retry that trades one mode for another (#140) Codex P1. Acceptance keyed on "any mode improved decisively", but taking the retry replaces the whole spectrum, not the modes that prompted it. So a candidate that rescued mode 0 while pushing a previously acceptable mode 3 above the failure threshold was accepted, handing back a new bad mode in place of an old one. The comparison now returns both verdicts and the caller requires improvement somewhere and regression nowhere. Regression deliberately mirrors improvement rather than introducing a second notion of acceptable: a mode has regressed when it ends up above the threshold and is worse there by the margin that would have counted as decisive the other way. The symmetry is what keeps rigid-body modes out of it, since their residual reads ~1 in both candidates and wobbles either way, and a bare alt > sym test would read that as a regression and block every rescue sitting beside a free-free mode. Tests take the exact residual vectors from the report, plus the two cases the rule has to keep apart: a mode six times worse but nowhere near the threshold is not a regression, and neither is rigid-body noise. --- CHANGELOG.md | 9 +++ src/pybmodes/fem/solver.py | 54 ++++++++++++----- tests/fem/test_ill_conditioned_mass.py | 83 ++++++++++++++++++++++++-- 3 files changed, 126 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2cfce89..590bbc0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,6 +42,15 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). 1.4× better while *splitting* a degenerate fore-aft / side-side pair the symmetric solver resolves exactly. + Acceptance requires the candidate to be **non-regressive** as well as + decisively better somewhere. Taking the retry replaces the whole + spectrum, not the modes that prompted it, so a candidate that rescues + one mode while pushing a previously acceptable one above the failure + threshold is refused — it would hand back a new bad mode in place of an + old one. Regression mirrors improvement with the same threshold and + factor, so rigid-body modes reading ~1 in both candidates are not + mistaken for it. + The comparison is made **per mode** rather than on the two maxima, so that rigid-body modes cannot distort it. Their backward error is a ratio of two near-zero quantities and reads ~1 in both candidates diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 14b2060..790e990 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -392,14 +392,17 @@ def solve_modes( if alt_vecs.size: _normalize_columns_l2(alt_vecs) alt_r = _modal_residuals(gk_s, gm_s, alt_vals, alt_vecs) - improved = ( - _decisively_improved_modes( + improved, regressed = ( + _compare_candidate_modes( sym_r, alt_r, alt_vals.size, eigvals.size, ) if ordering_sound - else np.zeros(0, dtype=bool) + else (np.zeros(0, dtype=bool), np.zeros(0, dtype=bool)) ) - if improved.any(): + # Accepting replaces the whole spectrum, not just the modes + # that prompted the retry, so a candidate that fixes one mode + # while ruining another is not an improvement to the result. + if improved.any() and not regressed.any(): idx = int(np.argmax(np.where(improved, sym_r[:improved.size], 0.0))) warnings.warn( f"the symmetric eigensolver returned " @@ -497,13 +500,29 @@ def _build_diagnostics( # one is a rigid-body mode: a free-free floating platform has up to six, # and an unrestrained DOF (a symmetric column's yaw) gives an exactly # zero one. -def _decisively_improved_modes( +def _compare_candidate_modes( sym_r: np.ndarray, alt_r: np.ndarray, n_alt: int, n_sym: int, -) -> np.ndarray: - """Which modes the general path solves decisively better, per mode. +) -> tuple[np.ndarray, np.ndarray]: + """Per-mode verdicts on the alternative: ``(improved, regressed)``. + + Both are needed because accepting the retry replaces the **whole** + spectrum, not the modes that prompted it. A candidate that fixes one + mode while ruining another is not an improvement to the result even + though it is an improvement to that mode, so "some mode got + decisively better" is only half the test; the other half is that no + mode got decisively worse. + + Regression mirrors improvement exactly, using the same threshold and + the same factor rather than a second notion of acceptable: a mode has + regressed when it ends up above the failure threshold *and* is worse + there by the margin that would have counted as decisive in the other + direction. The symmetry matters for rigid-body modes, whose residual + is ~1 in both candidates and wobbles a little either way — that is + noise, not a regression, and a bare ``alt > sym`` test would read it + as one and block every rescue that happens to sit beside them. The comparison has to be **per mode**, not on the two maxima. A rigid-body mode's backward error is a ratio of two near-zero @@ -522,19 +541,22 @@ def _decisively_improved_modes( retry preserves rigid-body modes for this reason), so equal indices describe the same mode. - Returns a boolean mask over the compared modes. Empty when the - alternative recovered fewer modes than the symmetric solve — losing a - mode is never an improvement, whatever the residuals say. + Returns two boolean masks over the compared modes, both empty when + the alternative recovered fewer modes than the symmetric solve — + losing a mode is never an improvement, whatever the residuals say. """ + empty = np.zeros(0, dtype=bool) if n_alt < n_sym: - return np.zeros(0, dtype=bool) + return empty, empty n = min(sym_r.size, alt_r.size) if n == 0: - return np.zeros(0, dtype=bool) - return ( - (sym_r[:n] > _SOLVER_OPTIONS.residual_retry_threshold) - & (alt_r[:n] < _SOLVER_OPTIONS.residual_retry_improvement * sym_r[:n]) - ) + return empty, empty + threshold = _SOLVER_OPTIONS.residual_retry_threshold + factor = _SOLVER_OPTIONS.residual_retry_improvement + sym, alt = sym_r[:n], alt_r[:n] + improved = (sym > threshold) & (alt < factor * sym) + regressed = (alt > threshold) & (sym < factor * alt) + return improved, regressed def _modal_residuals( diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 1530eb6..03664db 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -367,7 +367,7 @@ def test_the_maxima_rule_misses_what_the_per_mode_rule_catches(self): two rigid-body modes reading ~1 in both candidates, one elastic mode corrupted to 0.8 and fixed to 1e-9, one mode already exact. """ - from pybmodes.fem.solver import _decisively_improved_modes + from pybmodes.fem.solver import _compare_candidate_modes sym_r = np.array([1.0, 1.0, 0.8, 1.0e-12]) alt_r = np.array([1.0, 1.0, 1.0e-9, 1.0e-12]) @@ -378,15 +378,18 @@ def test_the_maxima_rule_misses_what_the_per_mode_rule_catches(self): # Per mode, the corrupted one is unmissable and the rigid ones # register as exactly what they are: no improvement either way. - improved = _decisively_improved_modes(sym_r, alt_r, 4, 4) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 4, 4) assert improved.tolist() == [False, False, True, False] + assert not regressed.any() def test_a_shorter_alternative_is_never_an_improvement(self): - from pybmodes.fem.solver import _decisively_improved_modes + from pybmodes.fem.solver import _compare_candidate_modes sym_r = np.array([1.0, 0.8, 1.0e-12]) alt_r = np.array([1.0e-9, 1.0e-9]) - assert not _decisively_improved_modes(sym_r, alt_r, 2, 3).any() + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 3) + assert not improved.any() + assert not regressed.any() def test_the_result_is_never_made_worse(self): """The portable guarantee when rigid modes and an ill-conditioned @@ -787,6 +790,78 @@ def boom(gk, gm, n_modes): assert diag.residual_fallback is False +class TestARetryThatTradesModesIsRefused: + """Accepting replaces the whole spectrum, so a candidate that fixes + one mode while ruining another is not an improvement to the result. + + The failing shape: the general solve is exact on the modes that + prompted the retry but pushes a previously acceptable mode above the + failure threshold. Judging on "any mode improved" would take it and + hand back a spectrum with a new bad mode in place of an old one. + """ + + def test_a_trade_is_not_an_improvement(self): + from pybmodes.fem.solver import _compare_candidate_modes + + # Mode 0 rescued decisively; mode 3 was fine and is now above the + # threshold and an order worse. + sym_r = np.array([0.52, 1.1e-3, 5.7e-6, 0.017]) + alt_r = np.array([1.1e-10, 3.0e-11, 3.5e-5, 0.26]) + + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 4, 4) + assert improved[0] + assert regressed[3] + # The rule the caller applies. + assert not (improved.any() and not regressed.any()) + + def test_a_mode_that_worsens_but_stays_acceptable_is_not_a_regression(self): + """Mode 2 above goes from 5.7e-6 to 3.5e-5 — six times worse and + entirely irrelevant, since it is nowhere near the threshold.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([0.52, 5.7e-6]) + alt_r = np.array([1.1e-10, 3.5e-5]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 2) + assert improved[0] + assert not regressed.any() + + def test_rigid_body_noise_is_not_a_regression(self): + """Residuals that read ~1 in both candidates wobble either way. + A bare ``alt > sym`` test would call that a regression and block + every rescue that happens to sit beside a free-free mode.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([1.0, 0.8]) + alt_r = np.array([1.0001, 1.0e-9]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 2) + assert improved[1] + assert not regressed.any() + + def test_end_to_end_a_trading_candidate_is_declined(self, monkeypatch): + import pybmodes.fem.solver as solvermod + + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + real = solvermod._general_spectrum_for_retry + + def trading(gk_, gm_, n_modes): + vals, vecs, sound = real(gk_, gm_, n_modes) + # Corrupt the last returned mode so it is decisively worse. + vecs = vecs.copy() + vecs[:, -1] = np.roll(vecs[:, -1], 1) + return vals, vecs, sound + + monkeypatch.setattr( + solvermod, "_general_spectrum_for_retry", trading, + ) + with warnings.catch_warnings(): + warnings.simplefilter("error") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" + + class TestTheModeCountWarningStaysHonest: """A shortfall is only newsworthy when modes were actually discarded. From 38c2ecce04f585b141340310362fdc3ebb829816 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 18:38:02 +0900 Subject: [PATCH 12/28] fix: treat any threshold crossing as a regression (#140) Codex P1, and a real hole in the mirroring argument from the last commit. Defining regression as "above the threshold and decisively worse" lets a mode slide from 0.09 to 0.8 unremarked, because 8.9x falls short of the tenfold margin, and the candidate then replaces the whole spectrum with a newly failing mode. I had used the factor to keep rigid-body noise out. Codex points out that is unnecessary for this case: rigid modes read ~1, so they start above the threshold and were never on the acceptable side to cross from. A crossing test cannot false-positive on them. Regression is now either half. A crossing from acceptable to failing counts however small, and an already-failing mode counts when it is decisively worse. The first catches 0.09 to 0.8, the second catches 0.2 to 5.0 where there is no crossing to see, and neither fires on rigid-body wobble or on a mode six times worse but nowhere near the threshold. --- CHANGELOG.md | 9 +++++--- src/pybmodes/fem/solver.py | 29 ++++++++++++++++++-------- tests/fem/test_ill_conditioned_mass.py | 24 +++++++++++++++++++++ 3 files changed, 50 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 590bbc0..aca060d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,9 +47,12 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). spectrum, not the modes that prompted it, so a candidate that rescues one mode while pushing a previously acceptable one above the failure threshold is refused — it would hand back a new bad mode in place of an - old one. Regression mirrors improvement with the same threshold and - factor, so rigid-body modes reading ~1 in both candidates are not - mistaken for it. + old one. A mode counts as regressed when the candidate leaves it + failing and either it was acceptable before, any crossing counting + however small, or it was already failing and is now decisively worse. + Rigid-body modes read ~1 in both candidates, so they were never on the + acceptable side to cross from and their noise is not mistaken for a + regression. The comparison is made **per mode** rather than on the two maxima, so that rigid-body modes cannot distort it. Their backward error is a diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 790e990..4631d26 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -515,14 +515,19 @@ def _compare_candidate_modes( decisively better" is only half the test; the other half is that no mode got decisively worse. - Regression mirrors improvement exactly, using the same threshold and - the same factor rather than a second notion of acceptable: a mode has - regressed when it ends up above the failure threshold *and* is worse - there by the margin that would have counted as decisive in the other - direction. The symmetry matters for rigid-body modes, whose residual - is ~1 in both candidates and wobbles a little either way — that is - noise, not a regression, and a bare ``alt > sym`` test would read it - as one and block every rescue that happens to sit beside them. + A mode has regressed when the candidate leaves it above the failure + threshold **and** either it was acceptable before — any crossing + counts, however small — or it was already failing and is now + decisively worse by the same factor that defines an improvement. + + The two halves cover different things and both are needed. Without + the crossing test a mode sliding from 0.09 to 0.8 escapes, since that + is less than the tenfold margin; without the decisive test a mode + already at 0.2 could be driven to 5.0 unremarked. And the crossing + test cannot mistake rigid-body noise for a regression, which is what + a bare ``alt > sym`` comparison would do: those modes read ~1 in both + candidates and wobble either way, so they were never on the + acceptable side of the threshold to cross from. The comparison has to be **per mode**, not on the two maxima. A rigid-body mode's backward error is a ratio of two near-zero @@ -555,7 +560,13 @@ def _compare_candidate_modes( factor = _SOLVER_OPTIONS.residual_retry_improvement sym, alt = sym_r[:n], alt_r[:n] improved = (sym > threshold) & (alt < factor * sym) - regressed = (alt > threshold) & (sym < factor * alt) + # A mode regresses when the candidate leaves it failing, and either + # it was acceptable before — any crossing of the threshold counts, + # however small — or it was already failing and is now decisively + # worse. The crossing test cannot mistake rigid-body noise for a + # regression, because those modes read ~1 and so were never on the + # acceptable side to cross from. + regressed = (alt > threshold) & ((sym <= threshold) | (sym < factor * alt)) return improved, regressed diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 03664db..61a38f1 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -814,6 +814,30 @@ def test_a_trade_is_not_an_improvement(self): # The rule the caller applies. assert not (improved.any() and not regressed.any()) + def test_a_small_crossing_of_the_threshold_is_still_a_regression(self): + """A mode sliding from just under the threshold to well over it + worsens by less than the decisive factor, so the mirrored test + alone lets it through. Any crossing counts.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([0.52, 0.09]) + alt_r = np.array([1.0e-10, 0.8]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 2) + assert improved[0] + assert regressed[1] + assert not (improved.any() and not regressed.any()) + + def test_an_already_failing_mode_driven_much_worse_is_a_regression(self): + """The other half: no crossing, because it was failing already, + but a decisive worsening all the same.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([0.52, 0.2]) + alt_r = np.array([1.0e-10, 5.0]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 2) + assert improved[0] + assert regressed[1] + def test_a_mode_that_worsens_but_stays_acceptable_is_not_a_regression(self): """Mode 2 above goes from 5.7e-6 to 3.5e-5 — six times worse and entirely irrelevant, since it is nowhere near the threshold.""" From 8e4e49f406703d91b2ba0797e049786ebed380f0 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 18:42:27 +0900 Subject: [PATCH 13/28] fix: close the sub-threshold trade, and bound what the rule tolerates Found by sweeping the (sym, alt) plane rather than sweeping inputs, after three consecutive findings landed in the decision rule while my input audit found nothing there. Gating the decisive-worsening test on the failure threshold itself left a mode free to go from machine precision to just under the bar. 1e-16 to 0.1 is fifteen orders of magnitude and was not flagged, so a candidate could rescue one mode while quietly degrading every other one to the edge of tolerance. Same class as the crossing case, one step below it. The worsening test is now gated a tenth lower, reusing the two constants already in play rather than adding a third. The sweep is now two property tests rather than a one-off. One asserts that improvement and regression are mutually exclusive for every mode, so the caller rule cannot read a contradiction. The other states the honest limit: some degradation must be tolerated or rigid-body wobble would block every rescue, but nothing unflagged can leave a mode worse than a tenth of the failure threshold. That is a bound rather than a claim of perfection, and it is checked across eighteen decades. --- CHANGELOG.md | 12 +++--- src/pybmodes/fem/solver.py | 58 +++++++++++++++++--------- tests/fem/test_ill_conditioned_mass.py | 50 ++++++++++++++++++++++ 3 files changed, 94 insertions(+), 26 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index aca060d..0f052e8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,12 +47,12 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). spectrum, not the modes that prompted it, so a candidate that rescues one mode while pushing a previously acceptable one above the failure threshold is refused — it would hand back a new bad mode in place of an - old one. A mode counts as regressed when the candidate leaves it - failing and either it was acceptable before, any crossing counting - however small, or it was already failing and is now decisively worse. - Rigid-body modes read ~1 in both candidates, so they were never on the - acceptable side to cross from and their noise is not mistaken for a - regression. + old one. A mode counts as regressed if it crossed the failure threshold + at any size, or if it worsened decisively while ending up somewhere + that could matter. Some degradation is always tolerated, or rigid-body + wobble would block every rescue, but the tolerated region is bounded: + nothing unflagged can leave a mode worse than a tenth of the failure + threshold, an order of magnitude inside tolerance. The comparison is made **per mode** rather than on the two maxima, so that rigid-body modes cannot distort it. Their backward error is a diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 4631d26..fa4b9c3 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -515,19 +515,30 @@ def _compare_candidate_modes( decisively better" is only half the test; the other half is that no mode got decisively worse. - A mode has regressed when the candidate leaves it above the failure - threshold **and** either it was acceptable before — any crossing - counts, however small — or it was already failing and is now - decisively worse by the same factor that defines an improvement. - - The two halves cover different things and both are needed. Without - the crossing test a mode sliding from 0.09 to 0.8 escapes, since that - is less than the tenfold margin; without the decisive test a mode - already at 0.2 could be driven to 5.0 unremarked. And the crossing - test cannot mistake rigid-body noise for a regression, which is what - a bare ``alt > sym`` comparison would do: those modes read ~1 in both - candidates and wobble either way, so they were never on the - acceptable side of the threshold to cross from. + A mode has regressed in either of two independent ways. It **crossed** + the failure threshold, having been acceptable and no longer being so, + at any size — landing on the wrong side of the bar is what the bar is + for. Or it **worsened decisively**, by the same factor that defines + an improvement, while ending up somewhere that could matter. + + Both are needed, and the bounds on each were arrived at by finding + the cases the other misses. + + - Without the crossing test, a mode sliding from 0.09 to 0.8 escapes: + 8.9x falls short of the tenfold margin. + - Without the worsening test, an already-failing mode can be driven + from 0.2 to 5.0 with no crossing to observe. + - Gating the worsening test on the threshold itself, rather than a + tenth of it, lets an exact mode be driven to just under the bar — + 1e-16 to 0.1 is fifteen orders and passed unflagged. + + What must *not* be flagged bounds it from the other side. Rigid-body + modes read ~1 in both candidates and wobble either way: they never + cross, because they were never acceptable, and they never worsen + decisively, because the wobble is small. A mode going from 5.7e-6 to + 3.5e-5 is six times worse and three orders below anything that + matters. A bare ``alt > sym`` comparison would flag both and block + nearly every legitimate rescue. The comparison has to be **per mode**, not on the two maxima. A rigid-body mode's backward error is a ratio of two near-zero @@ -560,13 +571,20 @@ def _compare_candidate_modes( factor = _SOLVER_OPTIONS.residual_retry_improvement sym, alt = sym_r[:n], alt_r[:n] improved = (sym > threshold) & (alt < factor * sym) - # A mode regresses when the candidate leaves it failing, and either - # it was acceptable before — any crossing of the threshold counts, - # however small — or it was already failing and is now decisively - # worse. The crossing test cannot mistake rigid-body noise for a - # regression, because those modes read ~1 and so were never on the - # acceptable side to cross from. - regressed = (alt > threshold) & ((sym <= threshold) | (sym < factor * alt)) + # Two independent ways to regress. + # + # A crossing: the mode was acceptable and is not any more. Size does + # not matter here — landing on the wrong side of the bar is the whole + # point of having one. + crossed = (alt > threshold) & (sym <= threshold) + # A decisive worsening that stays on the acceptable side. Gated at a + # tenth of the threshold so an exact mode driven to just under the + # bar still counts, while genuinely negligible churn does not: 1e-16 + # to 0.1 is fifteen orders and matters, 5.7e-6 to 3.5e-5 is six times + # and does not. Reuses the two constants already in play rather than + # introducing a third. + worsened = (alt > factor * threshold) & (sym < factor * alt) + regressed = crossed | worsened return improved, regressed diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 61a38f1..439bf7b 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -827,6 +827,56 @@ def test_a_small_crossing_of_the_threshold_is_still_a_regression(self): assert regressed[1] assert not (improved.any() and not regressed.any()) + def test_an_exact_mode_driven_to_just_under_the_bar_is_a_regression(self): + """Gating the worsening test on the threshold itself would let a + mode go from machine precision to 0.1 — fifteen orders — while + staying nominally acceptable. It is gated a tenth lower.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([0.52, 1.0e-16]) + alt_r = np.array([1.0e-10, 0.099]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 2) + assert improved[0] + assert regressed[1] + + def test_what_escapes_is_bounded_an_order_inside_tolerance(self): + """The honest limit of the rule, swept rather than argued. + + Some degradation is always tolerated, or rigid-body wobble and + harmless churn would block every rescue. What matters is that the + tolerated region is bounded: nothing unflagged can leave a mode + worse than a tenth of the failure threshold, which is an order of + magnitude inside tolerance. + """ + from pybmodes.fem.solver import _compare_candidate_modes + from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt + + bound = opt.residual_retry_improvement * opt.residual_retry_threshold + grid = np.logspace(-16, 2, 37) + for s in grid: + for a in grid: + _imp, reg = _compare_candidate_modes( + np.array([s]), np.array([a]), 1, 1, + ) + if a > 10.0 * s and not reg[0]: + assert a <= bound, ( + f"sym={s:.2e} -> alt={a:.2e} escaped unflagged " + f"above the {bound:.2e} bound" + ) + + def test_improved_and_regressed_are_mutually_exclusive(self): + """A single mode cannot be both, or the caller's rule would be + reading a contradiction.""" + from pybmodes.fem.solver import _compare_candidate_modes + + grid = np.logspace(-16, 2, 37) + for s in grid: + for a in grid: + imp, reg = _compare_candidate_modes( + np.array([s]), np.array([a]), 1, 1, + ) + assert not (imp[0] and reg[0]), f"sym={s:.2e} alt={a:.2e}" + def test_an_already_failing_mode_driven_much_worse_is_a_regression(self): """The other half: no crossing, because it was failing already, but a decisive worsening all the same.""" From 55ff557f1d41ab8afea6ab7b447db88b2166fda4 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 18:49:21 +0900 Subject: [PATCH 14/28] fix: make the acceptance rule match the guarantee it advertises Two from Codex on #140, and the second is a criticism of a claim I made rather than of the code. Rigid-mode noise could justify a swap. I had been asserting for several rounds that a rigid residual reads about 1 in both candidates. That is false: both sides of the ratio are roundoff, so it is unbounded, and 12.39 from eigh against 0.794 from eig on a healthy pencil is a tenfold win on pure noise, enough to replace the whole spectrum. Improvement now mirrors regression and requires the mode to cross the threshold into acceptable, not merely to get better. 0.794 is still failing, so it justifies nothing. The advertised bound was not enforced. I claimed nothing unflagged could leave a mode worse than a tenth of the threshold, and the property test only swept tenfold worsenings, so it never tested the claim. A mode sliding from 0.02 to 0.099 is under fivefold, lands at the edge of tolerance and was invisible. Regression now flags any worsening of an acceptable mode that lands above that bound, at any ratio. Both fixes converge on one rule rather than accumulating clauses. The threshold defines trustworthy, and verdicts are only issued on the trustworthy side: a mode that was acceptable is judged in both directions, a mode already failing is judged in neither. That is what makes rigid-body modes tractable without identifying them, which two earlier attempts showed cannot be done here, and it is stated as the deliberate hole it is rather than left to be discovered. max_residual still reports the untrustworthy modes. The property test now sweeps every pair on the acceptable side instead of only the tenfold ones, so it certifies the guarantee rather than assuming it. --- CHANGELOG.md | 14 +-- src/pybmodes/fem/solver.py | 94 ++++++++++++-------- tests/fem/test_ill_conditioned_mass.py | 114 +++++++++++++++++++------ 3 files changed, 153 insertions(+), 69 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0f052e8..24755ca 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,12 +47,14 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). spectrum, not the modes that prompted it, so a candidate that rescues one mode while pushing a previously acceptable one above the failure threshold is refused — it would hand back a new bad mode in place of an - old one. A mode counts as regressed if it crossed the failure threshold - at any size, or if it worsened decisively while ending up somewhere - that could matter. Some degradation is always tolerated, or rigid-body - wobble would block every rescue, but the tolerated region is bounded: - nothing unflagged can leave a mode worse than a tenth of the failure - threshold, an order of magnitude inside tolerance. + old one. The guarantee is one-sided and precise: a mode that was + acceptable can end up above a tenth of the failure threshold only by + having improved, never as collateral of another mode's rescue. A mode + already failing carries no verdict either way — above the threshold + neither candidate is trustworthy, and a rigid-body mode, whose residual + divides one roundoff quantity by another and has been measured at 12.4 + against 0.79 on a healthy model, lives entirely in that region. + `max_residual` still reports it. The comparison is made **per mode** rather than on the two maxima, so that rigid-body modes cannot distort it. Their backward error is a diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index fa4b9c3..1ea74bb 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -515,31 +515,38 @@ def _compare_candidate_modes( decisively better" is only half the test; the other half is that no mode got decisively worse. - A mode has regressed in either of two independent ways. It **crossed** - the failure threshold, having been acceptable and no longer being so, - at any size — landing on the wrong side of the bar is what the bar is - for. Or it **worsened decisively**, by the same factor that defines - an improvement, while ending up somewhere that could matter. - - Both are needed, and the bounds on each were arrived at by finding - the cases the other misses. - - - Without the crossing test, a mode sliding from 0.09 to 0.8 escapes: - 8.9x falls short of the tenfold margin. - - Without the worsening test, an already-failing mode can be driven - from 0.2 to 5.0 with no crossing to observe. - - Gating the worsening test on the threshold itself, rather than a - tenth of it, lets an exact mode be driven to just under the bar — - 1e-16 to 0.1 is fifteen orders and passed unflagged. - - What must *not* be flagged bounds it from the other side. Rigid-body - modes read ~1 in both candidates and wobble either way: they never - cross, because they were never acceptable, and they never worsen - decisively, because the wobble is small. A mode going from 5.7e-6 to - 3.5e-5 is six times worse and three orders below anything that - matters. A bare ``alt > sym`` comparison would flag both and block + A mode has regressed when it was acceptable, comes back worse, and + lands somewhere that could matter — above a tenth of the threshold. + Any worsening counts at any ratio: the tenfold margin belongs to the + improvement side, and requiring it here left a mode free to slide + from 0.02 to 0.099 unflagged, which is the edge of tolerance. + + The landing bound is what keeps this usable rather than paralysing: + a mode going from 5.7e-6 to 3.5e-5 is six times worse and three + orders below anything that matters, and flagging it would block nearly every legitimate rescue. + A mode already failing in the symmetric solve gets **no verdict** in + either direction. That is the single rule that makes rigid-body modes + tractable without identifying them, which two earlier attempts showed + cannot be done reliably here. Their residual divides one roundoff + quantity by another, so it is *not* dependably near 1 — the ratio is + unbounded and has been measured at 12.4 from one solver against 0.79 + from the other on a perfectly healthy pencil. Read as an improvement + that is a tenfold win on pure noise; read as a regression it would + veto every rescue that happens to sit beside a free-free mode. + Declining to judge the untrustworthy side avoids both, and + ``max_residual`` still reports the mode to the caller. + + The cost is a real case declined: a breakdown the alternative + improves a hundredfold but leaves failing anyway is not acted on. + Neither result is trustworthy there, so keeping the original and + reporting the backward error is the honest outcome. + + Together the two verdicts give the guarantee the caller relies on: a + mode that was acceptable can only end up above a tenth of the + threshold by having *improved*, never as collateral. + The comparison has to be **per mode**, not on the two maxima. A rigid-body mode's backward error is a ratio of two near-zero quantities and reads ~1 in *both* candidates however exact each is, @@ -570,21 +577,34 @@ def _compare_candidate_modes( threshold = _SOLVER_OPTIONS.residual_retry_threshold factor = _SOLVER_OPTIONS.residual_retry_improvement sym, alt = sym_r[:n], alt_r[:n] - improved = (sym > threshold) & (alt < factor * sym) - # Two independent ways to regress. + # Improvement is the mirror of regression: the mode must cross the + # threshold the *other* way, from failing to acceptable, and do so + # decisively. Requiring the crossing — rather than a factor alone — + # is what keeps a rigid-body mode from justifying a swap. Its + # residual divides one roundoff quantity by another, so it is not + # merely "~1 in both candidates" as it first appears: the ratio is + # unbounded and can read 12.4 from one solver and 0.79 from the + # other on a perfectly healthy pencil. That is a tenfold "win" on + # pure noise, and it used to be enough to replace the whole spectrum. + # Demanding that the candidate actually *resolve* the mode ignores + # it, because 0.79 is still a failing residual. + improved = (sym > threshold) & (alt <= threshold) & (alt < factor * sym) + # A mode that was acceptable must not come back materially worse. + # Any worsening counts, at any ratio — the tenfold margin belongs to + # the improvement side, and requiring it here left a mode free to + # slide from 0.02 to 0.099 unflagged. What bounds this instead is + # where the mode *lands*: below a tenth of the threshold the change + # cannot matter, which is what keeps harmless churn (5.7e-6 to + # 3.5e-5) from blocking every rescue. # - # A crossing: the mode was acceptable and is not any more. Size does - # not matter here — landing on the wrong side of the bar is the whole - # point of having one. - crossed = (alt > threshold) & (sym <= threshold) - # A decisive worsening that stays on the acceptable side. Gated at a - # tenth of the threshold so an exact mode driven to just under the - # bar still counts, while genuinely negligible churn does not: 1e-16 - # to 0.1 is fifteen orders and matters, 5.7e-6 to 3.5e-5 is six times - # and does not. Reuses the two constants already in play rather than - # introducing a third. - worsened = (alt > factor * threshold) & (sym < factor * alt) - regressed = crossed | worsened + # Modes already failing in the symmetric solve get no verdict at all. + # Neither value is trustworthy there, and a rigid-body mode — whose + # residual divides roundoff by roundoff and has been seen to read + # 12.4 against 0.79 on a healthy pencil — lives entirely in that + # region. Judging it would be judging noise, and doing so in this + # direction would let that noise veto every legitimate rescue. + # ``max_residual`` still reports such a mode to the caller. + regressed = (sym <= threshold) & (alt > factor * threshold) & (alt > sym) return improved, regressed diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 439bf7b..0899114 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -839,30 +839,57 @@ def test_an_exact_mode_driven_to_just_under_the_bar_is_a_regression(self): assert improved[0] assert regressed[1] - def test_what_escapes_is_bounded_an_order_inside_tolerance(self): - """The honest limit of the rule, swept rather than argued. - - Some degradation is always tolerated, or rigid-body wobble and - harmless churn would block every rescue. What matters is that the - tolerated region is bounded: nothing unflagged can leave a mode - worse than a tenth of the failure threshold, which is an order of - magnitude inside tolerance. + def test_an_acceptable_mode_cannot_be_degraded_as_collateral(self): + """The guarantee, checked over **every** pair rather than only + the tenfold ones. + + An earlier version of this test swept only ``alt > 10 * sym``, + which is why it certified a bound the rule did not actually hold: + a mode sliding from 0.02 to 0.099 is under fivefold and lands at + the edge of tolerance, and nothing looked at it. + + The claim is narrow and about the acceptable side only. A mode + that was acceptable can end up above a tenth of the threshold + solely by having improved, never as collateral of someone else's + rescue. """ from pybmodes.fem.solver import _compare_candidate_modes from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt - bound = opt.residual_retry_improvement * opt.residual_retry_threshold + t = opt.residual_retry_threshold + bound = opt.residual_retry_improvement * t grid = np.logspace(-16, 2, 37) for s in grid: + if s > t: + continue # not the acceptable side for a in grid: - _imp, reg = _compare_candidate_modes( + imp, reg = _compare_candidate_modes( np.array([s]), np.array([a]), 1, 1, ) - if a > 10.0 * s and not reg[0]: - assert a <= bound, ( - f"sym={s:.2e} -> alt={a:.2e} escaped unflagged " - f"above the {bound:.2e} bound" - ) + if reg[0] or imp[0]: + continue + assert a <= bound or a <= s, ( + f"sym={s:.2e} -> alt={a:.2e} escaped unflagged above " + f"the {bound:.2e} bound without improving" + ) + + def test_a_mode_already_failing_gets_no_verdict(self): + """The deliberate hole, stated so it is not mistaken for one. + + Above the threshold neither candidate is trustworthy, and a + rigid-body mode lives entirely there. Judging that region either + way turns roundoff into a decision. + """ + from pybmodes.fem.solver import _compare_candidate_modes + from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt + + t = opt.residual_retry_threshold + for sym, alt in [(0.2, 5.0), (12.39, 0.794), (0.794, 12.39), (1.0, 1.0)]: + assert sym > t + _imp, reg = _compare_candidate_modes( + np.array([sym]), np.array([alt]), 1, 1, + ) + assert not reg[0], f"sym={sym} alt={alt} should carry no verdict" def test_improved_and_regressed_are_mutually_exclusive(self): """A single mode cannot be both, or the caller's rule would be @@ -877,17 +904,6 @@ def test_improved_and_regressed_are_mutually_exclusive(self): ) assert not (imp[0] and reg[0]), f"sym={s:.2e} alt={a:.2e}" - def test_an_already_failing_mode_driven_much_worse_is_a_regression(self): - """The other half: no crossing, because it was failing already, - but a decisive worsening all the same.""" - from pybmodes.fem.solver import _compare_candidate_modes - - sym_r = np.array([0.52, 0.2]) - alt_r = np.array([1.0e-10, 5.0]) - improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 2) - assert improved[0] - assert regressed[1] - def test_a_mode_that_worsens_but_stays_acceptable_is_not_a_regression(self): """Mode 2 above goes from 5.7e-6 to 3.5e-5 — six times worse and entirely irrelevant, since it is nowhere near the threshold.""" @@ -911,6 +927,52 @@ def test_rigid_body_noise_is_not_a_regression(self): assert improved[1] assert not regressed.any() + def test_rigid_body_noise_cannot_justify_a_swap(self): + """The assumption that rigid residuals sit near 1 in both + candidates is false: both sides of the ratio are roundoff, so it + is unbounded. Measured at 12.39 from ``eigh`` against 0.794 from + ``eig`` on a healthy pencil — a tenfold "win" on pure noise. + + Requiring the candidate to *resolve* the mode rather than merely + improve it ignores that, since 0.794 is still failing. + """ + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([12.39, 3.0e-15, 2.0e-15, 1.0e-15]) + alt_r = np.array([0.794, 1.0e-15, 2.0e-15, 3.0e-15]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 4, 4) + assert not improved.any() + assert not regressed.any() + + def test_a_healthy_free_free_pencil_is_left_alone(self): + """End to end on the shape from the report: rank-deficient K, a + well-conditioned M, exact elastic modes.""" + rng = np.random.default_rng(19) + a = rng.normal(size=(4, 4)) + gm = a @ a.T + 4.0 * np.eye(4) + b = rng.normal(size=(4, 3)) + gk = b @ b.T + gk, gm = 0.5 * (gk + gk.T), 0.5 * (gm + gm.T) + with warnings.catch_warnings(): + warnings.simplefilter("error") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert eigvals.size == 4 + + def test_an_improvement_that_leaves_the_mode_failing_is_declined(self): + """The cost of the rule, stated rather than hidden: a hundredfold + gain that still ends above the threshold is not acted on, because + neither result is trustworthy there.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([50.0]) + alt_r = np.array([0.5]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 1, 1) + assert not improved.any() + assert not regressed.any() + def test_end_to_end_a_trading_candidate_is_declined(self, monkeypatch): import pybmodes.fem.solver as solvermod From 5c39b22b06e547c36b547613010e2ce25faed9de Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 19:02:28 +0900 Subject: [PATCH 15/28] fix: judge a rescue by the size of the win, not by where it lands Codex P2 on #140, the third distinct failure on this axis and the one that finally showed why the previous two were the wrong shape. Rigid roundoff can land below the failure threshold. 0.848 to 0.0762 on a healthy model looks exactly like a mode being resolved, so requiring the candidate to cross into acceptable, which was the last fix, is no more sufficient than the factor test before it. No absolute threshold can work here, because the value of a ratio of two roundoff quantities is arbitrary. Its ratio is not. Measuring both populations: genuine rescues improve by 1e5 to 1e10, worst observed 3.98e1 to 3.17e-4; rigid noise improves by 11x to 16x. Four orders of separation, so the improvement factor moves from 10x to 1000x and sits between them with about 60x margin on the noise side and 100x on the rescue side. Those numbers are now a test, so the constants cannot drift away from the evidence they came from. An absolute bar stays as a second condition for what the ratio cannot see, a wildly broken 1e6 against a candidate at 100, but at 1e-3 rather than the 1e-6 I first reached for. 1e-6 was tighter than the genuine point-mass rescue at 3.17e-4 and silently disabled it, which the existing test caught. The regression floor is now its own option rather than derived from the improvement factor. Deriving it meant retuning one silently retuned the other, and they answer different questions. Swept 600 healthy free-free pencils of varying size and rank deficiency: no spurious swaps. --- src/pybmodes/fem/solver.py | 60 ++++++++++++----------- src/pybmodes/options.py | 42 ++++++++++++---- tests/fem/test_ill_conditioned_mass.py | 66 +++++++++++++++++++++----- 3 files changed, 121 insertions(+), 47 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 1ea74bb..d667790 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -526,22 +526,25 @@ def _compare_candidate_modes( orders below anything that matters, and flagging it would block nearly every legitimate rescue. - A mode already failing in the symmetric solve gets **no verdict** in - either direction. That is the single rule that makes rigid-body modes - tractable without identifying them, which two earlier attempts showed - cannot be done reliably here. Their residual divides one roundoff - quantity by another, so it is *not* dependably near 1 — the ratio is - unbounded and has been measured at 12.4 from one solver against 0.79 - from the other on a perfectly healthy pencil. Read as an improvement - that is a tenfold win on pure noise; read as a regression it would - veto every rescue that happens to sit beside a free-free mode. - Declining to judge the untrustworthy side avoids both, and - ``max_residual`` still reports the mode to the caller. + A mode already failing in the symmetric solve gets **no regression + verdict**. Above the threshold neither candidate is trustworthy, and + judging that region would let noise veto every rescue that happens to + sit beside a free-free mode. ``max_residual`` still reports it. + + Rigid-body modes are kept out of the *acceptance* side by the + resolution bar rather than by being identified, which three attempts + established cannot be done reliably here — not by eigenvalue scale, + not by strain, and not by which side of the failure threshold the + residual happens to fall on. Dividing one roundoff quantity by + another produces an arbitrary number: 12.4, 0.79 and 0.076 have all + been measured on healthy models, and the last is *below* the failure + threshold, so it looked exactly like a mode being resolved. The one + thing roundoff reliably does not do is land near machine precision. The cost is a real case declined: a breakdown the alternative - improves a hundredfold but leaves failing anyway is not acted on. - Neither result is trustworthy there, so keeping the original and - reporting the backward error is the honest outcome. + improves substantially without resolving is not acted on. Neither + result is trustworthy there, so keeping the original and reporting + the backward error is the honest outcome. Together the two verdicts give the guarantee the caller relies on: a mode that was acceptable can only end up above a tenth of the @@ -577,18 +580,20 @@ def _compare_candidate_modes( threshold = _SOLVER_OPTIONS.residual_retry_threshold factor = _SOLVER_OPTIONS.residual_retry_improvement sym, alt = sym_r[:n], alt_r[:n] - # Improvement is the mirror of regression: the mode must cross the - # threshold the *other* way, from failing to acceptable, and do so - # decisively. Requiring the crossing — rather than a factor alone — - # is what keeps a rigid-body mode from justifying a swap. Its - # residual divides one roundoff quantity by another, so it is not - # merely "~1 in both candidates" as it first appears: the ratio is - # unbounded and can read 12.4 from one solver and 0.79 from the - # other on a perfectly healthy pencil. That is a tenfold "win" on - # pure noise, and it used to be enough to replace the whole spectrum. - # Demanding that the candidate actually *resolve* the mode ignores - # it, because 0.79 is still a failing residual. - improved = (sym > threshold) & (alt <= threshold) & (alt < factor * sym) + # What separates a rescue from noise is the *size* of the win, not + # which side of a line the candidate lands on. A rigid-body mode's + # residual divides one near-zero quantity by another, so its value is + # arbitrary — 12.4, 0.79 and 0.076 have all been measured on healthy + # models, and the last is below the failure threshold, so no absolute + # threshold can exclude it. Its *ratio*, though, stays around 11x to + # 16x, while a genuine rescue improves by 1e5 to 1e10. Four orders + # separate the two populations. + # + # The absolute bar is kept as a second condition for the case the + # ratio cannot see: a wildly broken 1e6 against a candidate at 100 + # clears any ratio while both remain garbage. + resolved = _SOLVER_OPTIONS.residual_retry_resolved + improved = (sym > threshold) & (alt <= resolved) & (alt < factor * sym) # A mode that was acceptable must not come back materially worse. # Any worsening counts, at any ratio — the tenfold margin belongs to # the improvement side, and requiring it here left a mode free to @@ -604,7 +609,8 @@ def _compare_candidate_modes( # region. Judging it would be judging noise, and doing so in this # direction would let that noise veto every legitimate rescue. # ``max_residual`` still reports such a mode to the caller. - regressed = (sym <= threshold) & (alt > factor * threshold) & (alt > sym) + floor = _SOLVER_OPTIONS.residual_regression_floor + regressed = (sym <= threshold) & (alt > floor) & (alt > sym) return improved, regressed diff --git a/src/pybmodes/options.py b/src/pybmodes/options.py index d738c0a..09ddfff 100644 --- a/src/pybmodes/options.py +++ b/src/pybmodes/options.py @@ -85,6 +85,18 @@ class SolverOptions: ~4e10) reaches ~2e-2 and is *not* meant to trigger; the degraded regime starts around 0.7. The default splits the last two gaps with roughly 5x either side. + residual_retry_resolved : float, default 1e-3 + Backward error the candidate must reach on a mode before that + mode can justify a swap. Guards the case the ratio alone cannot: + a wildly broken symmetric solve at 1e6 against a candidate at + 100 clears any ratio while both remain garbage. + residual_regression_floor : float, default 1e-2 + Where a worsened mode has to land before the worsening counts. + Below this the change cannot matter, which is what stops harmless + churn — 5.7e-6 to 3.5e-5 — from blocking every rescue. Kept + separate from ``residual_retry_improvement`` on purpose: they + answer different questions, and deriving one from the other + coupled two unrelated decisions. residual_retry_max_ndof : int, default 2000 Largest reduced system the retry will attempt. The retry is a dense ``eig``, whose cost grows as ``ngd^3`` with a much larger @@ -96,21 +108,35 @@ class SolverOptions: already available. A guard against a silent wrong answer should not be able to turn one into a silent hang, so above this size it declines and leaves the backward error to the diagnostics. - residual_retry_improvement : float, default 0.1 - How much better the general path's backward error must be before - its result is taken. The second, and more important, guard: on - that same land deck the general path is only ~1.4x better while + residual_retry_improvement : float, default 1e-3 + How much better the general path's backward error must be on a + mode before that mode can justify taking its result. The + load-bearing guard, and the one that separates a rescue from + noise. + + Two populations were measured while building this. Genuine + rescues improve by 1e5 to 1e10: the worst observed went from + 3.98e1 to 3.17e-4. Rigid-body roundoff, which divides one + near-zero quantity by another and so produces an arbitrary + number, improves by 11x to 16x: 0.848 to 0.0762 on a healthy + model, which no absolute threshold excludes because 0.0762 sits + below the failure line. Four orders separate the two, and the + default sits in the middle with roughly 60x margin on the noise + side and 100x on the rescue side. + + A marginal win is refused for a second reason as well: on the + bundled NREL 5MW land deck the general path is 1.4x better while *breaking* a physically real degenerate fore-aft / side-side pair the symmetric solver resolves exactly, which the downstream FA / - SS classifier depends on. A genuine breakdown is not marginal — - it improves by nine orders — so requiring a decisive win keeps - validated results untouched and still catches the real failure. + SS classifier depends on. """ sparse_ndof_threshold: int = 500 symmetry_rtol: float = 1.0e-12 residual_retry_threshold: float = 0.1 - residual_retry_improvement: float = 0.1 + residual_retry_improvement: float = 1.0e-3 + residual_retry_resolved: float = 1.0e-3 + residual_regression_floor: float = 1.0e-2 residual_retry_max_ndof: int = 2000 diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 0899114..d65e91a 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -857,7 +857,7 @@ def test_an_acceptable_mode_cannot_be_degraded_as_collateral(self): from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt t = opt.residual_retry_threshold - bound = opt.residual_retry_improvement * t + bound = opt.residual_regression_floor grid = np.logspace(-16, 2, 37) for s in grid: if s > t: @@ -927,23 +927,65 @@ def test_rigid_body_noise_is_not_a_regression(self): assert improved[1] assert not regressed.any() - def test_rigid_body_noise_cannot_justify_a_swap(self): - """The assumption that rigid residuals sit near 1 in both - candidates is false: both sides of the ratio are roundoff, so it - is unbounded. Measured at 12.39 from ``eigh`` against 0.794 from - ``eig`` on a healthy pencil — a tenfold "win" on pure noise. - - Requiring the candidate to *resolve* the mode rather than merely - improve it ignores that, since 0.794 is still failing. + @pytest.mark.parametrize("alt_rigid", [12.39, 0.794, 0.0762, 0.05]) + def test_rigid_body_noise_cannot_justify_a_swap(self, alt_rigid): + """Rigid residuals divide roundoff by roundoff, so the value is + arbitrary — 12.39, 0.794 and 0.0762 have all been measured on + healthy models, and the last two sit *below* the failure + threshold. No threshold can exclude them; a resolution bar near + machine precision can, because roundoff does not land there. """ from pybmodes.fem.solver import _compare_candidate_modes - sym_r = np.array([12.39, 3.0e-15, 2.0e-15, 1.0e-15]) - alt_r = np.array([0.794, 1.0e-15, 2.0e-15, 3.0e-15]) - improved, regressed = _compare_candidate_modes(sym_r, alt_r, 4, 4) + sym_r = np.array([0.848, 3.0e-15, 2.0e-15]) + alt_r = np.array([alt_rigid, 1.0e-15, 2.0e-15]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 3, 3) assert not improved.any() assert not regressed.any() + def test_the_two_populations_are_separated_by_the_ratio(self): + """The measurement the rule is calibrated on, kept as a test so + the constants cannot drift away from their evidence. + + The mode that justifies a swap is judged by *how much* the + candidate improves it, because a rigid residual's value is + arbitrary while its ratio is not. Genuine rescues improve by 1e5 + to 1e10; rigid roundoff by 11x to 16x. + """ + from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt + + rescues = [(1.56e0, 4.10e-10), (3.01e0, 1.96e-09), + (7.55e1, 4.93e-09), (3.98e1, 3.17e-04)] + noise = [(12.39, 0.794), (0.848, 0.0762)] + + worst_rescue = max(a / s for s, a in rescues) + best_noise = min(a / s for s, a in noise) + assert worst_rescue < opt.residual_retry_improvement < best_noise + # And the absolute bar admits every rescue. + assert max(a for _s, a in rescues) <= opt.residual_retry_resolved + + @pytest.mark.parametrize("seed", [14946, 7, 101, 2024, 55555]) + def test_healthy_free_free_pencils_never_trigger_a_swap(self, seed): + """The empirical half of the argument, over several draws rather + than one. A rank-deficient K with a well-conditioned M is a + healthy free-free model; whatever its null-mode roundoff happens + to read, it must never replace the spectrum.""" + rng = np.random.default_rng(seed) + n = 3 + a = rng.normal(size=(n, n)) + gm = a @ a.T + n * np.eye(n) + b = rng.normal(size=(n, n - 1)) + gk = b @ b.T + gk, gm = 0.5 * (gk + gk.T), 0.5 * (gm + gm.T) + assert np.linalg.cond(gm) < 100.0 # genuinely healthy + with warnings.catch_warnings(): + warnings.simplefilter("error") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=n, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert eigvals.size == n + def test_a_healthy_free_free_pencil_is_left_alone(self): """End to end on the shape from the report: rank-deficient K, a well-conditioned M, exact elastic modes.""" From 14f0b299a3a46114a32a38630babdf5dbf9bfd15 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 19:10:50 +0900 Subject: [PATCH 16/28] docs: bring the prose back in step with the rule it describes A coherence pass over a change whose design shifted several times. The comments still said the improvement margin was tenfold after it became a thousand, still described the regression floor as derived from that margin after it became its own option, and the CHANGELOG still said the retry is taken on an order-of-magnitude gain. The test module docstring is now organised around the three readings of the backward error that turned out to be wrong, since that is what almost every test in the file exists to pin: a large error is not evidence of a breakdown, a small one is not evidence of a rescue, and a better mode does not make a better spectrum. Whoever touches this next should meet the failure history before the code. --- CHANGELOG.md | 3 ++- src/pybmodes/fem/solver.py | 20 +++++++------- tests/fem/test_ill_conditioned_mass.py | 36 +++++++++++++++++--------- 3 files changed, 36 insertions(+), 23 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 24755ca..9b34607 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,7 +22,8 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). `solve_modes` now checks the backward error `||K x - λ M x|| / ||K x||` of a **dense** symmetric solve and, when it is large, redoes it through the general dense path, which factorises neither matrix. The retried - result is taken only when it is better by an order of magnitude, and a + result is taken only when it resolves a mode the symmetric solve had + failed, by a margin that separates a rescue from roundoff, and a `RuntimeWarning` names the swap. `SolverDiagnostics` gains `residual_fallback` recording it. diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index d667790..0fd3b07 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -516,10 +516,10 @@ def _compare_candidate_modes( mode got decisively worse. A mode has regressed when it was acceptable, comes back worse, and - lands somewhere that could matter — above a tenth of the threshold. - Any worsening counts at any ratio: the tenfold margin belongs to the - improvement side, and requiring it here left a mode free to slide - from 0.02 to 0.099 unflagged, which is the edge of tolerance. + lands above the regression floor. Any worsening counts at any ratio: + the improvement side's margin answers a different question — rescue + or noise — and borrowing it here left a mode free to slide from 0.02 + to 0.099 unflagged, which is the edge of tolerance. The landing bound is what keeps this usable rather than paralysing: a mode going from 5.7e-6 to 3.5e-5 is six times worse and three @@ -595,12 +595,12 @@ def _compare_candidate_modes( resolved = _SOLVER_OPTIONS.residual_retry_resolved improved = (sym > threshold) & (alt <= resolved) & (alt < factor * sym) # A mode that was acceptable must not come back materially worse. - # Any worsening counts, at any ratio — the tenfold margin belongs to - # the improvement side, and requiring it here left a mode free to - # slide from 0.02 to 0.099 unflagged. What bounds this instead is - # where the mode *lands*: below a tenth of the threshold the change - # cannot matter, which is what keeps harmless churn (5.7e-6 to - # 3.5e-5) from blocking every rescue. + # Any worsening counts, at any ratio: the improvement side's margin + # answers "is this a rescue or noise", which is a different question, + # and borrowing it here left a mode free to slide from 0.02 to 0.099 + # unflagged. What bounds this instead is where the mode *lands* — + # below the regression floor the change cannot matter, which is what + # keeps harmless churn (5.7e-6 to 3.5e-5) from blocking every rescue. # # Modes already failing in the symmetric solve get no verdict at all. # Neither value is trustworthy there, and a rigid-body mode — whose diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index d65e91a..0334f2e 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -9,16 +9,27 @@ symmetric path reported 0.103 Hz against a true 0.0436 Hz, a factor of 2.4, with no error and no warning. -The guard is the backward error ``||K x - lambda M x|| / ||K x||``. It -has two conditions, and the second matters more than the first: the -error must exceed the retry threshold, *and* the general path must beat -it decisively. Being above the threshold alone is not evidence of a -breakdown — the bundled NREL 5MW land deck sits at ~2e-2 because its -adapter leaves ``M`` at cond ~4e10, and there the general path is only -1.4x better while splitting a degenerate fore-aft / side-side pair the -symmetric solver resolves exactly. Swapping for that would churn a -validated frequency by 0.84 % and break the FA / SS classifier -downstream. A real breakdown improves by nine orders, not by a factor. +The guard is the backward error ``||K x - lambda M x|| / ||K x||``, and +almost every test here exists because some reading of it turned out to +be wrong. Three things it does *not* establish, each learned the hard +way and each pinned below: + +- **A large error is not evidence of a breakdown.** The bundled NREL 5MW + land deck sits at ~2e-2 because its adapter leaves ``M`` at cond ~4e10; + swapping there churns a validated frequency by 0.84 % and splits a + degenerate fore-aft / side-side pair the symmetric solver resolves + exactly, which the FA / SS classifier depends on. +- **A small error is not evidence of a rescue.** A rigid-body mode's + residual divides one roundoff quantity by another, so its value is + arbitrary: 12.4, 0.79 and 0.076 have all been measured on healthy + models, and the last sits *below* the failure threshold. No absolute + bar can separate that from a solved mode. +- **A better mode does not make a better spectrum.** Accepting replaces + every mode, so a candidate that rescues one while pushing another past + the threshold is a trade, not an improvement. + +What does separate the populations is the *size* of the win. Genuine +rescues improve by 1e5 to 1e10; rigid roundoff by 11x to 16x. Analytical reference: a cantilever whose beam mass is negligible next to a tip lump behaves as a spring-mass oscillator on the static tip @@ -841,7 +852,8 @@ def test_an_exact_mode_driven_to_just_under_the_bar_is_a_regression(self): def test_an_acceptable_mode_cannot_be_degraded_as_collateral(self): """The guarantee, checked over **every** pair rather than only - the tenfold ones. + the tenfold ones (an earlier version of this test swept only + those, which is how it certified a bound the rule did not hold). An earlier version of this test swept only ``alt > 10 * sym``, which is why it certified a bound the rule did not actually hold: @@ -849,7 +861,7 @@ def test_an_acceptable_mode_cannot_be_degraded_as_collateral(self): the edge of tolerance, and nothing looked at it. The claim is narrow and about the acceptable side only. A mode - that was acceptable can end up above a tenth of the threshold + that was acceptable can end up above the regression floor solely by having improved, never as collateral of someone else's rescue. """ From b22b93897119de1d5b284d970f076d092b9e2ba1 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 19:15:48 +0900 Subject: [PATCH 17/28] fix: attribute the retry warning to the cause that was measured Codex P2 on #140. The warning blamed a nearly singular mass matrix unconditionally, but the near-singular mass is the motivating case, not the only one. A stiffness spectrum spanning 1e-16 to 1 with M = I trips the same guard at cond(M) = 1, and the message then sends the reader to check a mass distribution that is perfectly fine, or fails a healthy solve in a warnings-as-errors environment with a wrong explanation. The cause is now measured before it is asserted. The condition number is computed only on this branch, which is rare, and skipped when the system is large enough for the estimate to cost something, where the text names both possibilities instead of guessing between them. The CHANGELOG also now states the actual acceptance rule, a thousandfold improvement and a resolved backward error, rather than the order-of-magnitude figure it inherited from an earlier revision. Codex flagged that too; the coherence pass had softened the wording without supplying the number. --- CHANGELOG.md | 13 ++++-- src/pybmodes/fem/solver.py | 60 +++++++++++++++++++++++--- tests/fem/test_ill_conditioned_mass.py | 38 ++++++++++++++++ 3 files changed, 101 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9b34607..567ef29 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,9 +22,16 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). `solve_modes` now checks the backward error `||K x - λ M x|| / ||K x||` of a **dense** symmetric solve and, when it is large, redoes it through the general dense path, which factorises neither matrix. The retried - result is taken only when it resolves a mode the symmetric solve had - failed, by a margin that separates a rescue from roundoff, and a - `RuntimeWarning` names the swap. `SolverDiagnostics` gains + result is taken only when it **resolves** a mode the symmetric solve + had failed — improving it by more than 1000× and reaching a backward + error of 1e-3 or better. Both conditions are needed because a + rigid-body mode's residual divides one roundoff quantity by another, + so its value is arbitrary (0.076, 0.79 and 12.4 have all been measured + on healthy models) while its improvement ratio stays near 10×, an + order of magnitude short of any real rescue. A `RuntimeWarning` names + the swap, and attributes it to the mass matrix only when the mass + conditioning supports that — a wide stiffness range trips the same + guard with a perfectly conditioned mass. `SolverDiagnostics` gains `residual_fallback` recording it. The **sparse** path is deliberately not retried and never sets diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 0fd3b07..deff0ef 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -409,13 +409,9 @@ def solve_modes( f"{int(improved.sum())} mode(s) that do not satisfy " f"K x = lambda M x — worst at index {idx}, backward " f"error {sym_r[idx]:.2e} against {alt_r[idx]:.2e} from " - f"the general dense path. Its Cholesky reduction of the " - f"mass matrix loses accuracy when that matrix is nearly " - f"singular, which a very light beam carrying a very " - f"heavy lump produces. The returned modes come from the " - f"general solve, which factorises neither matrix. Worth " - f"checking the mass distribution is the one you " - f"intended.", + f"the general dense path. The returned modes come from " + f"the general solve, which factorises neither matrix. " + + _retry_cause(gm_s), RuntimeWarning, stacklevel=2, ) @@ -500,6 +496,56 @@ def _build_diagnostics( # one is a rigid-body mode: a free-free floating platform has up to six, # and an unrestrained DOF (a symmetric column's yaw) gives an exactly # zero one. +# Above this the mass matrix is ill-conditioned enough for the Cholesky +# reduction to be the credible culprit; below it, something else in the +# pencil is. +_MASS_COND_ATTRIBUTION = 1.0e8 + + +def _retry_cause(gm: np.ndarray) -> str: + """The explanatory half of the retry warning, attributed honestly. + + The near-singular mass matrix is the *motivating* case, not the only + one: the symmetric reduction degrades on an ill-conditioned pencil + generally, and a stiffness spectrum spanning 1e-16 to 1 with ``M = I`` + triggers this guard while ``cond(M) = 1``. Naming the mass matrix + there would send the reader to check a mass distribution that is + perfectly fine. + + So the cause is measured before it is asserted. The condition number + is only computed on this branch, which is rare, and is skipped for a + system large enough for the O(n^3) estimate to matter — where the + text falls back to naming both possibilities. + """ + if gm.shape[0] > _COND_DENSE_MAX: + return ( + "This happens when the pencil is ill-conditioned — most often " + "a nearly singular mass matrix, from a very light beam " + "carrying a very heavy lump, but a very wide stiffness range " + "does it too. Worth checking the section properties for an " + "extreme mass or stiffness ratio." + ) + try: + cond = float(np.linalg.cond(gm)) + except np.linalg.LinAlgError: + cond = float("inf") + if cond > _MASS_COND_ATTRIBUTION: + return ( + f"The symmetric reduction goes through a Cholesky factor of " + f"the mass matrix, which is nearly singular here " + f"(cond = {cond:.1e}) — a very light beam carrying a very " + f"heavy lump does this. Worth checking the mass distribution " + f"is the one you intended." + ) + return ( + f"The mass matrix is well conditioned (cond = {cond:.1e}), so the " + f"reduction was defeated by the pencil rather than by the mass: a " + f"stiffness range wide enough to put a soft mode at the level of " + f"roundoff will do it. Worth checking the section properties for " + f"an extreme stiffness ratio." + ) + + def _compare_candidate_modes( sym_r: np.ndarray, alt_r: np.ndarray, diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 0334f2e..fab6634 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -663,6 +663,44 @@ def test_a_drop_inside_the_window_makes_the_ordering_unsound(self): assert vals.size == 2 +class TestTheWarningAttributesTheCauseHonestly: + """The near-singular mass matrix motivated this guard but is not the + only thing that trips it, and naming it unconditionally sends the + reader to check something that may be perfectly fine.""" + + def test_a_singular_mass_names_the_mass_matrix(self): + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with pytest.warns(RuntimeWarning, match="nearly singular here"): + solve_modes(gk, gm, n_modes=4) + + def test_a_well_conditioned_mass_does_not_blame_it(self): + """``M = I`` with a stiffness spectrum spanning 1e-16 to 1 trips + the guard at ``cond(M) = 1``.""" + rng = np.random.default_rng(474) + q, _r = np.linalg.qr(rng.normal(size=(3, 3))) + gk = q @ np.diag([1.0e-16, 1.0e-8, 1.0]) @ q.T + gk = 0.5 * (gk + gk.T) + gm = np.eye(3) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + _v, _x, diag = solve_modes(gk, gm, n_modes=3, + return_diagnostics=True) + text = " ".join(str(w.message) for w in caught) + if diag.residual_fallback: + assert "well conditioned" in text + assert "nearly singular" not in text + assert "stiffness ratio" in text + + def test_the_message_states_what_was_measured(self): + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + solve_modes(gk, gm, n_modes=4) + text = " ".join(str(w.message) for w in caught) + assert "do not satisfy" in text + assert "cond =" in text + + class TestTheRetryIsScopedToTheDensePath: """Only the dense symmetric path is retried, and that is about which matrix each routine factorises. From 282ed9be299ceab722efa19aa4ebd8324fa2466e Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 19:26:20 +0900 Subject: [PATCH 18/28] docs: rewrite the guard prose against the code instead of patching it Three more from Codex, all documentation drift, which makes six on this PR. Patching the reported phrase each time has plainly not worked, so this stops doing that. Reported: the test module still described eigsh as reducing through a Cholesky factor of the mass matrix, when it factorises K and that distinction is the whole reason sparse results are exempt from the retry. An orphaned comment still described a rigid-mode eigenvalue cutoff that three attempts abandoned and no code implements. The CHANGELOG enumerated three of the five new SolverOptions fields. Found by scanning rather than reported: the module docstring still cited nine orders where the measurement is 1e5 to 1e10, still described the regression floor as a tenth of the threshold after it became its own option, still referenced a preserve_full_spectrum flag replaced by a dedicated helper, and still said rigid classification was tried twice when it was three times. The guard section of the module docstring is rewritten rather than patched, organised around the rule and the reason each clause exists, since every one of them is there because a simpler version was wrong. A first attempt at auditing this mechanically only caught what I had already thought to look for, which is why four of the seven above needed a manual read. The lasting fix is the structure: one place states each rule, and it is written next to the constant it describes. --- CHANGELOG.md | 19 +-- src/pybmodes/fem/solver.py | 159 +++++++++++++------------ tests/fem/test_ill_conditioned_mass.py | 19 +-- 3 files changed, 108 insertions(+), 89 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 567ef29..4514d32 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -27,8 +27,8 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). error of 1e-3 or better. Both conditions are needed because a rigid-body mode's residual divides one roundoff quantity by another, so its value is arbitrary (0.076, 0.79 and 12.4 have all been measured - on healthy models) while its improvement ratio stays near 10×, an - order of magnitude short of any real rescue. A `RuntimeWarning` names + on healthy models) while its improvement ratio stays near 10x, four + orders short of the 1e5 to 1e10 a real rescue achieves. A `RuntimeWarning` names the swap, and attributes it to the mass matrix only when the mass conditioning supports that — a wide stiffness range trips the same guard with a perfectly conditioned mass. `SolverDiagnostics` gains @@ -56,8 +56,8 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). one mode while pushing a previously acceptable one above the failure threshold is refused — it would hand back a new bad mode in place of an old one. The guarantee is one-sided and precise: a mode that was - acceptable can end up above a tenth of the failure threshold only by - having improved, never as collateral of another mode's rescue. A mode + acceptable can end up above the regression floor only by having + improved, never as collateral of another mode's rescue. A mode already failing carries no verdict either way — above the threshold neither candidate is trustworthy, and a rigid-body mode, whose residual divides one roundoff quantity by another and has been measured at 12.4 @@ -92,9 +92,14 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). symmetric result and its diagnostics are kept rather than the whole solve failing. -- `SolverOptions` gains `residual_retry_threshold`, - `residual_retry_improvement` and `residual_retry_max_ndof` for the - conditions above. +- `SolverOptions` gains five fields for the conditions above, each + governing one of them: `residual_retry_threshold` (what counts as a + failing mode), `residual_retry_improvement` (how much better the + candidate must be for the win to be a rescue rather than roundoff), + `residual_retry_resolved` (the backward error it must actually reach), + `residual_regression_floor` (where a worsened mode has to land before + the worsening counts) and `residual_retry_max_ndof` (the size above + which the retry is not attempted at all). ## [1.18.0] — 2026-08-12 diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index deff0ef..d22adbb 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -36,72 +36,84 @@ the retry path when a symmetric solve comes back with a large backward error — see below. -The **dense** symmetric path reduces ``K x = λ M x`` through a Cholesky -factor of the mass matrix, and that reduction degrades once ``M`` is -nearly singular, which a very light beam carrying a very heavy lump -produces. The failure mode is silent: LAPACK returns confidently wrong -low modes rather than raising. (The sparse path factorises ``K`` -instead, so it is unaffected and is not retried.) -:func:`solve_modes` therefore checks the backward error of a dense -symmetric solve and, when it exceeds -:attr:`~pybmodes.options.SolverOptions.residual_retry_threshold`, tries -the general path as well — taking its result only if it is better by -:attr:`~pybmodes.options.SolverOptions.residual_retry_improvement`, and -warning when it does. - -That second condition is the load-bearing one, and it is what lets the -check be simple. A real deck can sit above the threshold without being -broken (the bundled NREL 5MW land tower reaches ~2e-2, its adapter -leaving ``M`` at cond ~4e10), and there the general path is only -marginally better while *splitting* the degenerate fore-aft / side-side -pair the symmetric solver resolves exactly — which the FA / SS -classifier downstream depends on. A true breakdown is not marginal: it -improves by nine orders of magnitude. - -The comparison is made **per mode** rather than on the two maxima, which -is what keeps rigid-body modes from distorting it. Their backward error -is a ratio of two near-zero quantities and reads ~1 in both candidates -however exact each is; on a maximum that puts a floor under the -alternative and hides a genuinely corrupted elastic mode alongside them, -while per mode they simply register as ~1 against ~1, i.e. no -improvement. Identifying such modes and excluding them was tried twice -and abandoned — neither their eigenvalue nor their strain separates them -reliably from a genuinely soft mode. - -The retry additionally runs with ``preserve_full_spectrum=True``, which -drops the sign filter the general path normally applies. ``eigh`` -filters nothing, so keeping it would return a *different set* of modes — -the same length, since the gap is backfilled from higher up — and the -per-index comparison would be reading two different spectra against each -other, able to accept a result that had quietly dropped a mode and -shifted every one above it. Both omissions are reachable: a free-free -model's zero-frequency modes, and the negative eigenvalues an indefinite -``K`` produces once ``run(gravity=...)`` loads a column past its -buckling weight. - -Both the measurement and the retry use the **symmetrised** matrices, the -ones the symmetric paths actually solve. The skew they discard is only -guaranteed small relative to ``max|K|``, which in a model with a wide -dynamic range can still be large relative to a soft mode's own -eigenvalue; judging an exact symmetric solve against the unsymmetrised -matrices would then read as a failure, and ``eig`` on those same -matrices would "win decisively" purely by answering a different +The residual retry +------------------ + +``scipy.linalg.eigh`` reduces ``K x = λ M x`` through a Cholesky factor +of the **mass** matrix, and that reduction degrades once ``M`` is nearly +singular — a very light beam carrying a very heavy lump. The failure is +silent: LAPACK returns confidently wrong low modes rather than raising. +On the case in ``tests/fem/test_ill_conditioned_mass.py`` it reported +0.103 Hz against a true 0.0436 Hz. + +:func:`solve_modes` therefore measures the backward error of a dense +symmetric solve and, above +:attr:`~pybmodes.options.SolverOptions.residual_retry_threshold`, solves +again through the general path and compares. Every other rule below +exists because some simpler version of that comparison was wrong. + +**Only the dense path.** ``eigsh(sigma=0, mode='normal')`` factorises +``K``, not ``M``, so the sparse path does not have this failure and is +never retried. That also avoids comparing two different mode sets: its +``which="LM"`` window selects the modes nearest zero in magnitude while +the retry selects the algebraically smallest. + +**Per mode, not on the maxima.** A rigid-body mode's residual divides one +roundoff quantity by another. Taking maxima lets that noise floor the +candidate's worst value and hide a genuinely corrupted elastic mode +beside it. + +**Judged by the size of the win.** No absolute bar separates a rescue +from rigid noise, because the noise value is arbitrary — 0.076, 0.79 and +12.4 have all been measured on healthy models, and the first is *below* +the failure threshold. Identifying such modes was tried three times and +abandoned: not by eigenvalue scale, which a rigid-only subset makes its +own reference; not by strain, which a genuinely soft mode also has +little of; and not by which side of the threshold the value falls on. +What does separate them is the ratio. Rescues improve by 1e5 to 1e10, +roundoff by 11x to 16x, so acceptance needs +:attr:`~pybmodes.options.SolverOptions.residual_retry_improvement` *and* +a candidate that reaches +:attr:`~pybmodes.options.SolverOptions.residual_retry_resolved`. + +**Non-regressive.** Accepting replaces the whole spectrum, so a candidate +that rescues one mode while pushing another past the threshold is a +trade, not an improvement. + +**The same matrices throughout.** Both symmetric paths symmetrise +internally, so the measurement and the retry use the symmetrised pair. +The tolerated skew is small only relative to ``max|K|``, which in a +wide-dynamic-range model can still swamp a soft mode's own eigenvalue; +judging an exact solve against the raw matrices reads as a failure, and +``eig`` on those same matrices then "wins" by answering a different question. -**What this guard does not promise.** It rescues the case it was built -for — a near-singular mass matrix, with no rigid-body modes — reliably -and identically on every platform. It is *safe* everywhere else but not -always *effective*: when rigid-body modes and a near-singular mass -matrix coincide, QZ may represent the theoretically real zero modes as -small complex-conjugate pairs that cannot be coerced back to real, and -where those land in the spectrum differs between LAPACK builds. When -they land inside the requested window the alternative's ordering cannot -be verified, so the retry declines and the symmetric result stands. -Declining is the deliberate choice: a guard added to stop a silent wrong -answer must never introduce one, and backfilling a dropped zero mode -with an elastic mode would do exactly that. The result in that situation -is no worse than without the guard, and ``max_residual`` still reports -the problem. +**The same spectrum throughout.** The retry keeps every real eigenvalue, +zeros and negatives included, and verifies that nothing was discarded +from inside the returned window. ``eigh`` filters nothing, so any filter +here would return a different set of the same length, backfilled from +higher up, and equal indices would stop meaning equal modes. Both +omissions are reachable: a free-free model's zero modes, and the negative +eigenvalues an indefinite ``K`` produces once ``run(gravity=...)`` loads +a column past its buckling weight. + +**It can always decline.** If the alternative raises on the same +defective pencil, or the system is larger than +:attr:`~pybmodes.options.SolverOptions.residual_retry_max_ndof`, or its +ordering cannot be verified, the symmetric result stands. + +What this does not promise +-------------------------- + +The rescue is reliable and platform-independent for the case it was built +for: a near-singular mass matrix with no rigid-body modes. Elsewhere it +is *safe* but not always *effective*. Where rigid-body modes and a +near-singular mass coincide, QZ may return the theoretically real zero +modes as complex-conjugate pairs, and where those land differs between +LAPACK builds; inside the requested window the ordering cannot be +verified and the retry declines. Declining is deliberate — a guard added +to stop a silent wrong answer must not be able to introduce one — and +``max_residual`` still reports the problem. Note on the user-spec mode choice: ``eigsh(..., sigma=0, mode='buckling')`` reduces to ``OP = K^-1 K = I`` for ``sigma=0``, @@ -492,10 +504,6 @@ def _build_diagnostics( ) -# A mode whose eigenvalue is below this fraction of the largest returned -# one is a rigid-body mode: a free-free floating platform has up to six, -# and an unrestrained DOF (a symmetric column's yaw) gives an exactly -# zero one. # Above this the mass matrix is ill-conditioned enough for the Cholesky # reduction to be the credible culprit; below it, something else in the # pencil is. @@ -593,18 +601,19 @@ def _compare_candidate_modes( the backward error is the honest outcome. Together the two verdicts give the guarantee the caller relies on: a - mode that was acceptable can only end up above a tenth of the - threshold by having *improved*, never as collateral. + mode that was acceptable can only end up above the regression floor + by having *improved*, never as collateral. The comparison has to be **per mode**, not on the two maxima. A rigid-body mode's backward error is a ratio of two near-zero quantities and reads ~1 in *both* candidates however exact each is, so it sets a floor under the alternative's maximum: with one present, - ``max(alt_r)`` stays near 1 and no amount of improvement elsewhere - can drive it below a tenth of ``max(sym_r)`` unless the symmetric - solve is worse than ~10. A free-free model with a genuinely corrupted - elastic mode at a backward error of ~0.8 would sail through, which is - exactly the breakdown this guard exists to catch. + ``max(alt_r)`` stays near 1, and no improvement elsewhere can drive + it below the required fraction of ``max(sym_r)`` unless the symmetric + solve is worse still by that same fraction inverted. A free-free + model with a genuinely corrupted elastic mode at a backward error of + ~0.8 would sail through, which is exactly the breakdown this guard + exists to catch. Comparing mode by mode removes the floor: the rigid modes contribute ~1 against ~1 and register as no improvement, while a corrupted diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index fab6634..5bf1b1f 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -1,13 +1,18 @@ """The symmetric eigensolvers degrade silently on a near-singular mass matrix, and the solver has to notice. -Both ``scipy.linalg.eigh`` and ``scipy.sparse.linalg.eigsh`` reduce -``K x = lambda M x`` through a Cholesky factor of the mass matrix. When -that matrix is nearly singular — a very light beam carrying a very heavy -lump — the reduction loses accuracy, and LAPACK returns confidently -wrong low modes rather than raising. On the case pinned below the dense -symmetric path reported 0.103 Hz against a true 0.0436 Hz, a factor of -2.4, with no error and no warning. +``scipy.linalg.eigh`` reduces ``K x = lambda M x`` through a Cholesky +factor of the **mass** matrix. When that matrix is nearly singular — a +very light beam carrying a very heavy lump — the reduction loses +accuracy, and LAPACK returns confidently wrong low modes rather than +raising. On the case pinned below the dense symmetric path reported +0.103 Hz against a true 0.0436 Hz, a factor of 2.4, with no error and no +warning. + +The sparse path is exempt and must stay exempt: ``eigsh(sigma=0, +mode='normal')`` factorises ``K`` instead, so a near-singular ``M`` does +not degrade it, and its ``which="LM"`` window selects a different set of +modes that must never be compared against the retry's by index. The guard is the backward error ``||K x - lambda M x|| / ||K x||``, and almost every test here exists because some reading of it turned out to From 430ec073dcdd4ed50043b39edb0752296ea07e60 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 19:33:02 +0900 Subject: [PATCH 19/28] docs: state each fact once and point at it, rather than restating it Two more from Codex, both in passages my last commit rewrote. That rewrite fixed the module docstring and left the inline comment beside the guard still attributing the mass-matrix failure to both symmetric paths, and a test docstring still recording two abandoned classifiers where the docstring above it now says three. So the rewrite was not the structural fix I claimed. Restating the same reasoning in the module docstring, the inline comments and the test docstrings means three copies that drift apart, and fixing whichever copy is reported keeps the other two stale. That is the actual mechanism behind eight of the findings on this pull request. The reasoning now lives in one place, the module docstring, and the other sites carry the operative fact plus a pointer. Less to keep in step, and a reader who wants the why has one place to look rather than three that may disagree. --- src/pybmodes/fem/solver.py | 40 +++++++------------------- tests/fem/test_ill_conditioned_mass.py | 28 +++++++----------- 2 files changed, 21 insertions(+), 47 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index d22adbb..2916169 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -329,15 +329,12 @@ def solve_modes( _normalize_columns_l2(eigvecs) - # Accuracy guarantee for the symmetric paths. Both ``eigh`` and - # ``eigsh`` reduce ``K x = λ M x`` through a Cholesky factor of one of - # the matrices, and that reduction degrades once the factored matrix - # is nearly singular — a very light beam carrying a very heavy lump - # does exactly that to ``M``. The failure is silent: LAPACK returns - # confidently wrong low modes rather than raising. The backward error - # catches it (healthy solves sit at ~1e-4 or below, degraded ones - # above 1), and the general path, which factorises neither matrix, - # stays exact there. + # The residual retry — see the module docstring for the rule and for + # why each of its clauses exists. In short: dense ``eigh`` reduces + # through a Cholesky factor of the *mass* matrix and fails silently + # when that is nearly singular, and the backward error is what + # catches it. + # # The matrices the returned modes actually solve. Both symmetric # paths symmetrise internally, so for them the diagnostics — and the # retry decision below — have to be measured against that pair, not @@ -346,26 +343,11 @@ def solve_modes( res_k, res_m = (0.5 * (gk + gk.T), 0.5 * (gm + gm.T)) if sym else (gk, gm) residual_fallback = False - # Only the *dense* symmetric path is retried, and that is a statement - # about which matrix each routine factorises rather than a - # convenience. ``eigh`` reduces through a Cholesky factor of the mass - # matrix, which is the one this guard exists for. ``eigsh(sigma=0, - # mode='normal')`` factorises ``K`` instead, so a near-singular ``M`` - # does not degrade it — the mesh sweep that motivated this work - # returns correct frequencies on exactly the meshes large enough to - # take the sparse path. - # - # Excluding it also removes a mismatch that would otherwise need - # separate handling: ``which="LM"`` on ``OP = K^-1 M`` selects the - # modes nearest zero *in magnitude*, while the retry selects the - # algebraically smallest. With negative eigenvalues present — a - # post-buckling ``run(gravity=...)`` column — those are different - # sets, and a per-index comparison between them would be pairing - # unrelated modes. - # - # The size ceiling matters only because a sparse solve that fails to - # converge falls back to the dense path at *any* size, where an - # unbounded ``eig`` could spend minutes on a result already in hand. + # Dense symmetric only — ``eigsh`` factorises ``K``, so it does not + # have this failure, and its mode window is a different set that must + # not be index-compared. Size-capped because a sparse solve that + # fails to converge falls back to dense at *any* size. Both reasons + # in full in the module docstring. if ( sym and path == "dense_symmetric" diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 5bf1b1f..850c39c 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -228,20 +228,12 @@ class TestRigidBodyModesAreNotMistakenForBreakdown: """A free-free model's zero-frequency modes must not corrupt the result. For a rigid-body mode ``K x ~ 0`` and ``lambda ~ 0``, so the relative - residual is a ratio of two near-zero quantities and reads ~1 however - exact the eigenpair is. Two attempts to *identify* such modes and - exclude them both failed — an eigenvalue-relative cutoff takes a - rigid-only subset's own noise as its scale, and a strain-relative one - cannot tell a rigid mode from a genuinely soft one (the 0.08 Hz lump - mode of a 1e10 N.m^2 beam carries less strain than a floating - platform's rigid modes do, and excluding it blinded the guard to a - case it had caught). - - So they are not identified at all. Both candidate solves are measured - the same way and the retry needs a decisive win, so a mode the metric - cannot speak to says the same nothing twice. The retry also preserves - zero eigenvalues, which is what makes a false positive merely wasteful - instead of destructive. + residual divides one roundoff quantity by another and its value is + arbitrary. Such modes are deliberately *not* identified — three + attempts to do so failed, which the module docstring of + :mod:`pybmodes.fem.solver` records. They are neutralised instead, by + judging the size of the win and by preserving zero eigenvalues so a + false positive is wasteful rather than destructive. """ def _free_free_with_a_zero_mode(self): @@ -312,11 +304,11 @@ def _six_rigid_dofs(self): @pytest.mark.parametrize("n_modes", [1, 3, 6]) def test_a_rigid_only_subset_keeps_its_modes(self, n_modes): - """The case that broke both classification attempts. + """The case that broke the eigenvalue-scale classifier. - Every requested mode is rigid-body, so the metric reads ~1 on all - of them and no reference scale drawn from the subset can say - otherwise. The retry may well run; what matters is that it cannot + Every requested mode is rigid-body, so no reference scale drawn + from the subset can say which of them is meaningful. The retry + may well run; what matters is that it cannot take modes away, because it now preserves zero eigenvalues and has to win decisively to be accepted at all. """ From acb57844de8cc36d0bc979b7997db33e92acd292 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 19:45:59 +0900 Subject: [PATCH 20/28] perf: stop the residual check allocating a copy of the matrices Codex P2. Measuring against the symmetrised pair was done by building it, unconditionally, on every solve. A large sparse solve therefore paid for two dense ngd-square allocations it never used, on exactly the systems the sparse path exists to keep cheap. I noted this cost in my own audit and judged it acceptable, which was wrong: at ngd = 1500 it is a 54 MB peak against 0.3 MB, a factor of 187. 0.5 (A + A.T) v equals 0.5 (A v + A.T v), so the basis change is now carried as a flag and applied through two thin products. Identical results to 1e-12, pinned by a test, along with one asserting the flag is not inert and one bounding the allocation below a quarter of a single matrix. The symmetrised pair is still built in one place, the retry, because eig needs matrices rather than products. That branch is rare and already size-capped. --- src/pybmodes/fem/solver.py | 57 ++++++++++++++++----- tests/fem/test_ill_conditioned_mass.py | 70 ++++++++++++++++++++++++++ 2 files changed, 114 insertions(+), 13 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 2916169..ec4ffdd 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -335,13 +335,12 @@ def solve_modes( # when that is nearly singular, and the backward error is what # catches it. # - # The matrices the returned modes actually solve. Both symmetric - # paths symmetrise internally, so for them the diagnostics — and the - # retry decision below — have to be measured against that pair, not - # against the raw one. Reporting the raw backward error would flag a - # correct solve as defective in telemetry meant to be auditable. - res_k, res_m = (0.5 * (gk + gk.T), 0.5 * (gm + gm.T)) if sym else (gk, gm) - + # Everything below measures against the matrices the returned modes + # actually solve — the symmetrised pair on a symmetric path, since + # both symmetrise internally. That is passed as a flag rather than by + # building the pair here: materialising it costs two dense ngd-square + # allocations, and a large sparse solve would pay for them on every + # call without ever needing them. residual_fallback = False # Dense symmetric only — ``eigsh`` factorises ``K``, so it does not # have this failure, and its mode window is a different set that must @@ -362,7 +361,11 @@ def solve_modes( # above the threshold, and ``eig`` on those same unsymmetrised # matrices would "win decisively" purely by answering a different # question — replacing a correct spectrum with the skew's. - gk_s, gm_s = res_k, res_m + # Only here is the symmetrised pair actually built: the retry + # feeds it to ``eig``, which needs matrices rather than products. + # Bounded by the size cap above, and this branch is rare. + gk_s = 0.5 * (gk + gk.T) + gm_s = 0.5 * (gm + gm.T) sym_r = _modal_residuals(gk_s, gm_s, eigvals, eigvecs) if sym_r.size and float(sym_r.max()) > _SOLVER_OPTIONS.residual_retry_threshold: try: @@ -448,7 +451,7 @@ def solve_modes( return eigvals, eigvecs diagnostics = _build_diagnostics( - res_k, res_m, eigvals, eigvecs, path=path, symmetric=sym, + gk, gm, eigvals, eigvecs, path=path, symmetric=sym, n_requested=n_modes, sparse_fallback=sparse_fallback, fallback_reason=fallback_reason, residual_fallback=residual_fallback, @@ -469,8 +472,17 @@ def _build_diagnostics( fallback_reason: str | None, residual_fallback: bool = False, ) -> SolverDiagnostics: - """Assemble a :class:`SolverDiagnostics` for a completed solve.""" - residuals = _modal_residuals(gk, gm, eigvals, eigvecs) + """Assemble a :class:`SolverDiagnostics` for a completed solve. + + ``symmetric`` selects the basis the residuals are measured on: a + symmetric path solved the symmetrised pair, so charging its modes for + the skew it was told to discard would report a correct solve as + defective. After a residual retry the modes came from ``eig`` on that + same symmetrised pair, so the basis is unchanged. + """ + residuals = _modal_residuals( + gk, gm, eigvals, eigvecs, symmetrise=symmetric, + ) cond = _mass_matrix_cond(gm, path) return SolverDiagnostics( path=path, @@ -651,18 +663,37 @@ def _compare_candidate_modes( return improved, regressed +def _apply(a: np.ndarray, v: np.ndarray, symmetrise: bool) -> np.ndarray: + """``A v``, or ``sym(A) v`` without ever forming ``sym(A)``. + + ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``, and the right-hand side + costs two thin products rather than a dense ngd-square allocation. + That matters because the residuals are computed on every solve that + asks for diagnostics, including the large sparse ones the dense + allocation would hurt most. + """ + if not symmetrise: + return np.asarray(a @ v) + return np.asarray(0.5 * (a @ v + a.T @ v)) + + def _modal_residuals( gk: np.ndarray, gm: np.ndarray, eigvals: np.ndarray, eigvecs: np.ndarray, + *, symmetrise: bool = False, ) -> np.ndarray: """Per-mode relative backward error ``||K x - λ M x|| / ||K x||``. The honest health metric for a generalised modal solve. Cheap (matrix-times-thin-matrix), so computed for every path. + + ``symmetrise`` measures against ``sym(A)`` instead, which is what a + symmetric path actually solved — see :func:`_apply` for why that is + done through the products rather than by building the pair. """ if eigvecs.size == 0: return np.empty(0, dtype=float) - kx = gk @ eigvecs # (ngd, k) - mx = gm @ eigvecs + kx = _apply(gk, eigvecs, symmetrise) # (ngd, k) + mx = _apply(gm, eigvecs, symmetrise) num = np.linalg.norm(kx - mx * eigvals[np.newaxis, :], axis=0) den = np.linalg.norm(kx, axis=0) return np.asarray(num / np.where(den > 0.0, den, 1.0), dtype=float) diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 850c39c..7d9d7c8 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -1130,6 +1130,76 @@ def test_a_genuine_shortfall_is_still_reported(self): solve_modes(gk, gm, n_modes=2) +class TestTheGuardDoesNotTaxEverySolve: + """Residuals are computed on every solve that asks for diagnostics, + so measuring them must not allocate a copy of the matrices. + + ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``, and only the second form + avoids a dense ngd-square pair. At ngd = 1500 that is a 54 MB peak + against 0.3 MB, and it falls on the large sparse solves the dense + copy would hurt most. + """ + + def _pair(self, n): + rng = np.random.default_rng(3) + a = rng.normal(size=(n, n)) + gm = a @ a.T + n * np.eye(n) + b = rng.normal(size=(n, n)) + gk = b @ b.T + return 0.5 * (gk + gk.T), 0.5 * (gm + gm.T) + + def test_symmetrised_products_match_materialising(self): + from pybmodes.fem.solver import _modal_residuals + + gk, gm = self._pair(40) + rng = np.random.default_rng(5) + v = np.linalg.qr(rng.normal(size=(40, 4)))[0] + w = np.linspace(1.0, 2.0, 4) + ks, ms = 0.5 * (gk + gk.T), 0.5 * (gm + gm.T) + assert np.allclose( + _modal_residuals(gk, gm, w, v, symmetrise=True), + _modal_residuals(ks, ms, w, v), + rtol=1.0e-12, atol=1.0e-15, + ) + + def test_an_asymmetric_pair_is_measured_unsymmetrised(self): + """The flag must actually change the basis, not be inert.""" + from pybmodes.fem.solver import _modal_residuals + + gk = np.array([[1.0, 0.9], [0.0, 2.0]]) + gm = np.eye(2) + rng = np.random.default_rng(9) + v = np.linalg.qr(rng.normal(size=(2, 2)))[0] + w = np.array([1.0, 2.0]) + raw = _modal_residuals(gk, gm, w, v) + sym = _modal_residuals(gk, gm, w, v, symmetrise=True) + assert not np.allclose(raw, sym) + + def test_measuring_does_not_allocate_a_matrix_copy(self): + import tracemalloc + + from pybmodes.fem.solver import _modal_residuals + + n = 400 + gk, gm = self._pair(n) + rng = np.random.default_rng(7) + v = np.linalg.qr(rng.normal(size=(n, 5)))[0] + w = np.linspace(1.0, 2.0, 5) + + tracemalloc.start() + base = tracemalloc.get_traced_memory()[0] + _modal_residuals(gk, gm, w, v, symmetrise=True) + peak = tracemalloc.get_traced_memory()[1] - base + tracemalloc.stop() + + # Well under a single dense copy, which is what materialising the + # symmetrised pair would have cost twice over. + assert peak < 0.25 * gk.nbytes, ( + f"peak {peak / 1e6:.1f} MB against a {gk.nbytes / 1e6:.1f} MB " + f"matrix — the symmetrised pair is being materialised" + ) + + class TestDiagnosticsContract: def test_residual_fallback_defaults_to_false(self): gk, gm = _cantilever_with_tip_lump(13, REALISTIC) From dfbb96c99d9cc748652a4c97982b05c27cfddc8e Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 20:05:34 +0900 Subject: [PATCH 21/28] perf: defer the symmetrised pair until the retry threshold is tripped Codex P2, and the same mistake as the previous commit one level in. I switched the diagnostics to thin products and left the sibling call site a hundred lines above still building the pair, before the threshold that decides whether the retry runs at all has even been tested. Every healthy dense solve up to the size cap paid two dense allocations to find out it did not need them. The pre-threshold check now goes through the products too. The pair is built inside the branch, where eig genuinely needs matrices. The test I first wrote for this asserted the whole solve allocates less than half a matrix, which is impossible and would have been a bad test even if it passed: eigh takes matrices, so _solve_dense_symmetric legitimately builds the pair it needs, and that dominates any peak measurement. It now asserts by identity that the check reads the caller arrays rather than a copy, which is the property actually at stake. --- src/pybmodes/fem/solver.py | 34 ++++++++++++++------------ tests/fem/test_ill_conditioned_mass.py | 32 ++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 15 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index ec4ffdd..6e4af76 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -352,22 +352,26 @@ def solve_modes( and path == "dense_symmetric" and ngd <= _SOLVER_OPTIONS.residual_retry_max_ndof ): - # Measure — and retry — against the matrices the symmetric paths - # actually solved. Both symmetrise internally, and the accepted - # skew is only guaranteed small relative to ``max|K|``: in a model - # with a wide dynamic range it can still be large relative to a - # soft mode's own eigenvalue. Judging an exact symmetric solve - # against the unsymmetrised matrices would then show a residual - # above the threshold, and ``eig`` on those same unsymmetrised - # matrices would "win decisively" purely by answering a different - # question — replacing a correct spectrum with the skew's. - # Only here is the symmetrised pair actually built: the retry - # feeds it to ``eig``, which needs matrices rather than products. - # Bounded by the size cap above, and this branch is rare. - gk_s = 0.5 * (gk + gk.T) - gm_s = 0.5 * (gm + gm.T) - sym_r = _modal_residuals(gk_s, gm_s, eigvals, eigvecs) + # Measured against the matrices the symmetric paths actually + # solved. Both symmetrise internally, and the accepted skew is + # only guaranteed small relative to ``max|K|``: in a model with a + # wide dynamic range it can still be large relative to a soft + # mode's own eigenvalue. Judging an exact symmetric solve against + # the unsymmetrised matrices would then show a residual above the + # threshold, and ``eig`` on those same unsymmetrised matrices + # would "win decisively" purely by answering a different question + # — replacing a correct spectrum with the skew's. + # + # Through the products, not the pair: this runs on every eligible + # solve, healthy ones included, and the threshold below has not + # been tested yet. + sym_r = _modal_residuals(gk, gm, eigvals, eigvecs, symmetrise=True) if sym_r.size and float(sym_r.max()) > _SOLVER_OPTIONS.residual_retry_threshold: + # Now the pair is worth building: ``eig`` needs matrices + # rather than products. Bounded by the size cap above, and + # this branch is rare. + gk_s = 0.5 * (gk + gk.T) + gm_s = 0.5 * (gm + gm.T) try: alt_vals, alt_vecs, ordering_sound = _general_spectrum_for_retry( gk_s, gm_s, n_modes, diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 7d9d7c8..faa738e 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -1175,6 +1175,38 @@ def test_an_asymmetric_pair_is_measured_unsymmetrised(self): sym = _modal_residuals(gk, gm, w, v, symmetrise=True) assert not np.allclose(raw, sym) + def test_the_threshold_check_reads_the_caller_s_own_matrices( + self, monkeypatch, + ): + """The check runs on every eligible solve, healthy ones included, + before the threshold has been tested. It must therefore read the + arrays it was given rather than a symmetrised copy of them. + + Asserted by identity rather than by measuring memory: the solve + legitimately allocates elsewhere — ``eigh`` takes matrices, so + ``_solve_dense_symmetric`` builds the pair it needs — which would + swamp any peak-usage threshold and make the test meaningless. + """ + import pybmodes.fem.solver as solvermod + + seen = [] + real = solvermod._modal_residuals + + def spy(k, m, vals, vecs, *, symmetrise=False): + seen.append((k is gk, m is gm, symmetrise)) + return real(k, m, vals, vecs, symmetrise=symmetrise) + + monkeypatch.setattr(solvermod, "_modal_residuals", spy) + gk, gm = self._pair(60) + _v, _x, diag = solve_modes(gk, gm, n_modes=4, return_diagnostics=True) + + assert diag.path == "dense_symmetric" + assert diag.residual_fallback is False + assert seen, "the threshold check did not run" + # Every call took the caller's arrays and asked for the + # symmetrised basis, rather than being handed a built pair. + assert all(is_gk and is_gm and sym for is_gk, is_gm, sym in seen), seen + def test_measuring_does_not_allocate_a_matrix_copy(self): import tracemalloc From e93347f6f47f1ac36326eaa4f046f5f48a77fddc Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 20:14:44 +0900 Subject: [PATCH 22/28] fix: choose the residual product route by block width The two forms of sym(A) @ V peak at 3nk and 2n^2 bytes of temporaries. Splitting the product was only ever cheaper for a narrow block, but n_modes=None is the public default and returns the whole spectrum, so the wide case was paying for two full matmuls to avoid one copy. Select the route by width. The crossover is measured rather than derived: a flop-count estimate put it at k = n/3 when the allocator puts it at 2n/3, which would have taken the dearer route across a third of the range. --- src/pybmodes/fem/solver.py | 32 ++++++++++++++---- tests/fem/test_ill_conditioned_mass.py | 46 ++++++++++++++++++++++++++ 2 files changed, 71 insertions(+), 7 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 6e4af76..4200095 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -667,17 +667,35 @@ def _compare_candidate_modes( return improved, regressed -def _apply(a: np.ndarray, v: np.ndarray, symmetrise: bool) -> np.ndarray: - """``A v``, or ``sym(A) v`` without ever forming ``sym(A)``. +# The two routes to ``sym(A) v`` peak at ``3 n k`` and ``2 n^2`` bytes of +# temporaries, so they cross over at ``k = 2 n / 3`` — measured, not +# derived: a first estimate of ``n / 3`` was wrong by a factor of two and +# would have taken the more expensive route across a third of the range. +_THIN_BLOCK_NUM, _THIN_BLOCK_DEN = 3, 2 + - ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``, and the right-hand side - costs two thin products rather than a dense ngd-square allocation. - That matters because the residuals are computed on every solve that - asks for diagnostics, including the large sparse ones the dense - allocation would hurt most. +def _apply(a: np.ndarray, v: np.ndarray, symmetrise: bool) -> np.ndarray: + """``A v``, or ``sym(A) v`` by whichever route is cheaper. + + ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``. The right-hand side avoids + a dense ngd-square allocation, which is what a modal solve usually + wants: residuals are computed on every solve that asks for + diagnostics, and a handful of modes out of a few thousand DOFs makes + those products very thin. + + It is not free, though, and the assumption fails at the other end. + ``n_modes=None`` is the public default and returns the whole + spectrum, making ``v`` square — the "products" are then two full + matrix multiplies whose temporaries exceed the single copy they were + avoiding. So the route is chosen by the block width rather than + assumed. """ if not symmetrise: return np.asarray(a @ v) + if v.ndim > 1 and ( + v.shape[1] * _THIN_BLOCK_NUM >= a.shape[0] * _THIN_BLOCK_DEN + ): + return np.asarray(0.5 * (a + a.T) @ v) return np.asarray(0.5 * (a @ v + a.T @ v)) diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index faa738e..6714848 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -1207,6 +1207,52 @@ def spy(k, m, vals, vecs, *, symmetrise=False): # symmetrised basis, rather than being handed a built pair. assert all(is_gk and is_gm and sym for is_gk, is_gm, sym in seen), seen + @pytest.mark.parametrize("k_over_n", [0.01, 0.2, 0.5, 0.9, 1.0]) + def test_whichever_route_is_taken_is_the_cheaper_one(self, k_over_n): + """The products are only thin while the block is. ``n_modes=None`` + is the public default and returns the whole spectrum, at which + point two full matmuls cost more than the single copy they were + avoiding, so the route is chosen by width. + + The crossover is measured rather than derived — a flop-count + estimate put it at ``k = n/3`` when it is ``2n/3``, which would + have taken the dearer route across a third of the range. + """ + import tracemalloc + + from pybmodes.fem.solver import _apply + + n = 400 + rng = np.random.default_rng(4) + a = rng.normal(size=(n, n)) + v = rng.normal(size=(n, max(1, round(k_over_n * n)))) + + def peak(fn): + tracemalloc.start() + base = tracemalloc.get_traced_memory()[0] + fn() + p = tracemalloc.get_traced_memory()[1] - base + tracemalloc.stop() + return p + + chosen = peak(lambda: _apply(a, v, True)) + products = peak(lambda: 0.5 * (a @ v + a.T @ v)) + materialise = peak(lambda: 0.5 * (a + a.T) @ v) + assert chosen <= 1.05 * min(products, materialise), ( + f"k/n={k_over_n}: took {chosen / 1e6:.2f} MB when " + f"{min(products, materialise) / 1e6:.2f} MB was available" + ) + + @pytest.mark.parametrize("k", [1, 7, 40]) + def test_both_routes_agree_numerically(self, k): + from pybmodes.fem.solver import _apply + + n = 40 + rng = np.random.default_rng(6) + a = rng.normal(size=(n, n)) + v = rng.normal(size=(n, k)) + assert np.allclose(_apply(a, v, True), 0.5 * (a + a.T) @ v) + def test_measuring_does_not_allocate_a_matrix_copy(self): import tracemalloc From cad2de11af3f82b8507d663f4e42000e6f2b9826 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 20:20:06 +0900 Subject: [PATCH 23/28] test: assert the route rule, not a measured allocator peak The width test compared three live tracemalloc peaks. That measures what the allocator did, not what the routes cost: a block freed by an earlier call gets reused by a later one, so the three time-share buffers differently depending on measurement order. It passed on Windows and failed on Linux at k/n = 0.9, where the route taken was correct and only the comparison figure had shifted. Lift the decision into _prefer_materialised and assert it against the cost model, which is deterministic everywhere. Adds explicit cover for the full-width case that prompted the rule and for the thin case it must not disturb, and widens the numerical-agreement test to straddle the crossover. --- src/pybmodes/fem/solver.py | 26 ++++++--- tests/fem/test_ill_conditioned_mass.py | 80 +++++++++++++++----------- 2 files changed, 65 insertions(+), 41 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 4200095..14c90c2 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -671,7 +671,18 @@ def _compare_candidate_modes( # temporaries, so they cross over at ``k = 2 n / 3`` — measured, not # derived: a first estimate of ``n / 3`` was wrong by a factor of two and # would have taken the more expensive route across a third of the range. -_THIN_BLOCK_NUM, _THIN_BLOCK_DEN = 3, 2 +def _prefer_materialised(n_rows: int, n_cols: int) -> bool: + """Is ``sym(A) v`` cheaper built than split into two products? + + The split form peaks at three ``n x k`` temporaries and the built one + at an ``n x n`` copy plus the ``n x k`` result, so they cross over at + ``3 n k = 2 n^2``, i.e. ``k = 2n/3``. + + Measured, not derived. A flop-count estimate put the crossover at + ``n/3``, which would have taken the dearer route across a third of + the range — precisely the middling widths this sees in practice. + """ + return n_cols * 3 >= n_rows * 2 def _apply(a: np.ndarray, v: np.ndarray, symmetrise: bool) -> np.ndarray: @@ -679,8 +690,8 @@ def _apply(a: np.ndarray, v: np.ndarray, symmetrise: bool) -> np.ndarray: ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``. The right-hand side avoids a dense ngd-square allocation, which is what a modal solve usually - wants: residuals are computed on every solve that asks for - diagnostics, and a handful of modes out of a few thousand DOFs makes + wants: residuals are measured on every eligible solve, healthy ones + included, and a handful of modes out of a few thousand DOFs makes those products very thin. It is not free, though, and the assumption fails at the other end. @@ -692,9 +703,7 @@ def _apply(a: np.ndarray, v: np.ndarray, symmetrise: bool) -> np.ndarray: """ if not symmetrise: return np.asarray(a @ v) - if v.ndim > 1 and ( - v.shape[1] * _THIN_BLOCK_NUM >= a.shape[0] * _THIN_BLOCK_DEN - ): + if v.ndim > 1 and _prefer_materialised(a.shape[0], v.shape[1]): return np.asarray(0.5 * (a + a.T) @ v) return np.asarray(0.5 * (a @ v + a.T @ v)) @@ -705,8 +714,9 @@ def _modal_residuals( ) -> np.ndarray: """Per-mode relative backward error ``||K x - λ M x|| / ||K x||``. - The honest health metric for a generalised modal solve. Cheap - (matrix-times-thin-matrix), so computed for every path. + The honest health metric for a generalised modal solve, and cheap + enough to compute on every path — a matrix against a block of + eigenvectors, usually a thin one. ``symmetrise`` measures against ``sym(A)`` instead, which is what a symmetric path actually solved — see :func:`_apply` for why that is diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 6714848..49d7368 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -1131,13 +1131,15 @@ def test_a_genuine_shortfall_is_still_reported(self): class TestTheGuardDoesNotTaxEverySolve: - """Residuals are computed on every solve that asks for diagnostics, - so measuring them must not allocate a copy of the matrices. - - ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``, and only the second form - avoids a dense ngd-square pair. At ngd = 1500 that is a 54 MB peak - against 0.3 MB, and it falls on the large sparse solves the dense - copy would hurt most. + """Residuals are measured on every eligible solve, healthy ones + included, so measuring them must not cost more than it saves. + + ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``, and the two differ only in + what they allocate. For a thin block the second form avoids a dense + ngd-square pair: at ngd = 1500 with six modes that is 0.3 MB against + 54 MB, and it falls on the large sparse solves the copy would hurt + most. For a wide one the comparison inverts, so the route is chosen + rather than assumed. """ def _pair(self, n): @@ -1207,44 +1209,56 @@ def spy(k, m, vals, vecs, *, symmetrise=False): # symmetrised basis, rather than being handed a built pair. assert all(is_gk and is_gm and sym for is_gk, is_gm, sym in seen), seen - @pytest.mark.parametrize("k_over_n", [0.01, 0.2, 0.5, 0.9, 1.0]) + @pytest.mark.parametrize("k_over_n", [0.01, 0.2, 0.5, 0.67, 0.9, 1.0]) def test_whichever_route_is_taken_is_the_cheaper_one(self, k_over_n): """The products are only thin while the block is. ``n_modes=None`` is the public default and returns the whole spectrum, at which point two full matmuls cost more than the single copy they were avoiding, so the route is chosen by width. - The crossover is measured rather than derived — a flop-count - estimate put it at ``k = n/3`` when it is ``2n/3``, which would - have taken the dearer route across a third of the range. + Asserted against the cost model rather than against a measured + peak. ``tracemalloc`` reports what the *allocator* did, and a + freed block from an earlier call can be reused by a later one, so + the same three routes time-share buffers differently depending on + the order they are measured in and on the platform. An earlier + version of this test compared three live measurements and failed + on Linux while passing on Windows, having measured reuse rather + than cost. + + The model itself is measured, in ``scripts`` runs recorded on the + pull request: split form ``3nk``, built form ``2n^2``, crossing + at ``k = 2n/3``. """ - import tracemalloc - - from pybmodes.fem.solver import _apply + from pybmodes.fem.solver import _prefer_materialised n = 400 - rng = np.random.default_rng(4) - a = rng.normal(size=(n, n)) - v = rng.normal(size=(n, max(1, round(k_over_n * n)))) - - def peak(fn): - tracemalloc.start() - base = tracemalloc.get_traced_memory()[0] - fn() - p = tracemalloc.get_traced_memory()[1] - base - tracemalloc.stop() - return p - - chosen = peak(lambda: _apply(a, v, True)) - products = peak(lambda: 0.5 * (a @ v + a.T @ v)) - materialise = peak(lambda: 0.5 * (a + a.T) @ v) - assert chosen <= 1.05 * min(products, materialise), ( - f"k/n={k_over_n}: took {chosen / 1e6:.2f} MB when " - f"{min(products, materialise) / 1e6:.2f} MB was available" + k = max(1, round(k_over_n * n)) + split, built = 3 * n * k, 2 * n * n + assert _prefer_materialised(n, k) == (built <= split), ( + f"k/n={k_over_n}: chose the " + f"{'built' if _prefer_materialised(n, k) else 'split'} form " + f"when split costs {split} and built costs {built}" ) - @pytest.mark.parametrize("k", [1, 7, 40]) + def test_a_full_width_block_is_not_treated_as_thin(self): + """The case that prompted the rule: ``n_modes=None`` returns the + whole spectrum, so ``v`` is square and the products are two full + matrix multiplies.""" + from pybmodes.fem.solver import _prefer_materialised + + assert _prefer_materialised(400, 400) + assert _prefer_materialised(2000, 2000) + + def test_a_few_modes_out_of_many_dofs_stays_thin(self): + from pybmodes.fem.solver import _prefer_materialised + + assert not _prefer_materialised(1500, 6) + assert not _prefer_materialised(400, 4) + + @pytest.mark.parametrize("k", [1, 7, 26, 27, 40]) def test_both_routes_agree_numerically(self, k): + """Whichever route is taken, the answer is ``sym(A) v``. ``k`` + spans both sides of the crossover at ``2n/3 = 26.7``.""" from pybmodes.fem.solver import _apply n = 40 From a34fd66c6facfb9fdbad03adea0bb5ea0461b0d3 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 20:22:12 +0900 Subject: [PATCH 24/28] docs: point the width test at the predicate rather than restating it --- tests/fem/test_ill_conditioned_mass.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 49d7368..e5bf127 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -1225,9 +1225,8 @@ def test_whichever_route_is_taken_is_the_cheaper_one(self, k_over_n): on Linux while passing on Windows, having measured reuse rather than cost. - The model itself is measured, in ``scripts`` runs recorded on the - pull request: split form ``3nk``, built form ``2n^2``, crossing - at ``k = 2n/3``. + The cost model asserted against — ``3nk`` split, ``2n^2`` built — + and where it came from are on :func:`_prefer_materialised`. """ from pybmodes.fem.solver import _prefer_materialised From 96040ec551a1934437f6bf65cecf2430b023db52 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 12 Aug 2026 20:23:46 +0900 Subject: [PATCH 25/28] docs: state the regression guarantee as the code actually makes it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The refusal is on a mode getting strictly worse, so an accepted retry proves an acceptable mode did not get worse — not that it improved. Below the regression floor it may move either way, which is what test_a_mode_that_worsens_but_stays_acceptable_is_not_a_regression pins. Say that rather than the slightly stronger claim. --- CHANGELOG.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4514d32..f55ab4a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -56,8 +56,10 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). one mode while pushing a previously acceptable one above the failure threshold is refused — it would hand back a new bad mode in place of an old one. The guarantee is one-sided and precise: a mode that was - acceptable can end up above the regression floor only by having - improved, never as collateral of another mode's rescue. A mode + acceptable can end up above the regression floor only by not having got + worse, never as collateral of another mode's rescue. Below that floor + it is free to move either way, which is deliberate — a residual already + that small is not a claim about accuracy worth defending. A mode already failing carries no verdict either way — above the threshold neither candidate is trustworthy, and a rigid-body mode, whose residual divides one roundoff quantity by another and has been measured at 12.4 From 022b06cc91b1ab4988c680aa52c0632011faff26 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Thu, 13 Aug 2026 11:27:11 +0900 Subject: [PATCH 26/28] perf: sweep the modal residuals in column blocks Routing the product form by width bounded _apply but not its caller. _modal_residuals held kx while building mx, then both while forming mx * eigvals and the difference, so a full-spectrum request kept four ngd-square arrays live purely to measure: 128 MB at the 2000-DOF retry cap, measured. The residual is per mode, so no step needs every mode at once. Sweep in column blocks and accumulate the two norm vectors. Peak is then set by the block width rather than by an argument the caller chooses, and the thin case stays a single pass computed exactly as before. At ngd = 2000 with every mode requested the peak drops from 128 MB to 10.3 MB, results agreeing to 1e-12. The same measurement shows no change at all for six modes. This also removes the need for the width test in _apply, since blocking guarantees it never sees a wide block. --- src/pybmodes/fem/solver.py | 73 ++++++------- tests/fem/test_ill_conditioned_mass.py | 135 +++++++++++++++++-------- 2 files changed, 128 insertions(+), 80 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 14c90c2..2b34ca6 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -671,40 +671,25 @@ def _compare_candidate_modes( # temporaries, so they cross over at ``k = 2 n / 3`` — measured, not # derived: a first estimate of ``n / 3`` was wrong by a factor of two and # would have taken the more expensive route across a third of the range. -def _prefer_materialised(n_rows: int, n_cols: int) -> bool: - """Is ``sym(A) v`` cheaper built than split into two products? - - The split form peaks at three ``n x k`` temporaries and the built one - at an ``n x n`` copy plus the ``n x k`` result, so they cross over at - ``3 n k = 2 n^2``, i.e. ``k = 2n/3``. - - Measured, not derived. A flop-count estimate put the crossover at - ``n/3``, which would have taken the dearer route across a third of - the range — precisely the middling widths this sees in practice. - """ - return n_cols * 3 >= n_rows * 2 +# Columns per pass of the residual sweep. Every temporary in the sweep +# is ``ngd x`` this, so the peak is bounded by it rather than by the +# number of modes the caller asked for. Wide enough that the BLAS call +# still amortises its own overhead. +_RESIDUAL_BLOCK = 128 def _apply(a: np.ndarray, v: np.ndarray, symmetrise: bool) -> np.ndarray: - """``A v``, or ``sym(A) v`` by whichever route is cheaper. - - ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``. The right-hand side avoids - a dense ngd-square allocation, which is what a modal solve usually - wants: residuals are measured on every eligible solve, healthy ones - included, and a handful of modes out of a few thousand DOFs makes - those products very thin. - - It is not free, though, and the assumption fails at the other end. - ``n_modes=None`` is the public default and returns the whole - spectrum, making ``v`` square — the "products" are then two full - matrix multiplies whose temporaries exceed the single copy they were - avoiding. So the route is chosen by the block width rather than - assumed. + """``A v``, or ``sym(A) v`` without building ``sym(A)``. + + ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``, and the right-hand side + avoids a dense ngd-square allocation. That only holds while ``v`` is + narrow — at full width the two products cost more than the copy they + were avoiding — so this is called on one column block at a time and + never sees a wide ``v``. The blocking, not a width test here, is what + keeps the choice safe; see :func:`_modal_residuals`. """ if not symmetrise: return np.asarray(a @ v) - if v.ndim > 1 and _prefer_materialised(a.shape[0], v.shape[1]): - return np.asarray(0.5 * (a + a.T) @ v) return np.asarray(0.5 * (a @ v + a.T @ v)) @@ -715,19 +700,35 @@ def _modal_residuals( """Per-mode relative backward error ``||K x - λ M x|| / ||K x||``. The honest health metric for a generalised modal solve, and cheap - enough to compute on every path — a matrix against a block of - eigenvectors, usually a thin one. + enough to compute on every path. ``symmetrise`` measures against ``sym(A)`` instead, which is what a - symmetric path actually solved — see :func:`_apply` for why that is - done through the products rather than by building the pair. + symmetric path actually solved. + + **Swept in column blocks.** The residual is per mode, so no step + needs every mode present at once, and holding them all would tie the + peak to a number the caller chooses: ``n_modes=None`` is the public + default and returns the whole spectrum, which at the retry size cap + would put four ngd-square arrays live at once — over 100 MB — purely + to measure. Blocking bounds every temporary at ``ngd`` by + :data:`_RESIDUAL_BLOCK` instead, and leaves the common thin case + (a handful of modes out of a few thousand DOFs) in a single pass, + computed exactly as before. """ if eigvecs.size == 0: return np.empty(0, dtype=float) - kx = _apply(gk, eigvecs, symmetrise) # (ngd, k) - mx = _apply(gm, eigvecs, symmetrise) - num = np.linalg.norm(kx - mx * eigvals[np.newaxis, :], axis=0) - den = np.linalg.norm(kx, axis=0) + n_modes_out = eigvecs.shape[1] + num = np.empty(n_modes_out, dtype=float) + den = np.empty(n_modes_out, dtype=float) + for lo in range(0, n_modes_out, _RESIDUAL_BLOCK): + hi = min(lo + _RESIDUAL_BLOCK, n_modes_out) + block = eigvecs[:, lo:hi] + kx = _apply(gk, block, symmetrise) # (ngd, <= block) + mx = _apply(gm, block, symmetrise) + den[lo:hi] = np.linalg.norm(kx, axis=0) + num[lo:hi] = np.linalg.norm( + kx - mx * eigvals[np.newaxis, lo:hi], axis=0, + ) return np.asarray(num / np.where(den > 0.0, den, 1.0), dtype=float) diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index e5bf127..a4376f1 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -1134,12 +1134,13 @@ class TestTheGuardDoesNotTaxEverySolve: """Residuals are measured on every eligible solve, healthy ones included, so measuring them must not cost more than it saves. - ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``, and the two differ only in - what they allocate. For a thin block the second form avoids a dense - ngd-square pair: at ngd = 1500 with six modes that is 0.3 MB against - 54 MB, and it falls on the large sparse solves the copy would hurt - most. For a wide one the comparison inverts, so the route is chosen - rather than assumed. + Two things keep it cheap. ``0.5 (A + A.T) v == 0.5 (A v + A.T v)`` + and only the second form avoids a dense ngd-square pair, which + matters most on the large sparse solves the copy would hurt. And the + sweep runs in column blocks, so its peak is set by the block width + rather than by how many modes the caller asked for — without that, + the identity above inverts at full width and ``n_modes=None`` is the + public default. """ def _pair(self, n): @@ -1209,55 +1210,101 @@ def spy(k, m, vals, vecs, *, symmetrise=False): # symmetrised basis, rather than being handed a built pair. assert all(is_gk and is_gm and sym for is_gk, is_gm, sym in seen), seen - @pytest.mark.parametrize("k_over_n", [0.01, 0.2, 0.5, 0.67, 0.9, 1.0]) - def test_whichever_route_is_taken_is_the_cheaper_one(self, k_over_n): - """The products are only thin while the block is. ``n_modes=None`` - is the public default and returns the whole spectrum, at which - point two full matmuls cost more than the single copy they were - avoiding, so the route is chosen by width. - - Asserted against the cost model rather than against a measured - peak. ``tracemalloc`` reports what the *allocator* did, and a - freed block from an earlier call can be reused by a later one, so - the same three routes time-share buffers differently depending on - the order they are measured in and on the platform. An earlier - version of this test compared three live measurements and failed - on Linux while passing on Windows, having measured reuse rather - than cost. - - The cost model asserted against — ``3nk`` split, ``2n^2`` built — - and where it came from are on :func:`_prefer_materialised`. + @pytest.mark.parametrize("k_over_n", [0.01, 0.2, 0.5, 1.0]) + def test_no_pass_of_the_sweep_ever_sees_a_wide_block( + self, k_over_n, monkeypatch, + ): + """The invariant the memory bound rests on. + + ``0.5 (A v + A.T v)`` is only cheaper than building ``sym(A)`` + while ``v`` is narrow, and nothing about the caller's ``n_modes`` + guarantees that — ``None`` is the public default and asks for the + whole spectrum. Blocking is what supplies the guarantee, so this + asserts the width every pass actually sees rather than a measured + peak. + + Measured, not asserted, and recorded here because it is why this + exists: at ngd = 2000 with every mode requested the unblocked + sweep peaked at 128 MB against 10.3 MB blocked. + + Deliberately not a ``tracemalloc`` comparison. That reports what + the allocator did, and a freed block from an earlier call gets + reused by a later one, so competing forms time-share buffers + differently by measurement order and by platform. An earlier + version of this test compared three live peaks and failed on + Linux while passing on Windows, having measured reuse. """ - from pybmodes.fem.solver import _prefer_materialised + from pybmodes.fem import solver - n = 400 + widths = [] + real = solver._apply + monkeypatch.setattr( + solver, "_apply", + lambda a, v, s: (widths.append(v.shape[1]), real(a, v, s))[1], + ) + + n = 200 + gk, gm = self._pair(n) k = max(1, round(k_over_n * n)) - split, built = 3 * n * k, 2 * n * n - assert _prefer_materialised(n, k) == (built <= split), ( - f"k/n={k_over_n}: chose the " - f"{'built' if _prefer_materialised(n, k) else 'split'} form " - f"when split costs {split} and built costs {built}" + rng = np.random.default_rng(12) + v = np.linalg.qr(rng.normal(size=(n, k)))[0] + solver._modal_residuals( + gk, gm, np.linspace(1.0, 2.0, k), v, symmetrise=True, + ) + + assert widths, "the sweep did not run" + assert max(widths) <= solver._RESIDUAL_BLOCK, ( + f"k={k}: a pass saw {max(widths)} columns, above the " + f"{solver._RESIDUAL_BLOCK}-column block" ) + assert sum(widths) == 2 * k, "every mode is measured exactly once" - def test_a_full_width_block_is_not_treated_as_thin(self): - """The case that prompted the rule: ``n_modes=None`` returns the - whole spectrum, so ``v`` is square and the products are two full - matrix multiplies.""" - from pybmodes.fem.solver import _prefer_materialised + def test_a_thin_request_still_runs_in_one_pass(self, monkeypatch): + """Blocking must not tax the common case — a handful of modes out + of a few thousand DOFs is one pass per matrix, exactly the two + products the unblocked form did.""" + from pybmodes.fem import solver - assert _prefer_materialised(400, 400) - assert _prefer_materialised(2000, 2000) + calls = [] + real = solver._apply + monkeypatch.setattr( + solver, "_apply", + lambda a, v, s: (calls.append(v.shape[1]), real(a, v, s))[1], + ) + + gk, gm = self._pair(400) + rng = np.random.default_rng(14) + v = np.linalg.qr(rng.normal(size=(400, 6)))[0] + solver._modal_residuals( + gk, gm, np.linspace(1.0, 2.0, 6), v, symmetrise=True, + ) + assert calls == [6, 6], f"expected one pass per matrix, got {calls}" - def test_a_few_modes_out_of_many_dofs_stays_thin(self): - from pybmodes.fem.solver import _prefer_materialised + @pytest.mark.parametrize("k", [1, 7, 128, 129, 260]) + def test_blocking_does_not_change_the_answer(self, k): + """Straddles the block boundary: one pass, one plus a remainder, + and several.""" + from pybmodes.fem.solver import _apply, _modal_residuals - assert not _prefer_materialised(1500, 6) - assert not _prefer_materialised(400, 4) + n = 300 + gk, gm = self._pair(n) + rng = np.random.default_rng(13) + v = np.linalg.qr(rng.normal(size=(n, k)))[0] + w = np.linspace(1.0, 2.0, k) + + kx, mx = _apply(gk, v, True), _apply(gm, v, True) + den = np.linalg.norm(kx, axis=0) + expected = np.linalg.norm( + kx - mx * w[np.newaxis, :], axis=0, + ) / np.where(den > 0.0, den, 1.0) + assert np.allclose( + _modal_residuals(gk, gm, w, v, symmetrise=True), expected, + rtol=1.0e-12, atol=1.0e-15, + ) @pytest.mark.parametrize("k", [1, 7, 26, 27, 40]) def test_both_routes_agree_numerically(self, k): - """Whichever route is taken, the answer is ``sym(A) v``. ``k`` - spans both sides of the crossover at ``2n/3 = 26.7``.""" + """Whichever route is taken, the answer is ``sym(A) v``.""" from pybmodes.fem.solver import _apply n = 40 From 4dd48df99f12e3c77b829b39748674f956b7c92a Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Thu, 13 Aug 2026 11:31:04 +0900 Subject: [PATCH 27/28] docs: justify the residual block width, and pin it as free 128 was a round number. Measured across widths on a 2000-DOF full spectrum it turns out to sit at the knee: 1.4 s at 16 and 32 columns where per-call BLAS overhead dominates, 0.76 s at 64, 0.42 s at 128, against a 0.38 s floor that 256 and above buy with two to twelve times the peak. Record that where the constant is defined. Add a test that the width changes nothing but the peak, so it can be retuned freely. The residual is per mode, so blocking only partitions independent columns. --- src/pybmodes/fem/solver.py | 10 ++++++++-- tests/fem/test_ill_conditioned_mass.py | 22 ++++++++++++++++++++++ 2 files changed, 30 insertions(+), 2 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 2b34ca6..a4cf80d 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -673,8 +673,14 @@ def _compare_candidate_modes( # would have taken the more expensive route across a third of the range. # Columns per pass of the residual sweep. Every temporary in the sweep # is ``ngd x`` this, so the peak is bounded by it rather than by the -# number of modes the caller asked for. Wide enough that the BLAS call -# still amortises its own overhead. +# number of modes the caller asked for. +# +# Only the peak depends on it — the result does not, since the residual +# is per mode and blocking merely partitions the columns. Chosen at the +# knee of the measured time curve: sweeping the full spectrum of a +# 2000-DOF system took 1.4 s at 16 and 32 columns, where per-call BLAS +# overhead dominates, then 0.76 s at 64 and 0.42 s at 128, against a +# 0.38 s floor that 256 and above buy with two to twelve times the peak. _RESIDUAL_BLOCK = 128 diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index a4376f1..e0d64a8 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -1302,6 +1302,28 @@ def test_blocking_does_not_change_the_answer(self, k): rtol=1.0e-12, atol=1.0e-15, ) + @pytest.mark.parametrize("block", [1, 3, 16, 64, 1000]) + def test_the_block_width_is_a_free_parameter(self, block, monkeypatch): + """It may be tuned for the peak-versus-overhead trade without + touching any answer. The residual is per mode, so blocking only + partitions independent columns — a block width that changed the + result would mean the columns were not independent after all.""" + from pybmodes.fem import solver + + n, k = 120, 37 + gk, gm = self._pair(n) + rng = np.random.default_rng(15) + v = np.linalg.qr(rng.normal(size=(n, k)))[0] + w = np.linspace(1.0, 2.0, k) + + monkeypatch.setattr(solver, "_RESIDUAL_BLOCK", 128) + reference = solver._modal_residuals(gk, gm, w, v, symmetrise=True) + monkeypatch.setattr(solver, "_RESIDUAL_BLOCK", block) + assert np.allclose( + solver._modal_residuals(gk, gm, w, v, symmetrise=True), + reference, rtol=1.0e-12, atol=1.0e-15, + ) + @pytest.mark.parametrize("k", [1, 7, 26, 27, 40]) def test_both_routes_agree_numerically(self, k): """Whichever route is taken, the answer is ``sym(A) v``.""" From 31b4d64b745ef904dfc80448b64b47e138a297cc Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Thu, 13 Aug 2026 11:51:41 +0900 Subject: [PATCH 28/28] fix: restore the width test, which blocking does not subsume I claimed blocking made the width test in _apply unnecessary. That is wrong below ngd = 192: the block is capped at 128 columns, which is narrow against a large system but not a small one, so a full-spectrum request there still hands _apply a block past the 2n/3 crossover. The two mechanisms answer different questions. Blocking bounds the peak as the mode count grows; the width test picks the cheaper route within whatever block it gets. Keep both. Capping the block relative to ngd was the other way to close this and is worse at both ends, measured: at ngd = 100 it splits one BLAS call into two, 206 us against 116 us, and at ngd = 1000 it lifts the peak from 5.1 MB to 21.3 MB. Restoring the width test costs nothing at large ngd, where the block is already narrow and it does not fire, and at ngd = 100 runs the sweep in 116 us against 178 us for the same peak. --- src/pybmodes/fem/solver.py | 33 +++++++++-- tests/fem/test_ill_conditioned_mass.py | 78 ++++++++++++++++++++++---- 2 files changed, 93 insertions(+), 18 deletions(-) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index a4cf80d..8e6a582 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -684,18 +684,39 @@ def _compare_candidate_modes( _RESIDUAL_BLOCK = 128 +def _prefer_materialised(n_rows: int, n_cols: int) -> bool: + """Is ``sym(A) v`` cheaper built than split into two products? + + The split form peaks at three ``n x k`` temporaries and the built one + at an ``n x n`` copy plus the ``n x k`` result, so they cross over at + ``3 n k = 2 n^2``, i.e. ``k = 2n/3``. + + Measured, not derived. A flop-count estimate put the crossover at + ``n/3``, which would have taken the dearer route across a third of + the range. + """ + return n_cols * 3 >= n_rows * 2 + + def _apply(a: np.ndarray, v: np.ndarray, symmetrise: bool) -> np.ndarray: - """``A v``, or ``sym(A) v`` without building ``sym(A)``. + """``A v``, or ``sym(A) v`` by whichever route is cheaper. ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``, and the right-hand side - avoids a dense ngd-square allocation. That only holds while ``v`` is - narrow — at full width the two products cost more than the copy they - were avoiding — so this is called on one column block at a time and - never sees a wide ``v``. The blocking, not a width test here, is what - keeps the choice safe; see :func:`_modal_residuals`. + avoids a dense ngd-square allocation — but only while ``v`` is + narrow. At full width the two products cost more than the copy they + were avoiding, so the route is tested rather than assumed. + + Blocking the caller's sweep does not remove the need for the test. + It bounds the block at :data:`_RESIDUAL_BLOCK` columns, which is + narrow relative to a large ``ngd`` but not to a small one: below + ``ngd = 192`` a full-spectrum request still hands this a block wider + than the crossover. Measured at ``ngd = 100``, testing the width + there runs the sweep in 116 us against 178 us for the same peak. """ if not symmetrise: return np.asarray(a @ v) + if v.ndim > 1 and _prefer_materialised(a.shape[0], v.shape[1]): + return np.asarray(0.5 * (a + a.T) @ v) return np.asarray(0.5 * (a @ v + a.T @ v)) diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index e0d64a8..7bd0b62 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -1130,6 +1130,13 @@ def test_a_genuine_shortfall_is_still_reported(self): solve_modes(gk, gm, n_modes=2) +def _prefer_materialised_for(n_rows: int, n_cols: int) -> bool: + """Thin wrapper so the width-rule cases read as arithmetic.""" + from pybmodes.fem.solver import _prefer_materialised + + return _prefer_materialised(n_rows, n_cols) + + class TestTheGuardDoesNotTaxEverySolve: """Residuals are measured on every eligible solve, healthy ones included, so measuring them must not cost more than it saves. @@ -1151,6 +1158,21 @@ def _pair(self, n): gk = b @ b.T return 0.5 * (gk + gk.T), 0.5 * (gm + gm.T) + def test_the_two_mechanisms_are_not_interchangeable(self): + """Blocking and the width test answer different questions, and + neither subsumes the other. Capping the block to the crossover + instead of testing the width was measured and is worse at both + ends: at ngd = 100 it splits one BLAS call into two for 206 us + against 116 us, and at ngd = 1000 it lifts the peak from 5.1 MB + to 21.3 MB.""" + from pybmodes.fem import solver + + # Large ngd: blocking is what bounds the peak, and the block it + # produces is already narrow, so the width test does not fire. + assert not _prefer_materialised_for(2000, solver._RESIDUAL_BLOCK) + # Small ngd: blocking barely bites, and the width test does. + assert _prefer_materialised_for(100, solver._RESIDUAL_BLOCK) + def test_symmetrised_products_match_materialising(self): from pybmodes.fem.solver import _modal_residuals @@ -1211,21 +1233,19 @@ def spy(k, m, vals, vecs, *, symmetrise=False): assert all(is_gk and is_gm and sym for is_gk, is_gm, sym in seen), seen @pytest.mark.parametrize("k_over_n", [0.01, 0.2, 0.5, 1.0]) - def test_no_pass_of_the_sweep_ever_sees_a_wide_block( + def test_no_pass_of_the_sweep_exceeds_the_block( self, k_over_n, monkeypatch, ): - """The invariant the memory bound rests on. + """The invariant the memory bound rests on: the peak follows the + block width, not the caller's ``n_modes`` — which defaults to + ``None`` and asks for the whole spectrum. Measured at ngd = 2000 + with every mode requested, the unblocked sweep peaked at 128 MB + against 10.3 MB blocked. - ``0.5 (A v + A.T v)`` is only cheaper than building ``sym(A)`` - while ``v`` is narrow, and nothing about the caller's ``n_modes`` - guarantees that — ``None`` is the public default and asks for the - whole spectrum. Blocking is what supplies the guarantee, so this - asserts the width every pass actually sees rather than a measured - peak. - - Measured, not asserted, and recorded here because it is why this - exists: at ngd = 2000 with every mode requested the unblocked - sweep peaked at 128 MB against 10.3 MB blocked. + This bounds the peak. It does *not* make the block narrow + relative to ``ngd``, which is a separate question settled by the + width test inside ``_apply`` — see + :meth:`test_a_small_system_still_gets_the_width_test`. Deliberately not a ``tracemalloc`` comparison. That reports what the allocator did, and a freed block from an earlier call gets @@ -1302,6 +1322,40 @@ def test_blocking_does_not_change_the_answer(self, k): rtol=1.0e-12, atol=1.0e-15, ) + @pytest.mark.parametrize("n", [30, 60, 100, 150]) + def test_a_small_system_still_gets_the_width_test(self, n, monkeypatch): + """Blocking bounds the block at 128 columns, which is narrow + against a large ``ngd`` but not a small one. Below ``ngd = 192`` + a full-spectrum request hands ``_apply`` a block wider than the + ``2n/3`` crossover, so the width test inside it is still load + bearing — dropping it cost 178 us against 116 us at ngd = 100, + for the same peak. + """ + from pybmodes.fem import solver + + assert n < 192, "the point of this case is that 128 is not narrow" + k = n # the n_modes=None default + assert min(solver._RESIDUAL_BLOCK, k) >= 2 * n / 3, ( + "this ngd no longer produces a wide block; pick a smaller one" + ) + assert _prefer_materialised_for(n, min(solver._RESIDUAL_BLOCK, k)) + + @pytest.mark.parametrize("n_rows,n_cols,built", [ + (400, 400, True), # full spectrum, the wide case + (2000, 2000, True), + (1500, 6, False), # a few modes out of many DOFs + (400, 4, False), + (100, 128, True), # small ngd: the block is not narrow + (400, 128, False), # large ngd: the same block is + ]) + def test_the_route_is_the_cheaper_of_the_two(self, n_rows, n_cols, built): + """Asserted against the cost model — ``3nk`` split against + ``2n^2`` built — rather than a measured peak, for the + platform-independence reason given above.""" + assert _prefer_materialised_for(n_rows, n_cols) is built + split, build = 3 * n_rows * n_cols, 2 * n_rows * n_rows + assert _prefer_materialised_for(n_rows, n_cols) == (build <= split) + @pytest.mark.parametrize("block", [1, 3, 16, 64, 1000]) def test_the_block_width_is_a_free_parameter(self, block, monkeypatch): """It may be tuned for the peak-versus-overhead trade without