feat(linalg): numpy.linalg in pure Go

[tannevaled]

Package linalg: numpy.linalg for ndarray arrays, in pure Go (CGO=0), over the package's packed SIMD GEMM.

What it has

Solve, Inv, Det, Slogdet, Cholesky (lower/upper), QR (reduced, complete, r, raw), SVD/SVDVals, Eigh/Eigvalsh (UPLO), Eig/Eigvals (always complex, as in numpy 2), Lstsq/LstsqRcond, Pinv/PinvRcond, MatrixRank/MatrixRankTol, Norm (all vector and matrix orders, axes, keepdims), Cond, MatrixPower, MultiDot.

• Stacked (..., M, N) inputs, one matrix per goroutine; Solve broadcasts like numpy's gufunc.
• Every dtype. As in numpy, the computation is in double; float32/complex64 results are rounded back, ints/bools promote to float64.
• Errors: ErrLinAlg (numpy's LinAlgError), which ErrSingular, ErrNotPosDef and ErrNoConvergence wrap; ErrShape; ErrArgs.

Algorithms (LAPACK's, generic over float64/complex128)

• Recursive LU (getrf2) and recursive Cholesky, so the work goes to the GEMM.
• Blocked Householder QR (compact WY).
• Hermitian eigenproblems: blocked tridiagonalization (latrd) with divide and conquer (Cuppen, Gu–Eisenstat) above order 25, and implicit QL/QR below that and for eigenvalues only.
• SVD: blocked Golub–Kahan bidiagonalization (labrd) with bidiagonal divide and conquer (as gesdd) above order 25, and Demmel–Kahan implicit QR (bdsqr) below that and for values only.
• General eigenproblems: gebal, Hessenberg reduction, Francis double-shift QR (single-shift for complex input), back-substitution with overflow scaling, gebak, and dgeev's normalization.

docs/linalg.md explains the choice of Golub–Kahan with QR over one-sided Jacobi: it has the same backward-error class as numpy's gesdd, at a fraction of Jacobi's sweeps.

Verification

• Fixture. 1,040 numpy 2.5.3 cases are embedded (linalg/testdata/cases.py). 914 are checked on their values and 124 must raise the same kind of error. The other 2 are logged as deliberate differences, listed in the docs.
• Tolerances. Unique results must agree within 32·n·ε·κ. Factors that are not unique are checked through their identities: orthogonality, residuals, triangularity, the eigenvalue order, dgeev's normalization and exact conjugate pairs. Eig's eigenvalues are checked against an empirical sensitivity stored in the fixture, which is how defective matrices are covered.
• Extra tests. Orders 100/131, parallel stacks, degenerate spectra, the hard Eig matrices, forced non-convergence, and unit tests of the 2×2 kernels.
• Negative controls.
  • Dropping the LU pivot swap fails 64 fixture cases (23 solve, 17 inv, 10 det, 10 slogdet, 4 matrix_power) plus the order-100 identity test.
  • Loosening the secular-equation stop test by 1e8 fails 4 eigh cases plus the bidiagonal divide-and-conquer test.
• Coverage. 100.0% under the CI gate (-race -coverpkg).
• Checks. go vet (native, 386, wasm), the s390x test build and gofmt are clean.

Speed against numpy (Accelerate) — numpy wins on single matrices

Medians of three interleaved rounds on a shared host (load 16–23 on 16 cores). Go time ÷ numpy time at n = 64 / 256 / 1024:

op	64	256	1024
solve	2.1	3.9	1.9
inv	3.7	3.8	1.7
cholesky	3.3	2.2	1.7
qr	1.8	1.5	1.5
svd	2.6	2.5	2.1
eigh	1.9	2.1	1.7
eig	1.0	2.8	7.1

• Most of the gap is the GEMM itself: 18 ms vs Accelerate's 3.8 ms at 1024².
• eig lacks LAPACK's multishift QR with aggressive early deflation.
• On stacks of small matrices this package is faster, because numpy loops over them on one core: 1.7–11× on inv/det/eigh/svd/eig over 10000×4×4 and 1000×32×32.

Not done

• Multishift QR with AED and a blocked Hessenberg reduction, for eig at large n.
• A small-matrix fast path.
• The hermitian= and rtol= options.
• tensorsolve/tensorinv/cross/trace/vecdot and friends.
• Eigh on NaN input returns ErrNoConvergence, where numpy returns NaNs.

🤖 Generated with Claude Code