feat(linalg): numpy.linalg in pure Go

README.md

@@ -32,7 +32,14 @@
 - **Indexing** — basic-indexing `Slice` views plus **boolean/fancy indexing**:
   `MaskSelect` (`a[mask]`), `Nonzero` (`flatnonzero`), and `Take`.
 - **Manipulation** — `Concatenate`/`Stack`/`VStack`/`HStack`.
-- **Linear algebra** — `MatMul`/`Dot`/`Inner`/`Outer`.
+- **Linear algebra** — `MatMul`/`Dot`/`Inner`/`Outer`; package
+  [`linalg`](linalg) is `numpy.linalg` in pure Go: `Solve`, `Inv`, `Det`,
+  `Slogdet`, `Cholesky`, `QR` (all four modes), `SVD`/`SVDVals`,
+  `Eigh`/`Eigvalsh`, `Eig`/`Eigvals`, `Lstsq`, `Pinv`, `MatrixRank`, `Norm`,
+  `Cond`, `MatrixPower`, `MultiDot` — stacked `(..., M, N)` inputs (run in
+  parallel), every dtype, LAPACK's algorithms (divide and conquer for the
+  SVD and Hermitian eigenproblems), checked against numpy 2.5.3 on 1,040
+  cases. See [docs/linalg.md](docs/linalg.md) for accuracy and speed.
 - **Fourier transforms** — package [`fft`](fft) is `numpy.fft` on arrays
   (`FFT`/`RFFT`/`FFT2`/`FFTN`/… with `n`/`s`, `axes` and `norm`, single
   precision kept), on [go-fft](https://github.com/go-fft/fft).
@@ -229,6 +236,8 @@ u, _ := random.Rand(3)                              // [0.5488135 0.71518937 0.6
 b, _ := ndarray.Arange(0, 6, 1)
 b, _ = b.Reshape(3, 2)
 prod, _ := m.MatMul(b)              // (2,3)·(3,2) -> (2,2)
+inv, _ := linalg.Inv(prod)          // np.linalg.inv
+left, sv, vh, _ := linalg.SVD(m, false) // np.linalg.svd(m, full_matrices=False)
 ```
 
 ## License

docs/bench_linalg_numpy.py

@@ -0,0 +1,32 @@
+# numpy side of docs/linalg.md: uv run --with numpy==2.5.3 python docs/bench_linalg_numpy.py
+import time, sys, numpy as np, numpy.linalg as la
+def mat(n, kind):
+    r = np.random.default_rng(n)
+    a = r.standard_normal((n, n))
+    if kind == "spd": return a @ a.T + n * np.eye(n)
+    if kind == "sym": return a + a.T
+    return a
+ops = [("solve", "", lambda a: la.solve(a, np.ones(a.shape[0]))), ("inv", "", la.inv),
+       ("cholesky", "spd", la.cholesky), ("qr", "", la.qr), ("svd", "", lambda a: la.svd(a, full_matrices=False)),
+       ("eigh", "sym", la.eigh), ("eig", "", la.eig)]
+sizes = [int(x) for x in sys.argv[1].split(",")] if len(sys.argv) > 1 else [64, 256, 1024]
+for name, kind, f in ops:
+    for n in sizes:
+        a = mat(n, kind)
+        f(a)
+        reps = max(1, int(0.3 / max(1e-6, (lambda t0: (f(a), time.perf_counter() - t0)[1])(time.perf_counter()))))
+        reps = min(reps, 200)
+        ts = []
+        for _ in range(reps):
+            t0 = time.perf_counter(); f(a); ts.append(time.perf_counter() - t0)
+        print(f"{name}/{n} {np.median(ts)*1e3:.3f} ms")
+
+# Stacks of small matrices.
+rs=np.random.default_rng(9)
+for cnt,n in [(10000,4),(1000,32)]:
+    a=rs.standard_normal((cnt,n,n)); s=a+np.swapaxes(a,-1,-2)
+    for name,f in [("inv",lambda: la.inv(a)),("det",lambda: la.det(a)),("eigh",lambda: la.eigh(s)),("svd",lambda: la.svd(a,full_matrices=False)),("eig",lambda: la.eig(a))]:
+        f(); ts=[]
+        for _ in range(7):
+            t=time.perf_counter(); f(); ts.append(time.perf_counter()-t)
+        print(f"{name}/{cnt}x{n} {np.median(ts)*1e3:.3f} ms")

docs/linalg.md

@@ -0,0 +1,256 @@
+# linalg — numpy.linalg in pure Go
+
+Package [`linalg`](../linalg) gives ndarray arrays the functions of
+`numpy.linalg`, written in Go without cgo or assembly of its own: the
+matrix products it relies on are the ndarray package's packed SIMD GEMM.
+
+| numpy.linalg | linalg |
+|---|---|
+| `solve(a, b)` | `Solve(a, b)` |
+| `inv(a)` | `Inv(a)` |
+| `det(a)`, `slogdet(a)` | `Det(a)`, `Slogdet(a)` |
+| `cholesky(a, upper=…)` | `Cholesky(a, Lower \| Upper)` |
+| `qr(a, mode)` | `QR(a, ModeReduced \| ModeComplete \| ModeR \| ModeRaw)` |
+| `svd(a, full_matrices)`, `svdvals(a)` | `SVD(a, full)`, `SVDVals(a)` |
+| `eigh(a, UPLO)`, `eigvalsh(a, UPLO)` | `Eigh(a, uplo)`, `Eigvalsh(a, uplo)` |
+| `eig(a)`, `eigvals(a)` | `Eig(a)`, `Eigvals(a)` |
+| `lstsq(a, b, rcond)` | `Lstsq(a, b)`, `LstsqRcond(a, b, rcond)` → `LstsqResult{X, Residuals, Rank, S}` |
+| `pinv(a, rcond)` | `Pinv(a)`, `PinvRcond(a, rcond)` |
+| `matrix_rank(A, tol)` | `MatrixRank(a)`, `MatrixRankTol(a, tol)` |
+| `norm(x, ord, axis, keepdims)` | `Norm(x, ord, axes, keepdims)`; `ord` is `Ord{}` (None), `P(p)`, `Fro`, `Nuc` |
+| `cond(x, p)` | `Cond(x, p)` |
+| `matrix_power(a, n)` | `MatrixPower(a, n)` |
+| `multi_dot(arrays)` | `MultiDot(arrays...)` |
+| `LinAlgError` | `ErrLinAlg`, wrapped by `ErrSingular`, `ErrNotPosDef`, `ErrNoConvergence` |
+
+Shape problems give `ErrShape` (numpy raises `LinAlgError` for some of them
+and `ValueError` for others), bad options `ErrArgs`.
+
+## Semantics
+
+- **Stacks.** The last two axes are the matrix; leading axes are a stack,
+  processed one matrix per goroutine. `Solve` broadcasts the stack
+  dimensions of `a` and `b` as numpy's gufunc does; a 1-D `b` is one
+  vector, any other `b` a stack of matrices (numpy 2's rule).
+- **Dtypes.** As in numpy, everything is computed in double precision:
+  float64 and complex128 directly, integers and bools promoted to float64,
+  float32 and complex64 widened. Results have numpy's dtypes: single
+  precision in, single precision out (rounded from the double result, which
+  is what numpy does too — its linalg gufuncs also run in double); singular
+  values, Hermitian eigenvalues, `logabsdet` and norms are real; `Eig` is
+  always complex (numpy 2). `MatrixPower` with n ≥ 0 keeps an integer dtype.
+- **Conventions.** LAPACK's: Householder reflectors with β = −sign(α)‖x‖, so
+  R's diagonal can be negative as in numpy; ascending Hermitian eigenvalues;
+  decreasing singular values; `Eig` vectors with unit norm and their largest
+  component real, conjugate vectors for conjugate eigenvalues of a real
+  matrix (dgeev's normalization).
+- **Not unique.** Signs and phases of eigenvectors and singular vectors, and
+  any orthonormal basis of a repeated eigenvalue's space, may differ from
+  numpy's. Compare `A·v = λ·v`, `A = U·S·Vh`, not the vectors.
+
+## Algorithms, and why
+
+All are LAPACK's, re-expressed for row-major storage. Each is written once,
+generic over `float64 | complex128`; Go 1.27 stencils the two and folds the
+type switches of the helpers (`conj`, `abs`, …), so the real code pays
+nothing for the complex case (measured: same speed as a hand-specialised
+loop).
+
+| operation | algorithm | where the time goes |
+|---|---|---|
+| `Solve`, `Inv`, `Det`, `Slogdet` | LU with partial pivoting, recursive (getrf2): halves the columns, so the trailing update and the triangular solves are GEMMs | GEMM |
+| `Cholesky` | recursive potrf: `L21 = A21·L11⁻ᴴ`, `A22 −= L21·L21ᴴ` | GEMM |
+| `QR` | blocked Householder (geqrf/orgqr, panels of 32, compact WY `I − V·T·Vᴴ`) | GEMM |
+| `Eigh` | blocked tridiagonalization (latrd, panels of 32: half the flops in GEMM, half in matrix-vector products) + divide and conquer (stedc: Cuppen's split, deflation, secular equation, Gu–Eisenstat's recomputed vector) above order 25, implicit QL/QR (steqr) below | tridiagonalization |
+| `Eigvalsh` | the same reduction + QL/QR without vectors | tridiagonalization |
+| `SVD` | blocked Golub–Kahan bidiagonalization (labrd) + divide and conquer on the bidiagonal (bdsdc: Gu and Eisenstat's arrow matrix and secular equation) above order 25, implicit zero-shift/shifted QR (bdsqr, Demmel–Kahan) below; wide matrices factor Aᴴ | bidiagonalization |
+| `SVDVals`, `MatrixRank`, `Cond` (2-norm), `Norm` (2, nuc) | the same reduction + bdsqr without vectors | bidiagonalization |
+| `Eig`, `Eigvals` | balancing (gebal: permutations and powers of two), Hessenberg reduction (gehd2), Francis double-shift QR (dlahqr; single-shift zlahqr for complex input), eigenvectors by back-substitution on the Schur form with overflow scaling (trevc), back-transformed by one GEMM, gebak, dgeev's normalization | QR iteration |
+| `Lstsq` | SVD, `x = V·S⁺·Uᴴ·b` over singular values above `rcond·σ₁`; the residuals from the Householder QR of `a` (the tail of `Qᴴb`, free of the cancellation in `b − a·x`) | SVD |
+| `Pinv` | SVD, `V·S⁺·Uᴴ` | SVD |
+| `Norm`, `MatrixPower`, `MultiDot` | direct; 2-norms rescaled against overflow; repeated squaring in numpy's order; numpy's matrix-chain dynamic programme | — |
+
+**The SVD: why Golub–Kahan with QR, and divide and conquer.** The choice
+was between Golub–Kahan bidiagonalization followed by the implicit-shift QR
+iteration, and one-sided Jacobi. One-sided Jacobi is the more accurate in a
+strong sense — it computes the small singular values of a matrix whose
+columns are badly scaled to high relative accuracy (Demmel–Veselić) — but it
+needs 5–10 sweeps of O(n³) rotations, several times the work, and it is not
+what numpy computes. Golub–Kahan reduction is backward stable: the computed
+singular values are those of `A + E` with `‖E‖ ≈ ε‖A‖`, so each is within
+about `ε·σ₁` of the truth — numpy's (gesdd's) accuracy, exactly. The QR
+iteration (bdsqr, with Demmel and Kahan's zero shift) then finds the
+singular values of the bidiagonal to high relative accuracy, so it adds
+nothing to that error. It is used at the leaves of the divide and conquer
+and whenever the vectors are not wanted. With vectors, the QR iteration's
+plane rotations cost O(n³) of non-vectorized work (5.8 s at n = 1024 in the
+first version); divide and conquer, LAPACK's gesdd algorithm and the one
+numpy calls, puts that work into GEMMs (0.26 s now). Its singular vectors
+stay orthogonal because the updating vector is recomputed from the computed
+roots (Gu and Eisenstat), and the secular equation is solved with each
+distance `dᵢ − σ` formed from the nearest pole, in the squared form
+`(dᵢ − σ)(dᵢ + σ)`.
+
+**Hermitian eigenproblems** follow the same pattern: the tridiagonal QL/QR
+iteration (steqr) is used without vectors and for orders up to 25, divide
+and conquer (stedc) above, which took n = 1024 from 3.6 s to 0.13 s.
+
+**General eigenproblems** use the classical double-shift QR (dlahqr), not
+LAPACK's multishift QR with aggressive early deflation (dlaqr0–5). That is
+the largest gap to numpy, at large orders (below).
+
+## Accuracy and verification
+
+**Against numpy.** `linalg/testdata/cases.py` generates 1,040 cases with
+numpy 2.5.3 (`uv run --with numpy python linalg/testdata/cases.py >
+linalg/testdata/cases.json`; the JSON is embedded in the test binary, so the
+emulated CI architectures read it from anywhere). They cover every
+function; float64, float32, complex128, complex64, int32/int64 and bool
+inputs; 0×0, 1×1, tall, wide, stacked and broadcast shapes; Hilbert
+matrices (κ up to 10¹⁶), graded and rank-deficient matrices, exactly
+singular ones, non-symmetric matrices with complex eigenvalues, Jordan
+blocks and other defective matrices, NaN and Inf inputs; and every error
+numpy raises on them (124 cases). The test checks, for each case:
+
+- errors: an error of the same kind (`ErrSingular`, `ErrNotPosDef`,
+  `ErrNoConvergence`, `ErrLinAlg` for non-finite input to `Eig`,
+  `ErrShape`, `ErrArgs`);
+- dtype and shape of every result, exactly;
+- unique values (solutions, inverses, determinants, Cholesky factors,
+  singular values, eigenvalues, least-squares solutions and residuals,
+  pseudo-inverses, ranks, norms, condition numbers, powers, products)
+  elementwise, within `32·n·ε·κ` of the largest entry of the result, with
+  `κ` the condition number of the problem (1 for singular values, Hermitian
+  eigenvalues and norms; `cond(A)` for solves and inverses; the effective
+  condition over the kept singular values for `Lstsq` and `Pinv`, squared for
+  the residuals), plus one float32 rounding for single-precision results;
+- factors that are not unique through their identities: `‖QᴴQ − I‖`,
+  `‖A − QR‖/‖A‖`, R upper triangular; `‖UᴴU − I‖`, `‖VhVhᴴ − I‖`,
+  `‖A − U·S·Vh‖/‖A‖`; `‖VᴴV − I‖`, `‖AV − VΛ‖/‖A‖` with ascending
+  eigenvalues; for `Eig`, `‖Av − λv‖/‖A‖` per vector, unit norm, largest
+  component real, conjugate pairs exactly conjugate — all within
+  `32·n·ε` (`ε` of the result's dtype);
+- general eigenvalues matched to numpy's nearest unused one, within 100
+  times how far numpy's own eigenvalues move when each entry is perturbed by
+  10⁻¹⁴ of itself (an empirical condition number, recorded in the fixture),
+  so a Jordan block's eigenvalue, accurate only to `ε^(1/4)`, is judged by
+  what is achievable.
+
+914 cases are value-checked and 124 are errors; the 2 deliberate departures
+below are logged.
+
+**Beyond the fixture** (`properties_test.go`, `kernels_test.go`): the same
+identities at orders 100 and 131 (above the blocking and divide-and-conquer
+thresholds), real and complex, tall and wide; a stack of 16 matrices run in
+parallel equals the matrices run one by one; repeated eigenvalues and
+singular values, rank-deficient and zero matrices (the deflation paths);
+cyclic permutations, Grcar and Frank matrices, orthogonal matrices,
+high-order Jordan blocks and coupled rotation blocks (exceptional shifts,
+overflow rescaling in the back-substitution); balancing at the limits of the
+exponent range; every iteration forced to fail (`ErrNoConvergence` from
+each function); and the 2×2 kernels (laev2, las2, lasv2, lanv2, lartg,
+solve2) checked through their defining identities on the inputs that reach
+their special cases.
+
+**Negative controls.** Dropping the row interchange in the LU leaf (no
+partial pivoting) fails 64 fixture cases — 23 `solve`, 17 `inv`, 10 `det`,
+10 `slogdet`, 4 `matrix_power` — and the order-100 identity test. Loosening
+the secular-equation solver's stopping test by 10⁸ fails the four order-34
+`eigh` cases (eigenvalues off by ~10⁻¹¹ relative) and the zero-diagonal
+bidiagonal divide-and-conquer test (`‖B − USVᵀ‖` of 10⁻⁸).
+
+**Coverage** is 100% of statements under the CI gate (`go test -race
+-coverpkg=./... ./...`).
+
+## Speed against numpy
+
+Median wall time of `BenchmarkLinalg` against numpy 2.5.3 (`la.solve(a,
+ones)`, `la.inv`, `la.cholesky` of `A·Aᵀ + nI`, `la.qr` reduced,
+`la.svd(full_matrices=False)`, `la.eigh` of `A + Aᵀ`, `la.eig`), n×n
+standard normal matrices, Apple M-series (16 cores, arm64), three
+interleaved rounds. numpy's wheel uses **Apple Accelerate**. The host was
+shared and heavily loaded (load average 16–23 on 16 cores during the runs),
+so absolute numbers are noisy; the ratios held across rounds.
+
+| op | n | Go (ms) | numpy (ms) | Go / numpy |
+|---|---|---|---|---|
+| solve | 64 | 0.046 | 0.022 | 2.1 |
+| solve | 256 | 0.99 | 0.25 | 3.9 |
+| solve | 1024 | 21.2 | 10.9 | 1.9 |
+| inv | 64 | 0.11 | 0.029 | 3.7 |
+| inv | 256 | 2.19 | 0.58 | 3.8 |
+| inv | 1024 | 39.4 | 23.7 | 1.7 |
+| cholesky | 64 | 0.030 | 0.009 | 3.3 |
+| cholesky | 256 | 0.54 | 0.24 | 2.2 |
+| cholesky | 1024 | 10.9 | 6.6 | 1.7 |
+| qr | 64 | 0.15 | 0.083 | 1.8 |
+| qr | 256 | 2.32 | 1.57 | 1.5 |
+| qr | 1024 | 54.7 | 37.3 | 1.5 |
+| svd | 64 | 0.78 | 0.30 | 2.6 |
+| svd | 256 | 13.3 | 5.25 | 2.5 |
+| svd | 1024 | 263 | 123 | 2.1 |
+| eigh | 64 | 0.41 | 0.22 | 1.9 |
+| eigh | 256 | 6.54 | 3.12 | 2.1 |
+| eigh | 1024 | 130 | 78 | 1.7 |
+| eig | 64 | 1.07 | 1.06 | 1.0 |
+| eig | 256 | 55 | 19.4 | 2.8 |
+| eig | 1024 | 5060 | 712 | 7.1 |
+
+**numpy is faster on single matrices**, by 1.5–4× (7× for `eig` at 1024).
+Most of that is the matrix product underneath: on the same host the ndarray
+GEMM takes 18 ms for a 1024² product where Accelerate (which uses the
+matrix coprocessor) takes 3.8 ms, 0.37 ms against 0.10 ms at 256. The
+blocked algorithms here spend most of their time in that GEMM, so their gap
+to numpy is smaller than the GEMM's. The exceptions:
+
+- `eig` at large n: the classical double-shift QR is O(10n³) of scalar
+  work; LAPACK's multishift QR with aggressive early deflation moves it into
+  GEMMs. This is the main follow-up.
+- small orders (64): fixed costs (conversion of the input, allocation,
+  the packing of a small GEMM) weigh more than the arithmetic.
+
+**Stacks of small matrices run in parallel** — numpy's gufuncs loop over
+them on one core — and there this package is faster (median of 3, same
+host):
+
+| op | stack | Go (ms) | numpy (ms) | numpy / Go |
+|---|---|---|---|---|
+| inv | 10000 × 4×4 | 0.94 | 3.83 | 4.1 |
+| det | 10000 × 4×4 | 0.81 | 2.12 | 2.6 |
+| eigh | 10000 × 4×4 | 3.0 | 13.8 | 4.6 |
+| svd | 10000 × 4×4 | 5.3 | 29.5 | 5.6 |
+| eig | 10000 × 4×4 | 4.9 | 55.7 | 11 |
+| inv | 1000 × 32×32 | 3.4 | 8.5 | 2.5 |
+| det | 1000 × 32×32 | 1.6 | 4.0 | 2.5 |
+| eigh | 1000 × 32×32 | 21.3 | 50.4 | 2.4 |
+| svd | 1000 × 32×32 | 44.3 | 75.9 | 1.7 |
+| eig | 1000 × 32×32 | 29.7 | 232 | 7.8 |
+
+Reproduce: `go test -run XXX -bench 'Linalg|Stacked' ./linalg`, and the
+numpy side with `uv run --with numpy==2.5.3 python docs/bench_linalg_numpy.py`
+(the same sizes and distributions; the values differ, being from another
+generator).
+
+## Differences from numpy
+
+| numpy | here | why |
+|---|---|---|
+| `det` is `sign·exp(logabsdet)` | the product of the pivots, accumulated as mantissa and exponent | exact to rounding, no spurious overflow or underflow of intermediate products, and no `|logdet|·ε` error from `exp` |
+| `norm` of a vector with entries near 10²⁰⁰ overflows to `inf` | the 2-norms are rescaled (`1.414e200` for `[1e200, 1e200]`) | the true value is representable |
+| float32 norms accumulate in float32 | in float64, then rounded | more accurate; within a float32 rounding of numpy |
+| `cond(singular, 'nuc')` raises `LinAlgError` (the SVD of its NaN "inverse" fails) | `+Inf`, as for every other order | a singular matrix has infinite condition number |
+| `eigh` of a matrix with a NaN returns NaN eigenvalues | `ErrNoConvergence` | the iteration cannot proceed; numpy's `svd` and `eig` refuse NaN too |
+| `LinAlgError` for a wrong rank or a non-square matrix | `ErrShape` | a shape error, not a numerical failure |
+| `svd(…, hermitian=True)`, `pinv(…, hermitian=True)`, `matrix_rank(…, hermitian=True)`, `matrix_rank(rtol=)`, `pinv(rtol=)` | not provided (`PinvRcond(a, max(M,N)·ε)` is `rtol=None`) | follow-ups |
+| `tensorsolve`, `tensorinv`, `cross`, `outer`, `trace`, `diagonal`, `matmul`, `vecdot`, `matrix_norm`, `vector_norm`, `matrix_transpose` | not in this package (`ndarray` has `MatMul`, `Outer`, `Transpose`; `Norm` covers both norms) | follow-ups |
+| `qr` modes `'full'` and `'economic'` (deprecated) | not provided | deprecated in numpy |
+
+## What is left
+
+- Multishift QR with aggressive early deflation for `Eig` (the 7× gap at
+  n = 1024), and a blocked Hessenberg reduction.
+- A faster small-matrix path (n ≤ 64): unblocked kernels without the GEMM
+  packing, fewer copies.
+- The `hermitian` options and `rtol` arguments above, and the tensor
+  functions.
+- `Eigh` on NaN input returning NaNs like numpy, if that is preferred.

internal/kernels/gemm.go

@@ -434,10 +434,15 @@ func dotRange(a, b []float64) float64 {
 // is pure unit-stride streaming — far cheaper than routing a single-column GEMM
 // through the packer. dst need not be pre-zeroed (each entry is written, not
 // accumulated). Rows are disjoint, so the parallel result equals the serial one.
-func MatVecP(dst, a, v []float64, m, k int) {
+func MatVecP(dst, a, v []float64, m, k int) { MatVecStridedP(dst, a, k, v, m, k) }
+
+// MatVecStridedP is MatVecP on rows lda elements apart: row i of a is
+// a[i*lda : i*lda+k], so a sub-block of a larger row-major matrix is
+// multiplied in place.
+func MatVecStridedP(dst, a []float64, lda int, v []float64, m, k int) {
 	body := func(lo, hi int) {
 		for i := lo; i < hi; i++ {
-			dst[i] = dotSIMD(a[i*k:i*k+k], v)
+			dst[i] = dotSIMD(a[i*lda:i*lda+k], v)
 		}
 	}
 	if m*k < matVecThreshold {