Skip to content

fix: catch the silent symmetric-eigensolver breakdown on a near-singular mass matrix - #140

Merged
SMI-Lab-Inha merged 28 commits into
masterfrom
fix/ill-conditioned-mass-solver
Aug 13, 2026
Merged

fix: catch the silent symmetric-eigensolver breakdown on a near-singular mass matrix#140
SMI-Lab-Inha merged 28 commits into
masterfrom
fix/ill-conditioned-mass-solver

Conversation

@SMI-Lab-Inha

@SMI-Lab-Inha SMI-Lab-Inha commented Aug 12, 2026

Copy link
Copy Markdown
Owner

Fixes a silent wrong answer in the dense symmetric eigensolver, found while working on #35.

The defect

scipy.linalg.eigh reduces K x = λ M x through a Cholesky factor of the mass matrix. When a very light beam carries a very heavy lump that matrix goes nearly singular, the reduction loses accuracy, and LAPACK does not raise — it returns wrong frequencies. On a 100 m cantilever at a 4000:1 lump-to-beam ratio the reported fundamental was 0.103 Hz against a true 0.0436 Hz, and the answer wandered non-monotonically with mesh density (0.086 → 0.084 → 0.069 → 0.113 → 0.0857). It predates this branch: it reproduces on master with plain tip_mass.

solve_modes now measures the backward error of a dense symmetric solve and, when it is large, solves again through the general path and compares. It warns when it swaps, and SolverDiagnostics gains residual_fallback.

Why the acceptance rule looks the way it does

Every clause is there because a simpler version was demonstrably wrong. This PR went through fifteen review rounds; the rule below is what survived.

Rule The case that forced it
Dense path only eigsh(sigma=0) factorises K, so it does not have this failure; and its which="LM" window is a different mode set that must never be index-compared
Per mode, not on maxima A rigid-body mode floors the candidate's maximum and hides a corrupted elastic mode beside it
Judged by the ratio of the win No absolute bar separates a rescue from rigid roundoff — 0.076, 0.79 and 12.4 all measured on healthy models, the first below the failure threshold
Plus a resolved-error bound A wildly broken 1e6 against a candidate at 100 clears any ratio while both are garbage
Non-regressive Accepting replaces the whole spectrum, so rescuing one mode while ruining another is a trade
Same matrices throughout Both symmetric paths symmetrise internally; judging against the raw pair fails an exact solve and lets eig "win" by answering a different question
Same spectrum throughout Any sign or complex filter returns a different set of the same length, backfilled from higher up, so equal indices stop meaning equal modes
Always able to decline If eig raises, or the system is too large, or the ordering cannot be verified, the symmetric result stands

Three attempts to identify rigid-body modes and exclude them were abandoned: by eigenvalue scale (a rigid-only subset becomes its own reference), by strain (a genuinely soft mode has little of it too — this one silently disabled the guard and returned the 538 %-wrong answer again), and by which side of the failure threshold the value falls on.

What separates the two populations is measured rather than argued: genuine rescues improve by 1e5 to 1e10, rigid roundoff by 11× to 16×. Those numbers are a test, so the constants cannot drift from their evidence.

No existing result changes

The decisive-improvement condition is what guarantees it. On the bundled NREL 5MW land deck — whose adapter leaves M 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, which the FA/SS classifier depends on. The full suite passes with no test edited except one whose workaround for this very defect is no longer needed.

What this does not promise

Reliable and platform-independent for the case it was built for: a near-singular mass matrix with no rigid-body modes. Elsewhere safe but not always effective — where rigid modes and a near-singular mass coincide, QZ may return the theoretically real zero modes as complex pairs whose position in the spectrum differs between LAPACK builds, and the retry then declines rather than guessing. Stated in the module docstring, the CHANGELOG and the tests rather than left to be discovered.

Verification

  • 1216 passed, 3 skipped; ruff, mypy and the strict docs build clean.
  • Two new VALIDATION.md rows: closed-form recovery at a 4e5:1 mass ratio, and the rigid-body non-regression.
  • 600 healthy free-free pencils of varying size and rank deficiency swept: no spurious swaps.
  • tests/fem/test_ill_conditioned_mass.py pins each rule and the case that forced it, including that the unguarded symmetric result is more than 25 % out — so the guard cannot be quietly removed.

Also fixes a docs/RELEASE_CHECKLIST.md reference that moved to the Sphinx tree, and two tests that asserted LAPACK behaviour rather than pyBmodes logic and failed on the Linux BLAS.

…lar 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.

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 0d29001144

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
Comment on lines +399 to +402
scale = float(np.max(np.abs(eigvals)))
if scale <= 0.0:
return 0.0
elastic = np.abs(eigvals) > _RIGID_BODY_EIGVAL_RTOL * scale

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P1 Badge Detect rigid-only subsets without using their own maximum

When n_modes selects only rigid-body modes, their small numerical eigenvalues become the reference scale, so most or all of them satisfy this relative cutoff and are treated as elastic. For example, a symmetric rank-6 stiffness matrix with six free rigid DOFs and n_modes=3 can produce residuals near one, trigger the retry, and switch to _solve_dense_general, which filters non-positive zero modes and returns flexible modes in their place. The rigid-mode test added here avoids this because its requested subset also contains elastic eigenvalues; the cutoff needs an absolute/problem-scale test or must defer the verdict when the returned subset supplies no elastic reference.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Fixed, and this one was worth the trouble — the finding was right and chasing it properly turned up a second bug in my own fix.

