linalg: det/slogdet, rank-k pinv, solve_triangular; einsum ellipsis; dsp welch - #93
Merged
Conversation
…dsp welch det and slogdet ride the on-ANE pivoted LU (permutation parity times the U diagonal; slogdet in log form for fp16-hostile magnitudes). pinv is the rank-k pseudoinverse from randomized_svd with a relative cutoff, assembled through the ANE gemm. solve_triangular substitutes on-engine through the routed accessor (no inverse formed; large entries stay finite, tested past the slice-x16 saturation threshold). einsum gains ellipsis support: '...' expands to explicit right-aligned batch letters before the v1 parser runs, implicit outputs follow numpy's ellipsis-dims-first rule, an explicit output omitting a nonempty '...' raises like numpy, and mismatched ellipsis dims reject (size-1 broadcast has no expand lowering). dsp gains welch: PSD from the on-ANE stft, matching scipy.signal.welch(detrend=False) for density and spectrum scalings; nperseg must be a power of two because the staged FFT pads to one. BiCGSTAB (#82) is deliberately NOT included: implemented three ways (unrolled fp16; fp64 host recurrence over ANE matvecs; plus best-iterate tracking and outer residual correction) it plateaus at relerr 0.11 (cond 10) / 0.80 (cond 100) regardless of iteration budget, and a noise-injection simulation reproduces the failure - the BiCG short recurrences cannot tolerate the ~5e-4 relative noise of an fp16 matvec, where GMRES's full orthogonalization can. gmres remains the nonsymmetric solver. 51 tests pass on-device (M5) across test_linalg.py, test_einsum_ellipsis.py, and test_dsp_welch.py.
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
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.
Fixes #79, fixes #80, fixes #81, fixes #83, fixes #84.
Implements five of the seeded starter issues (the sixth of this group, #82 BiCGSTAB, is deliberately absent — see below). #76 was taken by @AnayGarodia in #92, which merged first.
linalg
det/slogdet(linalg: det and slogdet #79) — on-ANE pivoted LU, then permutation parity × diag(U) host-side;slogdetin log form so fp16-hostile magnitudes stay finite. Tested vs numpy at cond 10/100 including a scaled matrix whose raw product is fp16-hostile.pinv(linalg: pinv (rank-k pseudoinverse) #80) — rank-k pseudoinverse fromrandomized_svdwith a relative cutoff, assembled through the ANE gemm. Tested vsnp.linalg.pinvon an exactly-rank-k matrix plus the Moore–Penrose identityA P A = A.solve_triangular(linalg: solve_triangular #81) — substitution on-engine via the routed accessor (no inverse formed), both triangles; tested past the slice-×16 saturation threshold (entries ~5000 stay finite).einsum (#83)
...expands to explicit right-aligned batch letters before the v1 parser runs. Implicit outputs follow numpy's ellipsis-dims-first rule; an explicit output that omits a nonempty...raises exactly like numpy (verified against numpy's own error); mismatched ellipsis dims reject (EinsumUnsupported) since size-1 broadcast has no expand lowering. Batched matmul'...ij,...jk->...ik'and friends verified on-device againstnp.einsum.dsp (#84)
welch: PSD composed from the on-ANEstft, matchingscipy.signal.welch(detrend=False)for bothdensityandspectrumscalings at 2% relative; the two injected tones stand >10× out of the noise floor.npersegmust be a power of two — the staged FFT pads to one, which would silently shift the bins otherwise.Why no BiCGSTAB (#82)
Implemented three ways and measured: fully-unrolled fp16 diverges; an fp64 host recurrence over ANE matvecs plateaus at relerr 0.11 (cond 10) / 0.80 (cond 100) regardless of iteration budget; best-iterate tracking + outer residual correction doesn't move it. A noise-injection simulation reproduces the failure exactly — the identical code converges to 1e-10 with exact matvecs and falls apart at the ~5e-4 relative noise an fp16 matvec carries. BiCG-family short recurrences require matvec accuracy on the order of the target accuracy (inexact-Krylov theory, now measured on this hardware); GMRES's full orthogonalization tolerates the noise, and
gmresremains the module's nonsymmetric solver. #82 should close as documented-infeasible rather than merge a solver that quietly returns 10–80% error.Verification (on-device, M5)
51 tests across
test_linalg.py(full re-run, nothing regressed),test_einsum_ellipsis.py,test_dsp_welch.py. Off-device selection unchanged; ruff, 2-space pylint, pyright clean.