cheatah
Source BigBrain
Source

stdlib/linalg/routines.hpp

1// Copyright (c) 2026 BigBrain LLC. MIT-licensed (see LICENSE).
2// Original work; see ACKNOWLEDGMENTS.md for the open-source ideas we build upon.
3#pragma once
5/**
6 * @file routines.hpp
7 * @brief cheatah `linalg` — numpy's linear-algebra API on ndarray, with
8 * SIMD-friendly contiguous kernels (a .purr program writes `linalg.solve(A, b)`).
9 *
10 * `import linalg` to use it (auto-links ndarray). Unit tests:
11 * stdlib/tests/linalg_routines_test.cpp; SIMD reporting tested in
12 * stdlib/tests/linalg_smoke_test.cpp. The suite runs under AddressSanitizer (the `asan`
13 * preset) and Valgrind (security/run-valgrind.sh) on every QA-gate run.
14 *
15 * Routines operate on `ndarray::NDArray` (2-D = matrix, 1-D = vector). They mirror
16 * https://numpy.org/doc/stable/reference/routines.linalg.html and are implemented
17 * in routines.cpp (LU w/ partial pivoting, Cholesky, Householder QR, one-sided
18 * Golub–Reinsch SVD, Householder-tridiagonal + QL symmetric eigen, Hessenberg + shifted-QR general
19 * eigen) at -O3 -march=native so the hot loops auto-vectorize.
20 *
21 * SIMD here is pure compiler auto-vectorization (no intrinsics). On a scalar build
22 * (no vector ISA) every routine still returns identical results, just slower — see
23 * simd.hpp's file comment for the full SIMD model, the no-SIMD behavior, and the
24 * compile-time-dispatch limitation.
25 *
26 * @note `n` below is the matrix dimension. The general eigensolvers `eig`/`eigvals`
27 * return a **complex** spectrum (@ref CNDArray) — a real matrix can have
28 * complex conjugate eigenvalue pairs — while the Hermitian solvers
29 * `eigh`/`eigvalsh` return a guaranteed-real spectrum. The LU/SVD-based scalar
30 * routines (`det`/`slogdet`/`cond`/`matrix_rank`) allocate scratch O(n²) for the
31 * factorization even though they return a scalar; the products and reductions
32 * (`dot`/`matmul`/`trace`/`norm`/…) read their operands in place — zero-copy
33 * when contiguous, packing a strided view once.
34 *
35 * @concurrency every routine is DELIBERATELY single-threaded (no hidden thread pool —
36 * the fastest-per-core contract); the caller composes parallelism across
37 * independent problems.
38 */
39#include <complex>
40#include <vector>
42#include "backend.hpp"
43#include "concepts.hpp"
44#include "enums.hpp"
45#include "ndarray.hpp"
47namespace cheatah::linalg {
49// The routines operate on cheatah::ndarray::NDArray, re-exported unqualified for brevity in the
50// signatures below. The directive below hides this re-export from the API doc generator so it
51// does not emit a phantom duplicate cheatah::linalg::NDArray class in the namespace/XML structure.
52/// @cond INTERNAL
53using ndarray::NDArray;
54/// @endcond
56/// A complex scalar (`std::complex<double>`) — the element type of @ref CNDArray and
57/// the return type of the complex inner products @ref dot / @ref vdot.
58using Cplx = std::complex<double>;
60/// A complex array (`basic_ndarray<std::complex<double>>`) — what the general
61/// eigensolvers return, since a real matrix can have complex eigenvalues. Prints
62/// element-wise as `a+bj` via @ref cheatah::ndarray::to_string.
63using CNDArray = ndarray::basic_ndarray<Cplx>;
65// ---- Matrix and vector products ----
66// Dot / vdot / inner — the scalar reductions — live in backend.hpp as the scalar-out kernel
67// pattern: an allocating front `T dot(a, b)` plus a `void dot(out, a, b)` kernel split by the
68// HostArray/DeviceArray concepts (one pair serving real, complex, host, and — via a device
69// extension — device operands). Same for `trace`.
70// Outer — both the allocating front `outer(a,b)` and the out-parameter kernel `outer(out,a,b)`
71// are the two-layer overload pair in backend.hpp (the matmul pattern).
72// Matmul — both the allocating front `matmul(a,b)` and the out-parameter kernel `matmul(out,a,b)`
73// are the two-layer `template <Field T, template<typename> class Array>` overloads in backend.hpp
74// (one pair serving real, complex, host, and — via a device extension — device operands).
76// ---- complex products (complex inner-product spaces) ----
77// dot / vdot / inner for complex operands are the SAME two-layer templates in backend.hpp,
78// instantiated at T = std::complex<double>; vdot's conjugation is an `if constexpr` branch. No
79// separate symbols. Complex matmul and conj_transpose are likewise the one generic template each
80// in backend.hpp, instantiating the complex element type — no separate complex symbols.
82/// @cond INTERNAL
83/// the HOST kernels of the factorization/solver seam. Every routine below is the
84/// backend.hpp pattern: an allocating front (inline in this header, `requires NumericArray`) that
85/// validates metadata, allocates via `Array<T>::uninitialized`, and calls the same-named out-param
86/// (or scalar-out) kernel UNQUALIFIED — so these HostArray kernels (defined + instantiated in
87/// routines.cpp) serve host arrays, and a device extension's `requires DeviceArray` overloads are
88/// found by ADL. Declared here, before the fronts, so the fronts' unqualified calls see them.
89/// Each kernel's full documentation sits with its declaration here; the allocating front that
90/// calls it is documented further down this header.
91/**
92 * Matrix power into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
93 * @ref matrix_power (the HOST kernel of the seam pattern; a device extension supplies its own
94 * `requires DeviceArray` overload, found by ADL).
95 * @param out destination; a contiguous n×n matrix, overwritten with Aⁿ.
96 * @param a square matrix.
97 * @param n exponent.
98 * @complexity O(n³·log|n|).
99 * @alloc reuses @p out; the binary-exponentiation products allocate their own scratch.
100 * @test LinalgRoutines.FactorizationOutReusesBuffer
101 */
102template <ndarray::Field T, template <typename> class Array>
103 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
104void matrix_power(Array<T>& out, const Array<T>& a, long long n);
105/**
106 * Cholesky factor into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
107 * @ref cholesky (the HOST kernel of the seam pattern; a device extension supplies its own
108 * `requires DeviceArray` overload, found by ADL).
109 * @param out destination; a contiguous n×n matrix, overwritten with the lower-triangular L.
110 * @param a square SPD matrix.
111 * @complexity O(n³).
112 * @alloc reuses @p out (the factor is computed into private scratch, then copied in).
113 * @test LinalgRoutines.FactorizationOutReusesBuffer
114 */
115template <ndarray::Field T, template <typename> class Array>
116 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
117void cholesky(Array<T>& out, const Array<T>& a);
118/**
119 * Reduced QR into the caller's buffers (outs FIRST) — the buffer-reuse overload of @ref qr,
120 * filling @p q and @p r instead of allocating a @ref QR (the HOST kernel of the seam pattern; a
121 * device extension supplies its own `requires DeviceArray` overload, found by ADL).
122 * @param q destination for the orthonormal factor; a contiguous m×n matrix, overwritten.
123 * @param r destination for the upper-triangular factor; a contiguous n×n matrix, overwritten.
124 * @param a m×n matrix.
125 * @complexity O(m²·n).
126 * @alloc reuses @p q and @p r (the factors are computed into private scratch, then copied in).
127 * @test LinalgRoutines.DecompositionOutReusesBuffer
128 */
129template <ndarray::Field T, template <typename> class Array>
130 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
131void qr(Array<T>& q, Array<T>& r, const Array<T>& a);
132/**
133 * Full SVD into the caller's buffers (outs FIRST) — the buffer-reuse overload of @ref svd,
134 * filling @p u, @p s and @p vh instead of allocating an @ref SVD (the HOST kernel of the seam
135 * pattern; a device extension supplies its own `requires DeviceArray` overload, found by ADL).
136 * @param u destination for the left singular vectors; a contiguous m×n matrix, overwritten.
137 * @param s destination for the singular values; a contiguous length-n vector, overwritten.
138 * @param vh destination for Vᵀ; a contiguous n×n matrix, overwritten.
139 * @param a m×n matrix (rows ≥ cols).
140 * @complexity iterative O(n³).
141 * @alloc reuses @p u, @p s, @p vh (the factors are computed into private scratch, then copied in).
142 * @test LinalgRoutines.DecompositionOutReusesBuffer
143 */
144template <ndarray::Field T, template <typename> class Array>
145 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
146void svd(Array<T>& u, Array<T>& s, Array<T>& vh, const Array<T>& a);
147/**
148 * Singular values into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
149 * @ref svdvals (the HOST kernel of the seam pattern; a device extension supplies its own
150 * `requires DeviceArray` overload, found by ADL).
151 * @param out destination; a contiguous length-min(m,n) vector, overwritten with the descending values.
152 * @param a m×n matrix.
153 * @complexity iterative O(n³).
154 * @alloc reuses @p out; the Golub–Reinsch reduction allocates its own scratch.
155 * @test LinalgRoutines.FactorizationOutReusesBuffer
156 */
157template <ndarray::Field T, template <typename> class Array>
158 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
159void svdvals(Array<T>& out, const Array<T>& a);
160/**
161 * General eigendecomposition into the caller's buffers (outs FIRST) — the buffer-reuse overload of
162 * @ref eig, filling @p values and @p vectors instead of allocating an @ref EigC (the HOST kernel
163 * of the seam pattern; a device extension supplies its own `requires DeviceArray` overload).
164 * @param values destination for the complex eigenvalues; a contiguous length-n vector, overwritten.
165 * @param vectors destination for the complex eigenvectors (columns); a contiguous n×n matrix, overwritten.
166 * @param a square matrix.
167 * @complexity iterative O(n³) for the eigenvalues; general (non-symmetric) eigenvectors
168 * add O(n⁴) (inverse iteration per eigenvalue — see @ref eig).
169 * @alloc reuses @p values and @p vectors (computed into private scratch, then copied in).
170 * @test LinalgRoutines.DecompositionOutReusesBuffer
171 */
172template <ndarray::Field T, template <typename> class Array>
173 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
174void eig(Array<ndarray::complex_of_t<T>>& values, Array<ndarray::complex_of_t<T>>& vectors,
175 const Array<T>& a);
176/**
177 * General eigenvalues into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
178 * @ref eigvals (the HOST kernel of the seam pattern; a device extension supplies its own
179 * `requires DeviceArray` overload, found by ADL).
180 * @param out destination; a contiguous length-n complex vector, overwritten with the descending spectrum.
181 * @param a square matrix.
182 * @complexity iterative O(n³).
183 * @alloc reuses @p out; the Hessenberg + shifted-QR iteration allocates its own scratch.
184 * @test LinalgRoutines.FactorizationOutReusesBuffer
185 */
186template <ndarray::Field T, template <typename> class Array>
187 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
188void eigvals(Array<ndarray::complex_of_t<T>>& out, const Array<T>& a);
189/**
190 * Symmetric/Hermitian eigendecomposition into the caller's buffers (outs FIRST) — the buffer-reuse
191 * overload of @ref eigh, filling @p values and @p vectors instead of allocating a result struct
192 * (the HOST kernel of the seam pattern; a device extension supplies its own `requires DeviceArray`
193 * overload). ONE two-layer kernel serving the real symmetric AND complex Hermitian paths: values
194 * are always the real spectrum, vectors match the input element.
195 * @param values destination for the real eigenvalues; a contiguous length-n vector, overwritten.
196 * @param vectors destination for the eigenvectors (columns); a contiguous n×n matrix, overwritten.
197 * @param a square symmetric (real) / Hermitian (complex) matrix.
198 * @complexity iterative O(n³).
199 * @alloc reuses @p values and @p vectors (computed into private scratch, then copied in).
200 * @test LinalgRoutines.DecompositionOutReusesBuffer
201 */
202template <ndarray::Field T, template <typename> class Array>
203 requires HostArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>
204void eigh(Array<ndarray::real_base_t<T>>& values, Array<T>& vectors, const Array<T>& a);
205/**
206 * Symmetric/Hermitian eigenvalues into the caller's buffer @p out (out FIRST) — the buffer-reuse
207 * overload of @ref eigvalsh (the HOST kernel of the seam pattern; a device extension supplies its
208 * own `requires DeviceArray` overload). One two-layer kernel: the complex Hermitian path is the
209 * same template at T = std::complex<double>, still writing the REAL spectrum.
210 * @param out destination; a contiguous length-n vector, overwritten with the descending eigenvalues.
211 * @param a square symmetric (real) / Hermitian (complex) matrix.
212 * @complexity iterative O(n³).
213 * @alloc reuses @p out; the tridiagonal-QL solver allocates its own scratch.
214 * @test LinalgRoutines.FactorizationOutReusesBuffer
215 */
216template <ndarray::Field T, template <typename> class Array>
217 requires HostArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>
218void eigvalsh(Array<ndarray::real_base_t<T>>& out, const Array<T>& a);
219/**
220 * Solve into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of @ref solve.
221 * @param out destination; a contiguous length-n vector, overwritten with the solution x.
222 * @param a square coefficient matrix.
223 * @param b right-hand-side vector.
224 * @complexity O(n³).
225 * @alloc reuses @p out; the LU factorization allocates its own scratch.
226 * @test LinalgRoutines.FactorizationOutReusesBuffer
227 */
228template <ndarray::Field T, template <typename> class Array>
229 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
230void solve(Array<T>& out, const Array<T>& a, const Array<T>& b);
231/**
232 * Least-squares solution into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
233 * @ref lstsq. Routes through the @ref matmul out-param so the final product is written into @p out
234 * with no allocation.
235 * @param out destination; a contiguous array of the solution's shape, overwritten.
236 * @param a m×n matrix.
237 * @param b right-hand side.
238 * @complexity iterative O(n³).
239 * @alloc reuses @p out; allocates the intermediate n×m pseudo-inverse and its SVD scratch.
240 * @test LinalgRoutines.FactorizationOutReusesBuffer
241 */
242template <ndarray::Field T, template <typename> class Array>
243 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
244void lstsq(Array<T>& out, const Array<T>& a, const Array<T>& b);
245/**
246 * Inverse into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of @ref inv.
247 * @param out destination; a contiguous n×n matrix, overwritten with A⁻¹.
248 * @param a square matrix.
249 * @complexity O(n³).
250 * @alloc reuses @p out; the LU factorization allocates its own scratch.
251 * @test LinalgRoutines.FactorizationOutReusesBuffer
252 */
253template <ndarray::Field T, template <typename> class Array>
254 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
255void inv(Array<T>& out, const Array<T>& a);
256/**
257 * Pseudo-inverse into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
258 * @ref pinv.
259 * @param out destination; a contiguous n×m matrix (for an m×n input), overwritten with the pseudo-inverse.
260 * @param a m×n matrix.
261 * @complexity iterative O(n³).
262 * @alloc reuses @p out; the Golub–Reinsch SVD allocates its own scratch.
263 * @test LinalgRoutines.FactorizationOutReusesBuffer
264 */
265template <ndarray::Field T, template <typename> class Array>
266 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
267void pinv(Array<T>& out, const Array<T>& a);
268template <ndarray::Field T, template <typename> class Array>
269 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
270void det(T& out, const Array<T>& a);
271template <ndarray::Field T, template <typename> class Array>
272 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
273void cond(T& out, const Array<T>& a);
274template <ndarray::Field T, template <typename> class Array>
275 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
276void matrix_rank(long long& out, const Array<T>& a);
277template <ndarray::Field T, template <typename> class Array>
278 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
279void norm(T& out, const Array<T>& a);
280/// @endcond
281/**
282 * Integer matrix power Aⁿ (negative n via @ref inv).
283 *
284 * Requires a square matrix (throws otherwise); n == 0 returns the identity, and
285 * negative n first inverts @p a via @ref inv (so it inherits @ref inv's
286 * singular-matrix behavior) before raising to |n|.
287 * @param a square matrix.
288 * @param n exponent.
289 * @return Aⁿ.
290 * @complexity O(n³·log|n|) by binary exponentiation.
291 * @alloc allocates a new NDArray result; the binary-exponentiation @ref matmul steps
292 * allocate their own intermediates.
293 * @test LinalgRoutines.MatrixPower
294 * @crtest LinalgCompileRun.MatrixPower
295 * @systest StdlibE2E.Linalg
296 */
297template <ndarray::Field T, template <typename> class Array>
298 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
299[[nodiscard]] Array<T> matrix_power(const Array<T>& a, long long n) {
300 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
301 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
302 Array<T> out = Array<T>::uninitialized({a.shape()[0], a.shape()[0]});
303 matrix_power(out, a, n);
304 return out;
306// ---- Decompositions ----
307/**
308 * Cholesky factor of a symmetric positive-definite matrix (throws otherwise).
309 *
310 * Requires a square matrix and computes L column by column reading only the
311 * lower triangle of @p a; if any pivot (the diagonal under the square root) is
312 * non-positive it throws "matrix is not positive-definite", which also catches
313 * non-SPD or non-symmetric input.
314 * @param a square SPD matrix.
315 * @return lower-triangular L with A = L·Lᵀ.
316 * @complexity O(n³).
317 * @alloc allocates a new NDArray result; the factor is computed into O(n²) private
318 * scratch, then copied in.
319 * @test LinalgRoutines.CholeskyAndQR
320 * @crtest LinalgCompileRun.Cholesky
321 * @systest StdlibE2E.Linalg
322 */
323template <ndarray::Field T, template <typename> class Array>
324 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
325[[nodiscard]] Array<T> cholesky(const Array<T>& a) { // lower-triangular L (A = L Lᵀ)
326 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
327 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
328 Array<T> out = Array<T>::uninitialized({a.shape()[0], a.shape()[0]});
329 cholesky(out, a);
330 return out;
332/** Result of qr(): A = q·r with orthonormal q and upper-triangular r. */
333/// The factors of a QR decomposition. Templated over the container type so a host `qr` yields
334/// `QR<NDArray>` and a device `qr` yields `QR<device_array<T>>`; `ArrT` defaults to `NDArray`
335/// (the host double result), so plain `QR` still names the common host type.
336template <class ArrT = NDArray>
337struct QR {
338 ArrT q; ///< Orthonormal columns, m×n (the Q in A = Q·R).
339 ArrT r; ///< Upper-triangular factor, n×n (the R in A = Q·R).
340};
341/**
342 * Reduced QR via Householder reflections (requires rows ≥ cols).
343 *
344 * Applies successive Householder reflectors to triangularize @p a, returning the
345 * thin/reduced factors; throws "qr requires rows >= cols" for wide matrices.
346 * Rank-deficient columns (zero pivot norm) are skipped, leaving the
347 * corresponding R entries zero.
348 * @param a m×n matrix.
349 * @return @ref QR with q (m×n, orthonormal cols) and r (n×n, upper-triangular).
350 * @complexity O(m²·n) — n reflectors each applied to the full m×m Q; O(n³) when square.
351 * @alloc allocates both members; the factorization works in O(m²) private scratch
352 * (a full m×m Q workspace), then copies the reduced factors in.
353 * @test LinalgRoutines.CholeskyAndQR
354 * @crtest LinalgCompileRun.Qr
355 * @systest StdlibE2E.Linalg
356 */
357template <ndarray::Field T, template <typename> class Array>
358 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
359[[nodiscard]] QR<Array<T>> qr(const Array<T>& a) {
360 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
361 const std::size_t m = a.shape()[0], n = a.shape()[1];
362 if (m < n) throw std::runtime_error("linalg: qr requires rows >= cols");
363 QR<Array<T>> out{Array<T>::uninitialized({m, n}), Array<T>::uninitialized({n, n})};
364 qr(out.q, out.r, a);
365 return out;
367/** Result of svd(): A = u·diag(s)·vh. */
368/// The factors of a singular value decomposition, templated over the container type (`ArrT`
369/// defaults to `NDArray`, the host double result, so plain `SVD` names the common host type).
370template <class ArrT = NDArray>
371struct SVD {
372 ArrT u; ///< Left singular vectors, m×n.
373 ArrT s; ///< Singular values in descending order (length n).
374 ArrT vh; ///< Right singular vectors transposed, n×n (the Vᵀ in A = u·diag(s)·Vᵀ).
375};
376/**
377 * Singular value decomposition (Golub–Reinsch; requires rows ≥ cols).
378 *
379 * Reduces @p a to upper-bidiagonal form by Householder reflections, then diagonalizes
380 * it with implicit-shift QR (accumulating U and V), and sorts the singular values
381 * descending — the world-standard dense SVD (what LAPACK's dgesvd reduces to). Throws
382 * "svd requires rows >= cols" for wide matrices (transpose first); singular values come
383 * out non-negative.
384 * @param a m×n matrix.
385 * @return @ref SVD with u (m×n), s (descending singular values), vh (n×n = Vᵀ).
386 * @complexity iterative O(m·n²); O(n³) when square.
387 * @alloc allocates all members; the Golub–Reinsch reduction allocates its own O(m·n) workspace.
388 * @test LinalgRoutines.SvdAndEigh
389 * @crtest LinalgCompileRun.Svd
390 * @systest StdlibE2E.Linalg
391 */
392template <ndarray::Field T, template <typename> class Array>
393 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
394[[nodiscard]] SVD<Array<T>> svd(const Array<T>& a) {
395 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
396 const std::size_t m = a.shape()[0], n = a.shape()[1];
397 if (m < n) throw std::runtime_error("linalg: svd requires rows >= cols (transpose otherwise)");
398 SVD<Array<T>> out{Array<T>::uninitialized({m, n}), Array<T>::uninitialized({n}),
399 Array<T>::uninitialized({n, n})};
400 svd(out.u, out.s, out.vh, a);
401 return out;
403/**
404 * Singular values only (≈ `numpy.linalg.svd(a, compute_uv=False)` / `svdvals`).
405 *
406 * Runs the same Golub–Reinsch reduction as @ref svd but takes the **values-only** fast
407 * path — it never accumulates U or V, and skips the U/V Givens rotations in the QR
408 * sweep — so it beats the full decomposition by a margin that grows with n (see the
409 * linalg benchmark page). Accepts any shape (singular values of `a` and `aᵀ` coincide).
410 * @param a m×n matrix.
411 * @return length-min(m,n) vector of singular values, descending.
412 * @complexity iterative O(n³), but a large constant factor below @ref svd.
413 * @alloc allocates a new NDArray result; the reduction still allocates its O(m·n)
414 * workspace (the values-only path skips the U/V accumulation work and result
415 * copies, not the working buffers).
416 * @test LinalgRoutines.SvdAndEigh
417 * @crtest LinalgCompileRun.Svdvals
418 * @systest StdlibE2E.Linalg
419 */
420template <ndarray::Field T, template <typename> class Array>
421 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
422[[nodiscard]] Array<T> svdvals(const Array<T>& a) {
423 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
424 const std::size_t m = a.shape()[0], n = a.shape()[1];
425 Array<T> out = Array<T>::uninitialized({m < n ? m : n});
426 svdvals(out, a);
427 return out;
430// ---- Matrix eigenvalues ----
431/** Result of eigh(): a real spectrum — column j of vectors is the eigenvector for values[j]. */
432template <class ArrT = NDArray>
433struct Eig {
434 ArrT values; ///< Eigenvalues (length n), real.
435 ArrT vectors; ///< Eigenvectors as columns: column j matches values[j] (empty if not computed).
436};
437/**
438 * Result of the general eig(): a **complex** spectrum, since a real matrix can have
439 * complex conjugate eigenvalue pairs. Column j of vectors is the eigenvector for values[j].
440 */
441template <class ArrT = CNDArray>
442struct EigC {
443 ArrT values; ///< Eigenvalues (length n), complex.
444 ArrT vectors; ///< Eigenvectors as columns: column j matches values[j].
445};
446/// Result of the complex Hermitian eigh(): **real** eigenvalues with **complex** eigenvectors —
447/// so the two members have DIFFERENT container/element types (`ValsT` real, `VecsT` complex).
448/// Defaults `NDArray`/`CNDArray` are the host result, so plain `EighC` names the common host type.
449template <class ValsT = NDArray, class VecsT = CNDArray>
450struct EighC {
451 ValsT values; ///< Eigenvalues (length n), real and descending.
452 VecsT vectors; ///< Eigenvectors as columns: column j matches values[j].
453};
454/// The result type of the unified @ref eigh: for a real element, `Eig<Array<T>>` (real values +
455/// vectors); for a complex element, `EighC<Array<real>, Array<T>>` (real values, complex vectors).
456/// eigh's return type differs per element, so it is expressed here (not `auto`) so the header
457/// declaration knows it without seeing the definition.
458template <ndarray::Field T, template <typename> class Array>
459using eigh_result_t = std::conditional_t<ndarray::is_complex_v<T>,
460 EighC<Array<ndarray::real_base_t<T>>, Array<T>>,
461 Eig<Array<T>>>;
462/**
463 * Eigen-decomposition of a general square matrix (**complex** spectrum and
464 * eigenvectors).
465 *
466 * For a symmetric @p a it delegates to @ref eigh (promoted to complex with zero
467 * imaginary part); otherwise it uses Hessenberg reduction + shifted QR for the
468 * eigenvalues, then **inverse iteration** for each eigenvector. A real matrix with a
469 * complex conjugate pair (e.g. a rotation) yields those complex eigenvalues and
470 * eigenvectors rather than throwing. Throws on a non-square matrix or if the QR
471 * iteration fails to converge.
472 * @param a square matrix.
473 * @return @ref EigC with complex values and matching complex eigenvector columns.
474 * @complexity iterative O(n³) via Hessenberg + shifted QR for the eigenvalues; the
475 * general (non-symmetric) eigenvectors add O(n⁴) — one inverse iteration, each
476 * with its own O(n³) complex LU factorization, per eigenvalue (a symmetric @p a
477 * stays O(n³) via @ref eigh, which accumulates the vectors in the QL sweep).
478 * @alloc allocates both members; plus O(n²) factorization scratch (on the non-symmetric
479 * path, a fresh complex n×n LU per eigenvalue).
480 * @test LinalgRoutines.GeneralEig
481 * @crtest LinalgCompileRun.Eig
482 * @systest StdlibE2E.Linalg
483 */
484template <ndarray::Field T, template <typename> class Array>
485 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
486[[nodiscard]] EigC<Array<ndarray::complex_of_t<T>>> eig(const Array<T>& a) { // general square matrix
487 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
488 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
489 const std::size_t n = a.shape()[0];
490 using C = Array<ndarray::complex_of_t<T>>;
491 EigC<C> out{C::uninitialized({n}), C::uninitialized({n, n})};
492 eig(out.values, out.vectors, a);
493 return out;
495/**
496 * Eigenvalues of a general square matrix (**complex**), descending.
497 *
498 * Routes symmetric input through tridiagonal QL and everything else through
499 * Hessenberg + shifted QR, then sorts the result descending (by real part, then by
500 * imaginary part). A real matrix with a complex conjugate pair yields those complex
501 * eigenvalues rather than throwing. Throws on a non-square matrix or non-convergence
502 * of the QR iteration.
503 * @param a square matrix.
504 * @return length-n complex vector of eigenvalues.
505 * @complexity iterative O(n³).
506 * @alloc allocates a new CNDArray result; the iteration allocates its own O(n²) scratch.
507 * @test LinalgRoutines.SvdAndEigh
508 * @crtest LinalgCompileRun.Eigvals
509 * @systest StdlibE2E.Linalg
510 */
511template <ndarray::Field T, template <typename> class Array>
512 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
513[[nodiscard]] Array<ndarray::complex_of_t<T>> eigvals(const Array<T>& a) {
514 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
515 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
516 Array<ndarray::complex_of_t<T>> out =
517 Array<ndarray::complex_of_t<T>>::uninitialized({a.shape()[0]});
518 eigvals(out, a);
519 return out;
521/**
522 * Eigen-decomposition of a symmetric matrix (Householder tridiagonalization + QL).
523 *
524 * Reduces @p a to tridiagonal form by Householder reflections, then diagonalizes it
525 * with implicit-shift QL, returning real eigenvalues sorted descending with matching
526 * eigenvector columns; it reads the full matrix and assumes symmetry rather than
527 * checking it, so asymmetric input yields meaningless results. Throws on a non-square
528 * matrix (or if the QL iteration fails to converge).
529 * @param a square symmetric matrix.
530 * @return @ref Eig (real input) or @ref EighC (complex Hermitian input): descending real values and matching eigenvector columns.
531 * @complexity iterative O(n³).
532 * @alloc allocates both members; the solver allocates its own O(n²) scratch (a complex
533 * Hermitian input first embeds into a 2n×2n real matrix).
534 * @test LinalgRoutines.SvdAndEigh
535 * @test LinalgRoutines.ComplexHermitianEigh
536 * @crtest LinalgCompileRun.Eigh
537 * @crtest LinalgCompileRun.EighComplex
538 * @systest StdlibE2E.Linalg
539 * @systest StdlibE2E.LinalgComplex
540 */
541template <ndarray::Field T, template <typename> class Array>
542 requires NumericArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>
543[[nodiscard]] eigh_result_t<T, Array> eigh(const Array<T>& a) { // symmetric / Hermitian
544 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
545 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
546 const std::size_t n = a.shape()[0];
547 eigh_result_t<T, Array> out{Array<ndarray::real_base_t<T>>::uninitialized({n}),
548 Array<T>::uninitialized({n, n})};
549 eigh(out.values, out.vectors, a);
550 return out;
552/**
553 * Eigenvalues of a symmetric matrix, descending (tridiagonal QL).
554 *
555 * Same tridiagonalization + QL as @ref eigh but **skips the eigenvector accumulation
556 * entirely**, so it runs well under `eigh`'s cost at large n (see the linalg benchmark page); assumes
557 * (does not verify) symmetry and throws on a non-square matrix.
558 * ONE two-layer template collapsing the former real and complex (Hermitian) overloads: a real
559 * element takes the symmetric path, a complex element the Hermitian path (`if constexpr`). The
560 * spectrum is always REAL, returned as `Array<real_base_t<T>>`.
561 * @tparam T the element type (`double` or `std::complex<double>`).
562 * @tparam Array the container.
563 * @param a square symmetric (real) / Hermitian (complex) matrix.
564 * @return length-n vector of real eigenvalues.
565 * @complexity iterative O(n³).
566 * @alloc allocates a new result; the solver allocates its own O(n²) scratch (a complex
567 * Hermitian input first embeds into a 2n×2n real matrix).
568 * @test LinalgRoutines.EigvalshSymmetric
569 * @test LinalgRoutines.ComplexHermitianEigh
570 * @crtest LinalgCompileRun.Eigvalsh
571 * @crtest LinalgCompileRun.EigvalshComplex
572 * @systest StdlibE2E.Linalg
573 * @systest StdlibE2E.LinalgComplex
574 */
575template <ndarray::Field T, template <typename> class Array>
576 requires NumericArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>
577[[nodiscard]] Array<ndarray::real_base_t<T>> eigvalsh(const Array<T>& a) {
578 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
579 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
580 Array<ndarray::real_base_t<T>> out =
581 Array<ndarray::real_base_t<T>>::uninitialized({a.shape()[0]});
582 eigvalsh(out, a);
583 return out;
585// EighC (real values + complex vectors) and the unified two-layer `eigh` are declared above with
586// the other eig-family structs; the complex Hermitian eigh is that template at T = complex<double>,
587// returning EighC<NDArray, CNDArray> via the if-constexpr Hermitian branch — and its buffer-reuse
588// form is the SAME two-layer eigh kernel above at T = complex<double> (values NDArray&, vectors
589// CNDArray&). Likewise complex eigvalsh is the eigvalsh kernel at T = complex<double>.
590// ---- Norms and other numbers ----
591/**
592 * Norm: L2 for vectors, Frobenius for matrices.
593 *
594 * Dispatches on rank: 1-D (or lower) inputs get the Euclidean L2 norm, 2-D
595 * inputs the Frobenius norm; either way it is the square root of the sum of
596 * squared entries. Two-layer over the element and container like every routine (the scalar-out
597 * kernel `norm(out, a)` is the host/device seam; host `double` is the shipped instantiation).
598 * @param a vector or matrix.
599 * @return √Σ xᵢ².
600 * @complexity O(n) for vectors / O(n²) for matrices.
601 * @alloc none for a contiguous operand (summed in place); a non-contiguous view packs
602 * once. Returns a double.
603 * @test LinalgRoutines.NormAndRank
604 * @crtest LinalgCompileRun.Norm
605 * @systest StdlibE2E.Linalg
606 */
607template <ndarray::Field T, template <typename> class Array>
608 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
609[[nodiscard]] T norm(const Array<T>& a) { // default: Frobenius / L2
610 T out;
611 norm(out, a);
612 return out;
614/**
615 * 2-norm condition number σ_max/σ_min (∞ if singular).
616 *
617 * Takes the ratio of largest to smallest singular value from a Golub–Reinsch SVD
618 * (transposing internally for wide matrices); returns +infinity when the
619 * smallest singular value is exactly zero (singular/rank-deficient).
620 * @param a matrix.
621 * @return condition number.
622 * @complexity iterative O(n³) via SVD.
623 * @alloc allocates scratch O(n²) for the factorization; returns a double.
624 * @test LinalgRoutines.SlogdetAndCond
625 * @crtest LinalgCompileRun.Cond
626 * @systest StdlibE2E.Linalg
627 */
628template <ndarray::Field T, template <typename> class Array>
629 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
630[[nodiscard]] T cond(const Array<T>& a) {
631 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
632 T out;
633 cond(out, a);
634 return out;
636/**
637 * Determinant via LU with partial pivoting.
638 *
639 * Computes the product of the LU pivots times the permutation sign; requires a
640 * square matrix (throws otherwise). A singular matrix yields a determinant of
641 * (or extremely near) zero rather than an error.
642 * @param a square matrix.
643 * @return det(A).
644 * @complexity O(n³).
645 * @alloc allocates scratch O(n²) for the factorization; returns a double.
646 * @test LinalgRoutines.SolveDetInv
647 * @crtest LinalgCompileRun.Det
648 * @systest StdlibE2E.Linalg
649 */
650template <ndarray::Field T, template <typename> class Array>
651 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
652[[nodiscard]] T det(const Array<T>& a) {
653 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
654 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
655 T out;
656 det(out, a);
657 return out;
659/**
660 * Numerical rank from SVD singular-value thresholding.
661 *
662 * Counts singular values above a tolerance scaled by the largest singular value
663 * and the matrix size (the standard numpy-style threshold); accepts any shape,
664 * transposing wide matrices internally.
665 * @param a matrix.
666 * @return rank.
667 * @complexity iterative O(n³) via SVD.
668 * @alloc allocates scratch O(n²) for the factorization.
669 * @test LinalgRoutines.NormAndRank
670 * @crtest LinalgCompileRun.MatrixRank
671 * @systest StdlibE2E.Linalg
672 */
673template <ndarray::Field T, template <typename> class Array>
674 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
675[[nodiscard]] long long matrix_rank(const Array<T>& a) {
676 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
677 long long out = 0;
678 matrix_rank(out, a);
679 return out;
681/** Result of slogdet(): det(A) = sign·exp(logabsdet). */
682struct SLogDet {
683 double sign; ///< Sign of the determinant (−1 or +1; a singular matrix keeps ±1, see logabsdet).
684 double logabsdet; ///< Natural log of |det(A)|, so det(A) = sign·exp(logabsdet).
685};
686/// @cond INTERNAL
687/// the scalar-out HOST kernel of the seam pattern (declared after @ref SLogDet,
688/// which its out-parameter needs; a device extension supplies its own DeviceArray overload).
689template <ndarray::Field T, template <typename> class Array>
690 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
691void slogdet(SLogDet& out, const Array<T>& a);
692/// @endcond
693/**
694 * Sign and log|det| via LU (overflow-safe determinant).
695 *
696 * Sums the logs of the absolute LU pivots (avoiding the over/underflow of a raw
697 * product) and tracks the sign from the pivot signs and permutation parity;
698 * requires a square matrix (throws otherwise). A singular matrix gives a hugely
699 * negative logabsdet rather than −infinity, since a zero pivot is nudged to a
700 * tiny value during factorization.
701 * @param a square matrix.
702 * @return @ref SLogDet.
703 * @complexity O(n³).
704 * @alloc allocates scratch O(n²) for the factorization (the struct members are plain doubles).
705 * @test LinalgRoutines.SlogdetAndCond
706 * @crtest LinalgCompileRun.Slogdet
707 * @systest StdlibE2E.Linalg
708 */
709template <ndarray::Field T, template <typename> class Array>
710 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
711[[nodiscard]] SLogDet slogdet(const Array<T>& a) {
712 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
713 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
714 SLogDet out{};
715 slogdet(out, a);
716 return out;
718// Trace — the allocating front `trace(a)` and the scalar-out kernel `trace(out, a)` are the
719// backend.hpp reduction pattern (host kernel here in routines.cpp; a device extension adds its
720// own DeviceArray overload, found by ADL).
722// ---- Solving equations and inverting matrices ----
723/**
724 * Solve A·x = b via LU with partial pivoting.
725 *
726 * Factorizes @p a once then does forward/back substitution against @p b;
727 * requires @p a square and @p b a vector of matching length (throws otherwise).
728 * A singular @p a does not throw but yields a garbage/overflowing solution
729 * (pivots are nudged off zero rather than detected).
730 * @param a square coefficient matrix.
731 * @param b right-hand-side vector.
732 * @return solution x.
733 * @complexity O(n³).
734 * @alloc allocates a new NDArray result; the LU factorization allocates its own O(n²) scratch.
735 * @test LinalgRoutines.SolveDetInv
736 * @crtest LinalgCompileRun.Solve
737 * @systest StdlibE2E.Linalg
738 */
739template <ndarray::Field T, template <typename> class Array>
740 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
741[[nodiscard]] Array<T> solve(const Array<T>& a, const Array<T>& b) { // A x = b
742 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
743 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
744 if (vector_len(b) != a.shape()[0]) throw std::runtime_error("linalg: solve dimension mismatch");
745 Array<T> out = Array<T>::uninitialized({a.shape()[0]});
746 solve(out, a, b);
747 return out;
749/**
750 * Least-squares solution min‖A·x − b‖ (computed as @ref pinv (a)·b).
751 *
752 * Forms the Moore–Penrose pseudo-inverse via SVD and multiplies it by @p b, so
753 * it handles over- and under-determined systems and returns the minimum-norm
754 * solution for rank-deficient @p a; @p b must be conformable for the
755 * @ref matmul step.
756 * @param a m×n matrix.
757 * @param b right-hand side.
758 * @return minimizing x.
759 * @complexity iterative O(n³) via SVD.
760 * @alloc allocates a new NDArray result; plus the intermediate n×m pseudo-inverse and
761 * its SVD scratch.
762 * @test LinalgRoutines.Lstsq
763 * @crtest LinalgCompileRun.Lstsq
764 * @systest StdlibE2E.Linalg
765 */
766template <ndarray::Field T, template <typename> class Array>
767 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
768[[nodiscard]] Array<T> lstsq(const Array<T>& a, const Array<T>& b) { // least-squares solution
769 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
770 if (b.ndim() != 2) throw std::runtime_error("linalg: matmul expects 2-D matrices");
771 if (a.shape()[0] != b.shape()[0])
772 throw std::runtime_error("linalg: matmul inner dimension mismatch");
773 Array<T> out = Array<T>::uninitialized({a.shape()[1], b.shape()[1]});
774 lstsq(out, a, b);
775 return out;
777/**
778 * Matrix inverse via LU with partial pivoting.
779 *
780 * Factorizes @p a once and back-solves against each identity column; requires a
781 * square matrix (throws otherwise). A singular @p a does not throw but produces
782 * garbage/overflowing entries since zero pivots are nudged rather than detected.
783 * @param a square matrix.
784 * @return A⁻¹.
785 * @complexity O(n³).
786 * @alloc allocates a new NDArray result; the LU factorization and the whole-identity
787 * back-solve allocate their own O(n²) scratch.
788 * @test LinalgRoutines.SolveDetInv
789 * @crtest LinalgCompileRun.Inv
790 * @systest StdlibE2E.Linalg
791 */
792template <ndarray::Field T, template <typename> class Array>
793 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
794[[nodiscard]] Array<T> inv(const Array<T>& a) {
795 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
796 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
797 Array<T> out = Array<T>::uninitialized({a.shape()[0], a.shape()[0]});
798 inv(out, a);
799 return out;
801/**
802 * Moore–Penrose pseudo-inverse via SVD (any shape).
803 *
804 * Computes V·diag(1/σ)·Uᵀ from a Golub–Reinsch SVD, transposing wide matrices
805 * internally so any shape works; singular values at or below a size-scaled
806 * tolerance are dropped (treated as zero) so it stays well-defined for
807 * rank-deficient input.
808 * @param a m×n matrix.
809 * @return n×m pseudo-inverse.
810 * @complexity iterative O(n³) via SVD.
811 * @alloc allocates a new NDArray result; the SVD and the assembly allocate their own
812 * O(m·n) scratch.
813 * @test LinalgRoutines.PinvCondRankOnWideMatrix
814 * @crtest LinalgCompileRun.Pinv
815 * @systest StdlibE2E.Linalg
816 */
817template <ndarray::Field T, template <typename> class Array>
818 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
819[[nodiscard]] Array<T> pinv(const Array<T>& a) { // Moore–Penrose pseudo-inverse
820 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
821 Array<T> out = Array<T>::uninitialized({a.shape()[1], a.shape()[0]});
822 pinv(out, a);
823 return out;
826} // namespace cheatah::linalg