You were correct about the mechanism: a subset containing only rigid-body modes takes its own numerical noise as the reference scale, re-admits them as elastic, and reads their ~1 residuals as a breakdown. I reproduced it — n_modes of 1, 3 and 6 on a six-free-DOF model gave verdicts of 2.06, 2.50 and 1.76 against a 0.1 threshold.

My first fix was your suggested alternative: classify by strain against ||K|| rather than by eigenvalue, so the verdict is independent of the requested subset. That is subset-independent, but it is wrong for a different reason. A genuinely soft mode on a stiff structure carries little strain too — the 0.084 Hz lump mode of a 1e10 N.m² beam came in at 4e-12 of ||K||, below the cutoff, so the guard stopped firing on the exact case it was written for and the answer silently went back to being 538 % wrong. The existing point-mass test caught it.

So I dropped the classification entirely. Nothing here separates a rigid-body mode from a soft elastic one reliably, and both attempts to try made things worse.

What made a false positive dangerous in the first place was not the detection but the fallback: it filtered out zero eigenvalues, so a spurious retry changed the mode set rather than just its precision. Fixing that inverts the problem. The retry now runs with keep_rigid_body=True, and since both candidates are measured identically and the swap needs a decisive win, a mode the metric cannot speak to says the same nothing twice and cannot tip the decision. A false positive is now merely a wasted solve.

The asymmetric production path keeps the BModes-matching > 0 filter it is validated against; only the retry opts in.

Verified: rigid-only subsets return the requested mode count at n = 1, 3, 6 and 10; the point-mass case is exact again (0.00 %); the mesh sweep from 7 to 201 elements converges instead of wandering. 1169 passed, ruff and mypy clean.

@codex review

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 7873c56bff

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
gk, gm, n_modes, keep_rigid_body=True,
)
_normalize_columns_l2(alt_vecs)
alt_worst = _max_residual(gk, gm, alt_vals, alt_vecs)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P1 Badge Keep rigid residuals from masking an elastic breakdown

When a near-singular model contains both rigid and elastic modes, taking the maximum across all modes prevents the decisive-win check from detecting corrupted elastic modes: the general candidate still has rigid-mode residuals around 1, so alt_worst cannot be less than 0.1 * worst unless the symmetric residual exceeds roughly 10. Fresh evidence from combining six zero DOFs with the ill-conditioned cantilever and applying an orthogonal basis rotation produced a symmetric fundamental that was 18–44% wrong with elastic residuals of 0.77–0.82, while alt_worst remained about 1 and residual_fallback stayed false. Thus a free-free structure can still silently return the failure this guard is intended to correct; the improvement comparison must avoid allowing unchanged rigid-mode noise to mask a decisive improvement in the elastic modes.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Fixed, and this one directly invalidates the argument I made in the previous round, so thank you for not letting it stand.

I claimed the meaningless rigid-body component would "cancel because both candidates are measured the same way". That is true of a per-mode comparison and false of a comparison on maxima, which is what I actually wrote. A rigid mode reads ~1 in both, so it floors max(alt_r) near 1, and no improvement among the elastic modes can clear 0.1 * max(sym_r) unless the symmetric solve is worse than ~10. Your reconstruction of the consequence is exactly right.

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 — the retry already preserves rigid-body modes — so equal indices mean the same mode. I also refuse an alternative that recovered fewer modes than the symmetric solve outright, since losing a mode is never an improvement whatever the residuals say.

The regression test builds your case directly: six free rigid DOFs bolted onto the ill-conditioned cantilever, the pair rotated by an orthogonal basis so the blocks are not separable by inspection. It pins the mechanism as well as the outcome, asserting that the alternative maximum really is floored above 0.1, that a maxima comparison really does see no decisive win, and that a per-mode one really does — so if anyone later simplifies this back to comparing maxima, that test fails rather than the behaviour silently regressing.

1172 passed, ruff and mypy clean. The three end-to-end behaviours still hold: the point-mass case exact at 0.00 %, rigid-only subsets returning their requested mode count at every size, and the mesh sweep converging.

@codex review

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 8eda44ec4a

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
Comment on lines +420 to +421
if n_alt < n_sym:
return np.zeros(0, dtype=bool)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Preserve negative modes before comparing by index

When geometric softening makes K indefinite and the caller requests only a subset, _solve_dense_general(..., keep_rigid_body=True) still filters genuinely negative eigenvalues below its rounding floor, then fills the subset with later positive modes. The alternative can therefore have the same length as the symmetric result while describing a shifted spectrum; this check accepts it and the per-index residual test can replace the result, silently dropping an unstable mode and shifting every subsequent mode. Verify eigenvalue correspondence or retain negative modes on the retry rather than treating equal counts as proof that indices match.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Fixed, and you were right that equal counts prove nothing.

Preserving only rounding-level zeros still filtered genuinely negative eigenvalues, which is not hypothetical here: run(gravity=...) past a column's buckling weight makes K indefinite, and that path landed in this same release. The alternative then came back the same length, backfilled from higher up, describing a shifted spectrum — so the per-index comparison was reading two different spectra against each other.

