linalg

NumPy-style linear algebra on ndarray, with SIMD-accelerated contiguous kernels.

For the small-and-hot regime — a 3-D direction, a 4×4 transform built and consumed millions of times a second — reach for the sibling fixarray module: the same mathematics with the shape moved into the type (vec3f/mat4f, allocation-free, and faster than GLM across their whole overlap). linalg is the home of the heavy, shape-generic numerics below.

The routines mirror numpy.linalg and operate on ndarray::NDArray (2-D = matrix, 1-D = vector). The general eigensolvers eig/eigvals return a complex spectrum (CNDArray) — a real matrix can have complex conjugate eigenvalue pairs, e.g. a rotation has ±i — while the Hermitian solvers eigh/eigvalsh return a guaranteed-real spectrum (the same split as numpy). Kernels are compiled at -O3 -march=native so the hot loops auto-vectorize.

Usage

import linalg              # auto-links ndarray

let x = linalg.solve(A, b) # A·x = b
let d = linalg.det(A)

Functions

Every routine below returns a fresh array. Each also has an allocation-free out-parameter overload (solve(out, A, b), svd(u, s, vh, A), …) in routines.hpp for hot loops that reuse one scratch buffer per call — same math, caller-owned storage.

Products

• dot / vdot / inner — vector dot product (Σ aᵢbᵢ).

• outer — outer product of two vectors → matrix.

• matmul — matrix multiply.

• matrix_power — integer matrix power Aⁿ (negative via inv).

• kron — Kronecker (block) product.

Decompositions

• cholesky — lower-triangular L with A = L·Lᵀ (SPD only).

• qr — reduced QR via Householder reflections.

• svd — singular value decomposition (Golub–Reinsch: bidiagonalization + implicit QR).

• svdvals — singular values only (the SVD fast path, without forming U/Vᵀ).

Eigenvalues

• eig / eigvals — general square matrix (complex spectrum + eigenvectors; Hessenberg + shifted QR for the values, inverse iteration for the vectors).

• eigh / eigvalsh — symmetric or complex Hermitian matrix (real spectrum; real eigenvectors for a real symmetric matrix, complex eigenvectors for a Hermitian one; Householder tridiagonalization + QL, via a real 2n embedding for Hermitian input).

