Saloon: symmetric eigendecomposition (eig) via cyclic Jacobi#47
Merged
Conversation
Phased plan for eigenvalues/eigenvectors on Lagoon arrays, the marquee consumer of %cplx. Phase A: symmetric/Hermitian via the Jacobi algorithm (real eigenvalues, orthonormal eigenvectors, pure %i754, only needs sqrt). Phase B: general real -> complex via Hessenberg + Francis double-shift QR (%cplx). Decision: non-symmetric input crashes (assert symmetry) rather than silently symmetrizing. Tests check invariants (A*v=lambda*v, V'V=I, reconstruction) vs NumPy eigh, not raw eigenvectors (sign/order/degeneracy). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Phase A1 of the eig design: real symmetric eigenvalues/eigenvectors on Lagoon %i754 arrays (f32/f64), the marquee Saloon capability. Adds a %linalg section to the sa core: +eig -> [vals=ray vecs=ray] (eigenvalues 1-D, eigenvectors as columns) +eigvals -> ray (eigenvalues only) +eigvecs -> ray (eigenvectors only) plus building blocks (bloq-dispatched scalar float helpers, symmetric check, off-diagonal/Frobenius norms, diag, Givens row/column rotations, sweep-once). Asserts squareness and exact symmetry (crashes otherwise, per design) so a nonsymmetric matrix can't silently get the wrong algorithm. Cyclic Jacobi needs only sqrt, converges quadratically; 60-sweep cap with a warning rather than a silent non-converged result. Also fixes a latent bug surfaced by building Saloon against the current Lagoon: trans-scalar/fun-scalar used exhaustive `?- kind`, but Lagoon's kind union has grown (%unum, %cplx). Switched both to `?+ kind !!` so unsupported kinds crash gracefully instead of failing to compile. Tests (tests/lib/saloon-eig, generated by tools/eig_check.py against numpy.linalg.eigh) cover the four standard sources: closed-form (1-D Laplacian), torture matrices (Rosser 8x8, Wilkinson W7+, Pascal, Hilbert), real-world SPD (cycle Laplacian, Gram), and seeded random SPD. Per matrix: A*V = V*diag(vals) (order/sign-free) and sort(computed) = sort(eigh). All 18 pass in ~20s. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This was referenced Jun 7, 2026
sigilante
added a commit
that referenced
this pull request
Jun 8, 2026
Extends +eig to complex Hermitian matrices via complex cyclic Jacobi. +eig now dispatches on kind: %cplx -> Hermitian path, else the symmetric path. The unitary rotation J that zeros a_pq has real diagonal c and complex off-diagonal b = s*(a_pq/|a_pq|) with J[q,p] = -conj(b); updates are A <- J^H*A*J, V <- V*J. Eigenvalues are real (%i754), eigenvectors unitary (%cplx) -- matching NumPy eigh. Works on @cs (bloq 6) and @cd (bloq 7); asserts exact Hermitian. Adds complex scalar helpers (cadd/csub/cmul/cdiv/cconj/cre/cpak/cabs-re over /lib/complex) and the Hermitian arms (hermitian, off-norm-h, frob-h, diag-real, rot-cols-h, rot-rows-h, sweep-herm, eig-herm). Durability: replaced the unbounded /lib/math Newton sqrt with an iteration-capped fsqt (50-step backstop, width-fixed feps). A sub-ULP rtol otherwise makes Newton oscillate between two adjacent floats forever, which OOM-crashed the ship while developing this. This also hardens the A1 path (regression-checked: all 18 symmetric tests still green). Tests (tests/lib/saloon-eigh, generated by tools/eigh_check.py vs numpy.linalg.eigvalsh) cover @cs 2x2/3x3/4x4 and @cd 3x3/4x4 Hermitian matrices; sort(computed eigenvalues) = sort(eigh). All 5 pass. Caveats documented in EIG-DESIGN.md: rtol width must match the component width (@rs for @cs, @rd for @cd); the exact Hermitian assert rejects +-0.0 mismatches in conjugate pairs. Note: this branch stacks on %cplx (PR #46) and eig A1 (PR #47); rebase once those land. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Implements Phase A1 of the Saloon eigendecomposition design (
saloon/EIG-DESIGN.md): real symmetric eigenvalues and eigenvectors on Lagoon%i754arrays (f32/f64) via the cyclic Jacobi algorithm. This is the marquee numerical-linear-algebra capability Saloon was built toward, and the motivating consumer of the%cplxwork.New arms in a
+| %linalgsection of thesacore:+eig[vals=ray vecs=ray]— eigenvalues (1-D), eigenvectors (columns)+eigvalsray— eigenvalues only+eigvecsray— eigenvectors onlyplus building blocks: bloq-dispatched scalar float helpers, exact-symmetry check, off-diagonal / Frobenius norms,
diag, Givens row/column rotations, and a single cyclic sweep.Design choices
sqrt, converges quadratically. 60-sweep cap with a warning rather than a silent non-converged result.rs/rdmath doors at the array'sbloq; per-entryget-item/set-itemrotations (O(n) per pair) instead of fullmmul.Latent fix
Building Saloon against the current Lagoon surfaced a pre-existing bug:
trans-scalar/fun-scalarused exhaustive?- kind, but Lagoon'skindunion has since grown (%unum,%cplx). Both switched to?+ kind !!so unsupported kinds crash gracefully instead of failing to compile.Tests
tests/lib/saloon-eigis generated bytools/eig_check.pyagainstnumpy.linalg.eigh, covering the four standard sources of symmetric test matrices:MᵀMPer matrix it checks the order/sign-free invariant
A·V = V·diag(vals)andsort(computed) = sort(eigh). All 18 pass in ~20s on a live ship.Follow-ons (not in this PR)
%cplx) Jacobi.🤖 Generated with Claude Code