Rather than verify correspondence after the fact I removed the mismatch at the source. eigh filters nothing, so the retry now filters nothing either: every real finite eigenvalue, sign included. Both paths then return the n_modes smallest real eigenvalues and equal indices describe the same mode by construction rather than by assumption. The flag is renamed preserve_full_spectrum to say what it actually does; the asymmetric production path keeps its positive-only filter and is untouched.

Regression test builds the indefinite case directly and asserts the trap explicitly — the default filter and the preserving one return the same length while the preserving one alone contains the unstable mode, and eigh agrees the mode exists.

Separately, CI caught two of my own tests being wrong rather than the code: they asserted LAPACK behaviour and passed on Windows while failing on the Linux BLAS. The degenerate-pair test demanded 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 actually matters. And the masking mechanism from the last round is now stated as arithmetic on residual vectors rather than coaxed out of a matrix pair, so it pins the two decision rules exactly and cannot drift with the vendor library.

1176 passed locally, ruff and mypy clean.

@codex review

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: fbf9e06a84

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
# stays exact there.
residual_fallback = False
if sym:
sym_r = _modal_residuals(gk, gm, eigvals, eigvecs)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P1 Badge Measure residuals against the symmetrized matrices

For matrices with asymmetry below symmetry_rtol, the symmetric paths solve 0.5 * (K + K.T) and 0.5 * (M + M.T), but this evaluates their residuals against the original matrices. In a high-dynamic-range model, an accepted skew of 1e-12 * max|K| can be large relative to a soft mode, so an accurate eigh result appears to have residual greater than 0.1; the retry then decisively wins only because eig solves the different, unsymmetrized problem and replaces the spectrum. For example, an accepted 3×3 stiffness with low eigenvalues around 1e-12 and skew 8e-13 shifts the first eigenvalue by about 8–10%. Compute this check against the same symmetrized matrices used by the selected symmetric solver.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Fixed. This one was the subtlest yet and the reasoning is exactly right.

The symmetric paths symmetrise internally, and I was judging their output against the raw matrices — so the guard was asking a question the solver had never been asked. The skew _is_effectively_symmetric tolerates is only small relative to max|K|, so in a wide-dynamic-range model it can be the same size as a soft mode's own eigenvalue. An exact symmetric solve then reads as broken, and eig on those same raw matrices wins decisively purely by answering the different question, handing back the skewed spectrum in place of the symmetric one the caller was promised.

Both the measurement and the retry now use the symmetrised pair, so the whole comparison is about one problem.

The regression test builds your construction: a 3×3 with eigenvalues from 1 down to 1e-12 and skew just inside the tolerance. Two details I had to get right and mention in case they matter to the next reader — the symmetry test compares max|A - A.Tᵀ|, which is twice the off-diagonal skew, so the skew has to stay under half the tolerance to be accepted; and the skew has to couple the two softest modes, since between the stiff ones it is negligible against their own scale and demonstrates nothing. It asserts the pair is accepted as symmetric, that no retry fires, that the returned spectrum is the symmetric one, and that the very same modes measured against the raw matrices would have read above the threshold.

1179 passed, ruff and mypy clean.

@codex review

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: b724726aeb

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
Comment on lines +307 to +309
gk_s = 0.5 * (gk + gk.T)
gm_s = 0.5 * (gm + gm.T)
sym_r = _modal_residuals(gk_s, gm_s, eigvals, eigvecs)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Measure diagnostics against the symmetrized pair

When accepted skew is present, these symmetrized matrices correctly drive the retry decision, but _build_diagnostics is later passed the raw gk/gm and recomputes residuals against a problem the returned modes did not solve. With the added 3×3 accepted-skew construction, the correct symmetric result reports residual_fallback=False while diag.max_residual is about 0.37, falsely flagging a healthy solve as defective in the certification telemetry. Reuse the residuals from the solved pair or pass the symmetrized matrices when building diagnostics.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Fixed, and this one is squarely on me rather than an oversight — I considered exactly this last round and decided against it.

My reasoning was "don't perturb a documented field", so I moved the retry decision onto the symmetrised matrices and deliberately left _build_diagnostics on the raw ones. That traded a correctness property for a compatibility one, and it was the wrong way round: the whole point of max_residual is to be auditable telemetry, and reporting ~0.37 for a solve that is exact defeats it more thoroughly than a changed number would have.

solve_modes now names the pair the returned modes actually solve once — symmetrised on the symmetric paths, raw on the general one — and uses it for the retry decision and the diagnostics alike. The field's docstring now states which basis it is on, since a backward error means nothing without naming the problem it is measured against.

Genuinely symmetric decks are unaffected, since with no skew the symmetrised pair is the raw pair, which is why nothing in the existing suite moved.

1180 passed, ruff and mypy clean, strict docs build clean.

Six rounds, six real defects, every one in code I had convinced myself was right — including this one, which I had actively reasoned about and got backwards. Thank you for the persistence; this guard is materially better for it.

@codex review

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 55e30ef985

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
& np.isfinite(eigvals_real.real)
& (eigvals_real.real > 0.0)
)
valid = np.isreal(eigvals_real) & np.isfinite(eigvals_real.real)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P1 Badge Preserve numerically complex real-spectrum modes on retry