All eigenvalue routines return the spectrum descending (note: numpy's eigvalsh returns ascending), with eigenvector columns reordered to match.

# A real rotation matrix has complex eigenvalues ±i:
let r = ndarray.reshape(ndarray.array([0.0, -1.0, 1.0, 0.0]), [2, 2])
io.print(ndarray.to_string(linalg.eigvals(r)))   # [0+1j, 0-1j]

Complex inner-product spaces

Vectors and matrices can be complex (ndarray.complex(re, im)), so the routines work over complex inner-product spaces — Hermitian operators, complex wavefunctions:

• dot — bilinear product Σ aᵢbᵢ (complex, no conjugation; matches numpy).

• vdot — conjugate-linear Hermitian inner product ⟨a, b⟩ = Σ conj(aᵢ)·bᵢ (conjugates the first argument; vdot(a, a) is the real ‖a‖²).

• matmul — complex matrix multiply.

• conj_transpose — conjugate transpose (Hermitian adjoint) Aᴴ.

import io
import ndarray
import linalg

# A Hermitian operator H = [[2, 1+i], [1-i, 3]] — real eigenvalues 4, 1:
let re = ndarray.array([2.0, 1.0, 1.0, 3.0])
let im = ndarray.array([0.0, 1.0, -1.0, 0.0])
let H = ndarray.reshape(ndarray.complex(re, im), [2, 2])
io.print(ndarray.to_string(linalg.eigvalsh(H)))   # [4, 1]

# Hermitian inner product ⟨a, a⟩ = ‖a‖² is real:
let a = ndarray.complex(ndarray.array([1.0, 3.0]), ndarray.array([2.0, -1.0]))
io.print(linalg.vdot(a, a))                        # 15+0j

Norms & numbers

• norm — L2 (vector) / Frobenius (matrix).

• cond — 2-norm condition number.

• det / slogdet — determinant (LU); slogdet is overflow-safe.

• matrix_rank — numerical rank from SVD thresholding.

• trace — sum of the main diagonal.

Solving & inverses

• solve — solve A·x = b via LU with partial pivoting.

• lstsq — least-squares solution min‖A·x − b‖.

• inv — matrix inverse (LU).

• pinv — Moore–Penrose pseudo-inverse (SVD, any shape).

SIMD

SIMD here is pure compiler auto-vectorization (no intrinsics): the kernels are contiguous, unit-stride loops compiled at -O3 -march=native. These functions only report the build's capability:

• simd_features — instruction sets this build targets (e.g. AVX2;FMA, NEON, scalar).

• simd_lane_doubles — widest SIMD lane width in doubles.

On a build with no SIMD every routine returns identical results, just scalar / slower — SIMD is never a correctness dependency. The full model (and the compile-time-dispatch limitation) is documented in simd.hpp.

Performance vs NumPy

linalg goes head-to-head with NumPy, whose array ops dispatch to BLAS/LAPACK — hand-tuned, vectorized, often multi-threaded Fortran. The scripts/numpy_compare.py harness feeds the same fixed-seed, well-conditioned matrix to both, runs the same op many times with the result consumed, and checks the answers agree. Each function's Performance row above carries its own measurement; the full size-dependence:

What we compared against — read this before the numbers. NumPy's absolute speed, and where every crossover lands, depends far more on which BLAS it links than on which NumPy version it is. The generated stamp below records the resolved library rather than the version string, because "NumPy 1.26.4" does not identify a measurement and the library does.

On this machine that resolves to libblas.so.3 → the reference implementation, not OpenBLAS and not MKL. That matters for how much these wins are worth: reference BLAS is the unoptimized baseline, so a cheatah win over it is a win over untuned Fortran, not over the tuned kernels most NumPy installs actually use. A tuned BLAS would narrow the large-n rows and push the crossovers lower. The small-n wins are a different claim and survive either way — they come from cheatah having no per-call Python and dispatch overhead to pay, which no choice of BLAS changes.

The Eigen comparison is a separate measurement, in the native Google Benchmark harness (tests/benchmarks/eigen_compare_bench.cpp), where cheatah and Eigen 3.4 are both compiled C++ timed identically on one thread — an apples-to-apples per-core comparison.

Two harnesses, two tables — deliberately. These numbers used to share one table, with the cheatah column measured by numpy_compare.py (separate processes) sitting beside an Eigen column measured by Google Benchmark (compiled C++, one process), and a prose warning not to read across. A warning is a weaker fix than a structure: a reader who divides two adjacent columns gets a wrong answer no matter how the caption is worded. Each harness now publishes only what it measured, and neither table contains a column it did not produce.

Both are generated — see docs/performance.md for the methodology (striated, interleaved, medians with dispersion) and the reference machine.

vs NumPy

vs Eigen

Tally (a difference counts only above 1.15x AND 0.25 ns) — vs eigen: 6 faster / 1 parity / 1 slower.

• Loss vs eigen: BM_dot/64 — cheatah 10.18 ns vs 6.69 ns (1.52x slower)

After a focused optimization round (hunting a few recurring mistakes across every routine — a heap allocation in a hot predicate, single-accumulator reductions, a column-stride QR walk, and a result buffer that was zero-filled and then thrown away):

• vs NumPy/LAPACK, cheatah now wins across nearly the whole library — products, the LU family, the SVD, the symmetric eigensolver, outer, qr, kron, and large norm (the last four previously lost). The handful still behind (svdvals/cond/ matrix_rank ~1.1×) are SVD-threshold queries where LAPACK's bidiagonal solver edges it.

• vs Eigen 3.4 on one core, cheatah matches or beats it on the bulk — inv (≈1.7×), svd (≈1.6×), det/matmul/trace (≈1.3–1.4×), outer (≈1.2–1.3×, now that the result buffer is built uninitialized and moved in zero-copy), solve/eigvalsh/eigh/ norm (≈1.1–1.2×). Eigen still leads on its blocked BLAS-3 kernels — qr (1.3–1.8×), cholesky (1.2×), and the eigh eigenvector path (1.1×) — which we flag honestly rather than hide.

cheatah does all of this on one core, by design — single-threaded is a feature, not a shortfall (no hidden threads, no contention, nothing to tune). Both NumPy's and Eigen's remaining edges are the very large or blocked dense problems where threaded/BLAS-3 kernels spread the work — a different operating point. See the Performance guide for the full rationale.

Per-function docs (parameters, complexity, heap behavior) are in routines.hpp. Tested in ../tests/linalg_routines_test.cpp and ../tests/linalg_smoke_test.cpp; ASan + Valgrind clean via the QA gate (security/run-valgrind.sh).