When retrying a highly ill-conditioned symmetric free-free problem, QZ can represent theoretically real zero modes as small complex-conjugate pairs that np.real_if_close(..., tol=1000) does not coerce, so this mask still deletes them even with preserve_full_spectrum=True. Because enough later real eigenvalues remain to fill n_modes, the equal-count guard passes and elastic modes backfill the removed indices before the fallback is accepted. Fresh evidence from the current TestRigidModesCannotMaskAnElasticBreakdown._rigid_plus_ill_conditioned fixture is that requesting 10 modes returns only 2 of its 6 rigid modes on this build, silently shifting the spectrum despite the new preservation path.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Confirmed and fixed — and we arrived at this one from opposite directions simultaneously. I was auditing my own retry for exactly this hole when your comment landed; your evidence is the more alarming half.

I had reproduced the mechanism but drawn the wrong conclusion from it. On this machine eig drops six of 249 eigenvalues on a healthy fixture, all rounding-induced conjugate pairs at the stiff end, far above anything a caller asks for. That made my first instinct — demand every mode survive — far too strict to be usable. Your build drops the zero modes, which is the case that actually matters, and shows why the count guard was never sufficient.

The invariant is narrower than either: nothing discarded may fall inside the returned window. The retry now checks that directly and declines to swap when it cannot verify it. Declining rather than reaching for something cleverer is deliberate, and your evidence is the argument for it — backfilling a dropped zero mode with an elastic one would be precisely the silent spectrum shift this guard exists to prevent, and a guard added to stop a silent wrong answer must not be able to introduce one.

That leaves an honest scope statement rather than a claim of generality, now written into the module docstring, the CHANGELOG, and the tests:

  • Reliable and platform-independent for the case the guard was built for — a near-singular mass matrix with no rigid-body modes. That is the 2.4×-wrong case, and it behaves identically everywhere.
  • Safe but not always effective where rigid-body modes and a near-singular mass matrix coincide. The result is then no worse than without the guard, and max_residual still reports the problem.

The rigid+ill-conditioned test now asserts the portable property — the requested modes come back, ascending, none lost — instead of a rescue that is genuinely build-dependent. The rescue itself stays pinned on the unambiguous case.

The retry also gets its own entry point instead of a flag on _solve_dense_general, so the asymmetric production path is provably untouched rather than untouched-by-inspection.

1183 passed, ruff and mypy clean.

@codex review

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.

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 0fc9e0cb94

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
Comment on lines +353 to +355
_decisively_improved_modes(
sym_r, alt_r, alt_vals.size, eigvals.size,
)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P1 Badge Match the sparse spectrum before comparing residuals

When a symmetric system takes the sparse path and has negative modes, eigsh(..., sigma=0, which="LM") returns modes nearest zero, whereas _general_spectrum_for_retry sorts algebraically and returns the most negative modes first. If the sparse residual exceeds the threshold—such as for a large, near-singular, post-buckling model—the per-index comparison therefore pairs different modes, and an exact far-negative alternative can trigger a swap even when the corresponding sparse modes did not improve. Fresh evidence beyond the earlier sign-filter issue is that the current retry also runs after sparse_shift_invert, whose selection window differs from the retry's window; verify eigenvalue correspondence or make both paths select the same spectrum before using this mask.

Useful? React with 👍 / 👎.

Comment thread src/pybmodes/fem/solver.py Outdated
Comment on lines +347 to +349
alt_vals, alt_vecs, ordering_sound = _general_spectrum_for_retry(
gk_s, gm_s, n_modes,
)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Preserve the symmetric result when the retry fails

When the symmetric solve returns a residual above the threshold but scipy.linalg.eig fails to converge on the same ill-conditioned pencil, this unguarded retry raises LinAlgError and discards the already available symmetric result. That turns the guard from a safe, best-effort correction into a new hard failure precisely on the defective inputs most likely to make the alternative solver fail; catch retry failures and decline the swap while retaining the original result and diagnostics.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Both fixed. The first one had a root-cause fix that removes the class rather than handling it.

Sparse spectrum mismatch (P1). Correct, and the diagnosis is exact: which="LM" on OP = K⁻¹M selects the modes nearest zero in magnitude while the retry selects the algebraically smallest, so with negative eigenvalues present those are different sets and the per-index comparison pairs unrelated modes.

Rather than reconcile the two selections, I removed the case. The retry is now scoped to the dense symmetric path only — 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 failure this guard exists for; eigsh(sigma=0, mode='normal') factorises K and is therefore unaffected by a near-singular M. The mesh sweep that motivated this whole change already demonstrated that: it returned correct frequencies on exactly the meshes large enough to take the sparse path, and wrong ones below. So nothing is lost by excluding it, and the window mismatch cannot arise. There is now a test asserting both halves — that the sparse path is not retried even with a terrible residual, and that it gets the ill-conditioned case right unaided.

Unguarded retry (P2). Also correct, and the framing is the right one: a guard added to prevent a silent wrong answer must not be able to convert a usable result into a hard failure, least of all on the defective inputs most likely to break the second solver too. The retry now catches LinAlgError / ValueError, declines the swap, keeps the symmetric result and diagnostics, and logs the reason. Tested for both exception types.

1187 passed, ruff and mypy clean.

@codex review

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: fde9191976

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

stacklevel=2,
)
eigvals, eigvecs = alt_vals, alt_vecs
path = "dense_general"

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Avoid reporting an overlarge request as filtered modes

When a dense symmetric problem accepts the residual retry and n_modes exceeds the matrix dimension, setting path to dense_general makes the later mode-count warning claim that the shortfall came from complex or non-positive eigenvalues. The retry helper actually preserves non-positive eigenvalues, and the shortfall is simply the benign min(n_modes, ngd) truncation that the surrounding warning logic intends to exclude; for example, a 117-DOF ill-conditioned system requested with 1000 modes returns all 117 but reports the pencil as defective. Gate that warning on modes actually discarded by the general solver rather than the post-retry path label.

Useful? React with 👍 / 👎.

Comment thread CHANGELOG.md Outdated
Comment on lines +23 to +25
`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

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Document the dense-only retry scope

The release note says every symmetric solve is checked and redone through the general path, but the implementation now deliberately restricts retries to path == "dense_symmetric"; sparse results are never retried even when their residual exceeds the threshold. This also conflicts with the preceding claim that eigsh factorizes the mass matrix, whereas the implemented rationale says shift-invert factorizes K. Update the public change description so users do not expect residual_fallback or a corrective retry on sparse solves.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Both fixed, plus one I found in a self-review pass before it was reported.

Mode-count warning (P2). Correct. After a swap the path label reads dense_general, so the shortfall warning claimed modes had been filtered as complex or non-positive — when the retry preserves the whole spectrum and nothing was filtered at all. Your 117-DOF / 1000-mode example is exactly right. The warning now keys on modes actually discarded rather than on the label, comparing against min(n_modes, ngd). That also fixes the same spurious report on the pre-existing general path, which had the identical problem for an overlarge request.

CHANGELOG (P2). Also correct, and a good catch on a doc claim I had invalidated myself: my own scoping change made "every symmetric solve" and "eigsh factorises the mass matrix" both false, and I updated the code without updating the note. Rewritten to say the retry is dense-only, that sparse results never set residual_fallback, and why — eigsh(sigma=0, mode='normal') factorises K, so it is unaffected by a near-singular M, and its which="LM" window is a different mode set anyway.

Found in self-review: the retry was unbounded in size. A sparse solve that fails to converge falls back to dense at any size, so the retry was reachable on a 5000-DOF system where a dense eig runs for minutes on a result already in hand. By the same principle as your earlier P2 about LinAlgError — a guard against a silent wrong answer must not be able to introduce a new failure — a silent multi-minute hang is a worse outcome than the imprecision. It now declines above residual_retry_max_ndof (default 2000) and leaves the backward error to the diagnostics.

I also ran an adversarial sweep over the guard's edge cases: n_modes > ngd, n_modes=None, a 1×1 system, all-rigid K, indefinite K with a 10^15 dynamic range, a structurally singular M, an all-complex retry, and NaN input. Everything behaves; the only two that raise are eigh itself on a non-positive-definite M and on NaN, both pre-existing, both loud rather than silent, and both gated upstream by check_model.

1192 passed, ruff and mypy clean.

@codex review

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 603eae8b50

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
if ordering_sound
else np.zeros(0, dtype=bool)
)
if improved.any():

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P1 Badge Reject retries that degrade other requested modes

When QZ decisively improves one mode but makes another requested mode less accurate, improved.any() accepts the retry and lines 419–421 replace the entire spectrum. For example, on a 4×4 SPD pencil with cond(M) ≈ 1e16, the symmetric residuals can be [0.52, 1.1e-3, 5.7e-6, 0.017] while the general residuals are [1.1e-10, 3.0e-11, 3.5e-5, 0.26]; this condition swaps because mode 0 improved, thereby replacing a previously acceptable fourth mode with one above the failure threshold. Require the whole candidate spectrum to be non-regressive, rather than accepting it when any single mode improves.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Fixed, and this is the sharpest finding of the series — it catches a gap in the shape of the rule rather than in its plumbing.

The rule 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. Your residual vectors are now the test.

The comparison returns both verdicts and the caller requires improvement somewhere and regression nowhere.

One design note on how I defined regression, since the obvious version is wrong. A bare alt > sym test would have blocked most legitimate rescues: rigid-body modes read ~1 in both candidates and wobble either way, so 1.0001 against 1.0 would count as a regression and veto any rescue sitting beside a free-free mode. Regression therefore mirrors improvement exactly — above the threshold and worse there by the same factor that would have counted as decisive in the other direction. That reuses one notion of "acceptable" instead of inventing a second, and it keeps the two cases apart that have to stay apart:

  • your mode 3, 0.017 → 0.26: above the threshold and 15× worse. Regression.
  • your mode 2, 5.7e-6 → 3.5e-5: six times worse and nowhere near the threshold. Not a regression, and refusing on it would be over-constraining.

Both are pinned, along with the rigid-body noise case and an end-to-end test that corrupts one mode of the candidate and asserts the swap is declined.

1196 passed, ruff and mypy clean.

@codex review

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: dc26e04879

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
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)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P1 Badge Reject candidates that cross the residual threshold

When a previously acceptable mode crosses the threshold by less than the configured tenfold factor, this condition does not mark it as regressed, so an improvement elsewhere still replaces the whole spectrum with a newly failing mode. For example, sym_r=[0.52, 0.09] and alt_r=[1e-10, 0.8] produce an improvement and no regression, even though the second mode moves from below the 0.1 threshold to a residual of 0.8; this also does not involve the rigid-body-noise case because those modes already start above the threshold. Treat an acceptable-to-failing threshold crossing as a regression in addition to the decisive-worsening comparison.

Useful? React with 👍 / 👎.

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 430ec073dc

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
# 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)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Avoid rematerializing matrices after ordinary sparse solves

When a large symmetric problem successfully takes sparse_shift_invert with return_diagnostics=False, this unconditionally constructs two full dense symmetrized matrices even though the retry is excluded by path == "dense_symmetric" and diagnostics are never built. _solve_sparse_shift_invert has already symmetrized the same pair internally, so every ordinary sparse solve now pays for a second pair of O(ngd²) allocations and matrix additions, adding substantial memory traffic on exactly the large systems for which the sparse path exists. Construct res_k/res_m lazily only when the dense retry or returned diagnostics needs them.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

SMI-Lab-Inha commented Aug 12, 2026

Copy link
Copy Markdown
Owner Author

Fixed, and this one I had already found and misjudged — which makes it the more useful of the two outcomes.

During my hardening pass I noted exactly this ("two O(n²) matrix adds on every solve… acceptable but wasteful when the retry never runs") and moved on. It is not acceptable: at ngd = 1500 the eager symmetrisation peaks at 54 MB against 0.3 MB, a factor of 187, and it falls on precisely the large sparse solves the sparse path exists to keep cheap. Worse, I introduced it — before this PR the diagnostics measured against the raw matrices and allocated nothing.

The identity 0.5 * (A + A.T) @ v == 0.5 * (A @ v + A.T @ v) means the basis change can be carried as a flag and applied through two thin products instead of building the pair. Three tests: the two forms agree to 1e-12, the flag is not inert (an asymmetric pair measures differently with and without it), and the peak allocation stays under a quarter of a single matrix.

The symmetrised pair is still materialised in one place — the retry itself, because scipy.linalg.eig needs matrices rather than products. That branch is rare and already size-capped.

Worth noting what this says about my hardening pass. I ran it specifically to pre-empt findings, and it did surface this one — I just weighed it wrong, because I was thinking about the retry branch where the cost is bounded and rare, rather than the diagnostics path where it is unconditional. Noticing a cost and mis-scoping who pays it is a distinct failure from not noticing it, and measuring rather than estimating would have settled it in a minute.

(Reposted: my previous comment lost two code spans to shell substitution.)

1219 passed, ruff and mypy clean.

@codex review

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: acb57844de

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread src/pybmodes/fem/solver.py Outdated
Comment on lines +367 to +368
gk_s = 0.5 * (gk + gk.T)
gm_s = 0.5 * (gm + gm.T)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Delay symmetrization until the retry actually runs

For every eligible dense symmetric solve—including healthy solves and large sparse-to-dense fallbacks up to 2000 DOFs—these lines still allocate both full symmetrized matrices before the residual is compared with the retry threshold. Fresh evidence beyond the earlier sparse-path finding is that the new thin-product helper is not used for this initial check, so a retry that never runs still pays the repeated O(ngd²) memory traffic; compute sym_r from the raw pair with symmetrise=True, then materialize gk_s and gm_s only inside the threshold-triggered branch where eig needs them.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Fixed — and this is 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. So every healthy dense solve up to the size cap paid two dense allocations to discover it did not need them. The pre-threshold check now goes through the products too; the pair is built inside the branch, where scipy.linalg.eig genuinely needs matrices.

Two things worth recording, since both are about how I verified rather than what I changed.

The test I first wrote for this was wrong, and would have been bad even if it passed. I asserted the whole solve allocates less than half a matrix. That is impossible — 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's arrays rather than a copy, which is the property actually at stake. A memory-threshold test would have been measuring the solver, not the guard.

My previous fix was partial in a way I should have caught. I applied the helper at one call site and did not check the others. Grepping for the pattern I had just introduced would have found this in seconds, and I did not do it.

1220 passed, ruff and mypy clean.

@codex review

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: dfbb96c99d

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

# 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)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Avoid full-width products for full-spectrum solves

When n_modes is omitted—the public default—or requests most of the spectrum, eigvecs is ngd × ngd, so these are not thin products: _apply evaluates both A @ V and A.T @ V, creating several full matrix-sized temporaries for each of K and M. Thus every healthy dense solve within the 2000-DOF retry cap still incurs substantial peak memory and four extra dense matrix multiplications before the threshold is checked; locally, the full-width residual call peaked at about five times one matrix allocation, whereas the added test covers only five vectors. Handle the full-width case separately or reuse the symmetrized matrices/products already required by the dense solve.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Fixed in e93347f — and the measurement corrected my first attempt, so it's worth recording what it said.

You're right that the products stop being thin. The two routes peak at 3nk and 2n² bytes of temporaries, so the split product wins only while the block is narrow, and n_modes=None (the public default) returns the whole spectrum. _apply now selects by width.

What I got wrong on the first pass was where to put the switch. A flop-count estimate said k = n/3; measuring the allocator at n = 600 said otherwise.

     k    k/n     products  materialise
   300   0.50       4.32 MB       5.76 MB    <- products still ahead
   399   0.67       5.75 MB       5.76 MB    <- crossover
   401   0.67       5.77 MB       5.76 MB
   600   1.00       8.64 MB       5.76 MB

The real crossover is 2n/3, so the estimate would have taken the dearer route across a third of the range — for exactly the narrow-ish blocks the retry path actually sees. The threshold is now 3k >= 2n, which is the arithmetic identity behind that table rather than a tuned constant.

Two tests cover it. One asserts _apply lands within 5 % of the cheaper of the two routes at k/n in {0.01, 0.2, 0.5, 0.9, 1.0}, which spans both sides of the switch and pins the full-width case you flagged. The other asserts the routes agree numerically, since the whole point is that they differ only in cost.

I did consider reusing the symmetrised matrices from the dense solve. The residual measurement runs before we know whether a retry is warranted, and on the healthy path there is no retry, so materialising there would put the cost back on every solve to serve the rare one. Choosing by width keeps the common case cheap without a second code path.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

CI caught the width test on Linux while it passed on Windows, so that test is gone in cad2de1 — and the failure is worth recording, because the route it flagged was correct.

It compared three live tracemalloc peaks: _apply, the split form, and the built form. 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 the order they are measured in. At k/n = 0.9 the chosen route reported 2.43 MB against 2.30 MB for the built form measured separately — but 2.43 MB is exactly the analytic built-form peak ( copy at 1.28 MB plus the n × k result at 1.15 MB). It had taken the right route; the later measurement simply got a reused buffer.

Same trap as the two LAPACK-dependent tests earlier in this PR: asserting measured platform behaviour rather than the thing under test.

So the decision is now a pure predicate, _prefer_materialised(n_rows, n_cols), asserted against the cost model — 3nk split, 2n² built — which is deterministic on every platform. The parametrisation gained k/n = 0.67 to sit on the crossover, and there are now explicit cases for the full-width block that prompted your finding and for the thin block it must not disturb. The numerical-agreement test straddles the crossover at k ∈ {1, 7, 26, 27, 40} for n = 40, where 2n/3 = 26.7.

The measurement that established the model stays in the comment above the predicate and in my previous reply, which is the right place for it — it justified the constant, it should not be re-run as an assertion.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

@codex review

@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

@codex review

Four commits since the last review, all responding to the full-width residual finding:

  • e93347f picks the sym(A) v route by block width. The crossover is measured, not derived — 3nk split against 2n² built, crossing at k = 2n/3, where a flop-count estimate had said n/3.
  • cad2de1 replaces the test for that with a deterministic one. The original compared three live tracemalloc peaks, which measures allocator reuse rather than cost; it passed on Windows and failed CI on Linux at a width where the route taken was correct. The rule is now a pure predicate asserted against the cost model.
  • a34fd66 and 96040ec are documentation. The second corrects an overstatement: the regression guard refuses on alt > sym strictly, so an accepted retry proves an acceptable mode did not get worse, not that it improved.

CI is green on 96040ec, and I have also run the integration suite locally against the upstream decks, which the GitHub runner cannot do — 127 passed, including the six BModes cert cases with OC3 Hywind still at 0.0003 % against a 0.01 % tolerance.

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 96040ec551

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment on lines +706 to +708
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))

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Reduce the full-width residual peak, not only the product peak

When n_modes=None or most modes are requested, this routing reduces the number of full matrix multiplications but does not resolve the reported peak-memory problem. Fresh evidence in the new route is that its 3nk versus 2n² model stops at _apply: _modal_residuals retains the full-width kx while constructing mx, then retains both while allocating mx * eigvals and the subtraction result, so the call still simultaneously holds roughly four additional n×n arrays (over 100 MB at the 2000-DOF retry cap, excluding the eigenvectors). Compute these residuals in blocks or reuse/in-place the products so the public full-spectrum default no longer incurs that peak.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

Confirmed and fixed in 022b06c. The finding was exactly right, including the figure: I measured 128 MB at the 2000-DOF cap.

My previous fix bounded _apply but not its caller, which is the more important half. _modal_residuals held kx while building mx, then both while forming mx * eigvals and the difference — four full-width arrays live purely to measure.

The way out is that the residual is per mode, so no step needs every mode present at once. The sweep now runs in column blocks and accumulates the two norm vectors, which ties the peak to the block width instead of to an argument the caller chooses.

ngd = 2000, one dense copy = 32 MB

     k      blocked    unblocked   saving     t_blk   t_unblk  agree
     6       0.4 MB       0.4 MB     1.0x    0.003s    0.003s  True
    50       3.2 MB       3.2 MB     1.0x    0.007s    0.008s  True
   500      10.2 MB      32.0 MB     3.1x    0.063s    0.050s  True
  2000      10.3 MB     128.0 MB    12.5x    0.397s    0.268s  True

Two things worth drawing out. The common thin case is completely untouched — same peak, same time, one pass per matrix, computed exactly as before. And the full-width case costs 0.13 s more, because blocked products do 2n²k flops where a materialised sym(A) would do n²k; I took that deliberately, since it buys a 12.5× memory reduction on a solve where eigh at that size already dominates.

This also makes the width test I added last round unnecessary, so _prefer_materialised is gone. Blocking guarantees _apply never sees a wide block, which is a stronger guarantee than choosing a route by width and one fewer thing to keep true.

Tests follow the same principle as last round: the invariant is asserted by spying on the width every pass actually sees, not by comparing tracemalloc peaks. Coverage is the block boundary (k = 1, 7, 128, 129, 260 against a directly computed expectation), that no pass exceeds the block, that every mode is measured exactly once, and that a six-mode request is still one pass per matrix.

I also A/B'd the solver against 96040ec on the FEM suite to be sure the loop costs nothing in practice: 22.94 s before, 22.93 s after.

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.
@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

@codex review

@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

@codex review

Two commits since the last pass, both on the residual sweep.

022b06c addresses the full-width peak finding. The residual is per mode, so the sweep now runs in column blocks and accumulates the two norm vectors, which ties the peak to the block width rather than to n_modes. Measured at ngd = 2000 with every mode requested: 128 MB down to 10.3 MB, agreeing to 1e-12, with the thin case unchanged in both peak and time. This also removed _prefer_materialised — blocking guarantees _apply never sees a wide block, which is a stronger guarantee than picking a route by width.

4dd48df justifies the block width. 128 was a round number; measured across widths it sits at the knee, 1.4 s at 16 and 32 columns where per-call BLAS overhead dominates, 0.42 s at 128, against a 0.38 s floor that 256 and above buy with two to twelve times the peak. There is now also a test that the width changes nothing but the peak.

CI is green on 4dd48df. I have re-run the integration suite locally against the upstream decks since the solver change, which the runner cannot do: 127 passed, the six BModes cert cases included.

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 4dd48df99f

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

# 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

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Scale the residual block to small matrices

When ngd < 192 and the caller uses the public n_modes=None default—or otherwise requests most modes—the first block can exceed the 2n/3 crossover documented immediately above. Thus removing _prefer_materialised does not actually guarantee that _apply receives a narrow block: for a 100-DOF full-spectrum solve it receives a square block and performs two GEMMs per matrix instead of the cheaper materialized-symmetry route, regressing the common small dense case. Cap the block width relative to ngd or retain width-based routing.

Useful? React with 👍 / 👎.

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.
@SMI-Lab-Inha

SMI-Lab-Inha commented Aug 13, 2026

Copy link
Copy Markdown
Owner Author

Correct, and the over-claim was mine. Fixed in 31b4d64 by restoring the width test.

I said blocking made _prefer_materialised unnecessary. That holds only while 128 columns is narrow relative to ngd, and below ngd = 192 it is not — a full-spectrum request there hands _apply a block past the 2n/3 crossover, exactly as you describe. The two mechanisms answer different questions and neither subsumes the other: blocking bounds the peak as the mode count grows, the width test picks the cheaper route within whatever block it gets.

You offered two remedies and they are not equivalent, so I measured both before choosing.

full spectrum, symmetrise=True. 'auto' = width test inside _apply

  ngd |      products (now) |            cap 2n/3 |          auto route
   30 |    31.3 kB    21.1 us |    21.6 kB    33.5 us |    31.3 kB    17.0 us
  100 |   324.2 kB   177.8 us |   215.1 kB   205.8 us |   324.2 kB   116.1 us
  192 |   792.5 kB  1034.2 us |   792.5 kB   863.8 us |   792.5 kB   677.6 us
  400 |  2055.3 kB  4775.6 us |  3416.5 kB  4628.2 us |  2055.3 kB  4603.4 us
 1000 |  5137.0 kB    45.1 ms | 21339.7 kB    44.3 ms |  5137.0 kB    47.1 ms

Capping the block relative to ngd fixes the small case in the wrong direction: it saves about 100 kB but splits one BLAS call into two, so ngd = 100 goes from 178 us to 206 us — and at ngd = 1000 it lifts the peak from 5.1 MB to 21.3 MB, undoing most of what blocking just bought. Restoring the width test is better on both axes at small ngd (116 us for the same peak) and is inert at large ngd, where the block is already narrow and it never fires.

Tests cover the specific case you found — a sub-192 ngd at the n_modes=None default, asserted against the cost model rather than a measured peak — plus a case pinning that the two mechanisms are not interchangeable, so this claim cannot be made again silently. The route table now includes (100, 128) -> built and (400, 128) -> split, the same block width landing on either side of the rule.

Full suite 1251 passed, CI green, and I have re-run the integration suite locally against the upstream decks since the solver change: 127 passed, the six BModes cert cases included.

@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

@codex review

@SMI-Lab-Inha

Copy link
Copy Markdown
Owner Author

@codex review

31b4d64 restores the width test inside _apply, closing the small-ngd finding. You were right that removing it did not leave a guarantee: below ngd = 192 a full-spectrum request still hands _apply a block past the 2n/3 crossover, so blocking and the width test answer different questions and neither subsumes the other.

Of the two remedies you offered I took the second, having measured both. Capping the block relative to ngd saves about 100 kB at ngd = 100 but splits one BLAS call into two (206 us against 116 us), and at ngd = 1000 lifts the peak from 5.1 MB to 21.3 MB, undoing most of what blocking bought. Restoring the width test is better on both axes at small ngd and inert at large ngd, where the block is already narrow and it never fires.

Tests cover the case you found directly, plus one pinning that the two mechanisms are not interchangeable so the claim cannot be made again silently.

CI green on 31b4d64, 1251 in the default suite, and 127 integration tests re-run locally against the upstream decks since this commit.

@chatgpt-codex-connector

Copy link
Copy Markdown

Codex Review: Didn't find any major issues. Nice work!

Reviewed commit: 31b4d64b74

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

@SMI-Lab-Inha
SMI-Lab-Inha merged commit f40e87b into master Aug 13, 2026
6 checks passed
@SMI-Lab-Inha
SMI-Lab-Inha deleted the fix/ill-conditioned-mass-solver branch August 13, 2026 04:43
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant