Linear algebra
Matrix products, the inverse, Cholesky, QR, the symmetric eigenproblem, the SVD and a positive-definite solver, all written over einsum.
This module gives you the everyday matrix operations — transpose, products, dot products, the identity —
and the classical decompositions: the inverse, the Cholesky factorisation, QR, the eigenvalues and
eigenvectors of a symmetric matrix, the singular value decomposition, and a solver for symmetric
positive-definite systems. Every one of them is written in Hilbert itself over einsum, triangle and
elementwise operations, and every decomposition is an iteration compiled to a bounded loop, so the whole
module runs as an ONNX graph like any other formula.
Conventions. A matrix is a rank-2 tensor n×p, a vector a rank-1 tensor. Decompositions return one
matrix or vector each: where the textbook result is a pair (Q and R, U, Σ and V) the module offers one
definition per factor. A definition whose result has a different rank from its argument (a determinant,
a diagonal, a vector of eigenvalues) should be called from C# through the Tensor overload, by casting
the argument to Tensor.
Accuracy. Because every loop has a fixed number of turns, these are not drop-in replacements for LAPACK
on hard problems. inverse runs forty Newton–Schulz turns: it is accurate to machine precision for
condition numbers up to about \(10^4\), and visibly wrong (errors of several percent) near
\(5 \times 10^5\). eigenvectors runs forty turns of the unshifted QR algorithm, which is accurate when
neighbouring eigenvalues differ clearly in magnitude and loses digits when two are nearly equal.
solve is conjugate gradient and is correct only for symmetric positive-definite matrices; on any other
matrix it returns a wrong answer without warning. The SVD is taken from the eigenproblem of
\(A^\top A\), which squares the condition number, so the small singular values of a badly conditioned
matrix are less accurate than a dedicated SVD would give. cholesky and determinant require a
symmetric positive-definite matrix and return NaN otherwise.
open Prelude.LinearAlgebra
Shapes and positions
A contraction says what to multiply and sum, never where a value sits, so the module builds positions as ordinary tensors. Summing the upper triangle of a sheet of ones down its columns gives \(1, 2, \ldots, p\); less one, that is the vector of column positions \(0, \ldots, p-1\), and “the k-th column” becomes an equality test against it.
identity
identity(a) · a square matrix n×n → identity matrix n×n
Let \(A \in \mathbb{R}^{n \times n}\) (the argument a). Then, for \(i, j = 1, \dots, n\),
The identity matrix with the same shape as a. Only the shape of a is used, never its values.
Intuition The matrix that does nothing: it maps every vector to itself. Hilbert has no way to write “an n×n matrix” without a tensor to take the size from, so you pass the matrix you are working with and only its shape is read.
Use it when you add a ridge or damping term (\(A + \lambda I\)), form a projection complement (\(I - P\)), start an iteration from “no transformation”, or check a result (mul(a, inverse(a)) should be close to it). The argument should be square.
regularisationlinear-systemsshape-handling
def identity(a) ≔ triangle(unit(a), 0) − triangle(unit(a), 1)
zeros
zeros(a) · a tensor → zero tensor shaped like a
Let \(A\) (the argument a) be a tensor of any shape with finite entries. Then
the tensor of the shape of \(A\) whose every entry is \(0\).
A tensor of zeros shaped like a. In a document where Arithmetic’s laws are in force, a − a is reduced
to a plain zero with no shape; that is harmless where the zero is added to a shaped value (as solve
does), but when you need a shaped zero inside a contraction use zeroed from Kernels.
Intuition An empty accumulator with the right shape: subtracting a tensor from itself keeps its extent and clears its values. It is the starting point of loops that add something in every turn.
Use it when a loop needs an initial sum, factor or estimate — cholesky starts \(L\) from zeros(a) and solve starts \(x\) from zeros(b). Where Arithmetic is open, cancels turns it into an unshaped zero, which is fine when it is only added to; inside a contraction use zeroed from Kernels instead. Entries must be finite, since infinity minus infinity is NaN.
shape-handlinglinear-systems
def zeros(a) ≔ a − a
Building blocks
unit(a)— Let \(A\) be a tensor of any shape. Then \(\operatorname{unit}(A)\) has the shape of \(A\) and every entry equal to \(A_{ij}^0 = 1\). A tensor of ones shaped likea.upper(a)— Let \(A \in \mathbb{R}^{n \times p}\). Then \(\operatorname{upper}(A)_{ij} = [\,j \ge i\,]\) for \(i = 1, \dots, n\) and \(j = 1, \dots, p\). The upper triangle (diagonal included) of a sheet of ones shaped likea.index(a)— Let \(A \in \mathbb{R}^{n \times p}\) with \(n \ge p\). Then \(\operatorname{index}(A)_j = \sum_{i=1}^{n} [\,j \ge i\,] - 1 = j - 1\) for \(j = 1, \dots, p\). The column positions \(0, 1, \dots, p-1\) of a matrix, as a vector.onehot(a, k)— Let \(A \in \mathbb{R}^{n \times p}\) with \(n \ge p\), and let \(k \in \{0, 1, \dots, p-1\}\) be a position counted from zero. Then \(\operatorname{onehot}(A, k)_j = [\,j - 1 = k\,]\), that is \(\operatorname{onehot}(A, k) = e_{k+1} \in \mathbb{R}^p\). A vector with 1 at positionkand 0 elsewhere.from(a, k)— Let \(A \in \mathbb{R}^{n \times p}\) with \(n \ge p\), and let \(k\) be a position counted from zero. Then \(\operatorname{from}(A, k)_j = [\,j - 1 \ge k\,]\) for \(j = 1, \dots, p\). A vector with 1 at every position at or afterkand 0 before it.
Contractions
The basic products, each a single einsum. These are the operations the rest of the module, and most
of the Prelude, is written with.
transpose
transpose(a) · a matrix n×p → matrix p×n
Let \(A \in \mathbb{R}^{n \times p}\) (the argument a). Then
The transpose \(A^\top\).
Intuition Flip the matrix over its diagonal, so rows become columns. For a data matrix it switches between “one row per observation” and “one row per feature”; for a linear map it gives the adjoint, which carries errors backwards through the map.
Use it when you form \(A^\top A\) or \(A^\top b\) (normal equations, gradients of a linear model), change the orientation of data, or undo an orthogonal matrix, whose transpose is its inverse. It is defined for rank-2 tensors only; a vector needs no transpose here, since dot and outer fix the orientation.
linear-systemsregressionshape-handling
def transpose(a) ≔ einsum("ij->ji", a)
mul
mul(a, b) · a matrix n×m · b matrix m×p → matrix n×p
Let \(A \in \mathbb{R}^{n \times m}\) and \(B \in \mathbb{R}^{m \times p}\) (the arguments a and b). Then
The matrix product \(AB\). Note that a · b in Hilbert is the elementwise product; use mul for the
matrix product.
Intuition Doing one linear map after another: \(AB\) first applies \(B\), then \(A\). Each entry is the dot product of a row of a with a column of b, which is why the inner sizes must match.
Use it when you compose transformations (rotations, change of basis, chained layers), form covariance-like products, or check an inverse. Cost is \(O(nmp)\), cubic for square matrices. The most common bug is writing a · b, which multiplies entry by entry; for a matrix times a vector use apply, and for \(A^\top A\) use inner.
linear-systemscoordinate-transforms
def mul(a, b) ≔ einsum("ij,jk->ik", a, b)
apply
apply(a, x) · a matrix n×p · x vector p → vector n
Let \(A \in \mathbb{R}^{n \times p}\) (the argument a) and \(x \in \mathbb{R}^p\). Then
The matrix–vector product \(Ax\).
Intuition Push a vector through a linear map. Read row-wise, each output is a weighted sum of the inputs; read column-wise, the result is a mix of the columns of a with weights x.
Use it when you predict with a linear model (\(\hat y = X\beta\), one prediction per row), step linear dynamics (\(x_{t+1} = A x_t\)), or evaluate a residual \(Ax - b\). It costs \(O(np)\), cheaper than building a column matrix and calling mul. For a map applied from the right (row vector times matrix), use apply(transpose(a), x).
regressionlinear-systemssimulation
def apply(a, x) ≔ einsum("ij,j->i", a, x)
outer
outer(x, y) · x vector n · y vector p → matrix n×p
Let \(x \in \mathbb{R}^n\) and \(y \in \mathbb{R}^p\). Then
\[\operatorname{outer}(x, y)_{ij} = (x y^{\top})_{ij} = x_i y_j .\]The outer product \(xy^\top\).
Intuition Every pairwise product of two vectors laid out as a table. The result is the simplest matrix there is: rank one, mapping any input onto the direction of x, scaled by how much the input lines up with y.
Use it when you build a rank-one update (as cholesky subtracts \(ff^\top\) each turn), a covariance contribution \((x-\mu)(x-\mu)^\top\), a projector \(uu^\top\) for a unit vector, or a table of all pairwise interactions. Memory is n×p, so avoid it for very long vectors when all you need is a scalar — that is dot.
spreadlinear-systemsgeometry
def outer(x, y) ≔ einsum("i,j->ij", x, y)
dot
dot(x, y) · x vector n · y vector n → scalar
Let \(x, y \in \mathbb{R}^n\). Then
\[\operatorname{dot}(x, y) = x^{\top} y = \sum_{i=1}^{n} x_i y_i .\]The inner product of two vectors of the same length.
Intuition How much two vectors point the same way, scaled by their lengths: \(\lVert x \rVert \lVert y \rVert \cos\theta\). It is positive when they agree, zero when they are perpendicular and negative when they oppose. With one argument a weight vector, it is simply a weighted sum.
Use it when you compute a weighted total (portfolio return, a linear score), project onto a direction, test orthogonality, or build a similarity (divide by both lengths for the cosine). Both vectors must have the same length. For all pairs of rows at once, use inner or mul instead of looping.
similarity-searchgeometrydata-summary
def dot(x, y) ≔ einsum("i,i->", x, y)
length
length(x) · x vector n → scalar
Let \(x \in \mathbb{R}^n\). Then
\[\operatorname{length}(x) = \lVert x \rVert = \sqrt{x^{\top} x} = \Big( \sum_{i=1}^{n} x_i^2 \Big)^{1/2} .\]The Euclidean norm of a vector.
Intuition Pythagoras in \(n\) dimensions: the straight-line size of a vector. Because squares are summed, one large component dominates, which is why it is the natural measure of error for Gaussian noise.
Use it when you measure a distance (length(x − y)), normalise to a unit vector (x / length(x), guarding a zero vector), clip a gradient, or test convergence by the size of a residual. It is for vectors; a matrix’s size is better described by its singular values.
geometryfeature-scalingsimilarity-search
def length(x) ≔ sqrt(dot(x, x))
diag
diag(a) · a square matrix n×n → vector n
Let \(A \in \mathbb{R}^{n \times n}\) (the argument a). Then
The diagonal of a, as a vector.
Intuition The entries where row and column agree. In a covariance matrix they are each variable’s own variance; in a Gram matrix, each column’s squared length; in a triangular factor, the pivots whose product is the determinant. Since einsum cannot address a position, the identity serves as a mask.
Use it when you read variances or standard errors off a covariance matrix, take a trace (einsum("i->", diag(a))), read eigenvalues off a diagonalised matrix (as eigenvalues does), or build a diagonal preconditioner (as solve does). The argument must be square.
spreaduncertainty-quantificationlinear-systems
def diag(a) ≔ einsum("ij->j", a · identity(a))
inner
inner(a) · a matrix n×p → matrix p×p
Let \(A \in \mathbb{R}^{n \times p}\) (the argument a). Then
The Gram matrix \(A^\top A\).
Intuition A table of how every column of a lines up with every other column: entry \((j, k)\) is the dot product of columns \(j\) and \(k\). For centred data, dividing by \(n - 1\) gives the covariance matrix. It is symmetric and positive semi-definite by construction.
Use it when you set up least-squares normal equations (\(X^\top X \beta = X^\top y\), then solve), compute a covariance, or take the SVD (as rightVectors does). Cost is \(O(np^2)\). Forming it squares the condition number, so badly scaled columns lose accuracy; scale or centre features first.
regressionspreaddimensionality-reduction
def inner(a) ≔ mul(transpose(a), a)
similar
similar(a, q) · a matrix n×n · q matrix n×n → matrix n×n
Let \(A, Q \in \mathbb{R}^{n \times n}\) (the arguments a and q). Then
The similarity transform \(Q^\top A Q\). When q is orthogonal this leaves the eigenvalues of a
unchanged.
Intuition The same linear map described in a different coordinate system, whose axes are the columns of q. When those axes are orthonormal, nothing intrinsic changes — the eigenvalues stay put — so this is how a matrix is turned toward a simpler, more diagonal view.
Use it when you express a covariance in a rotated basis (in its eigenvector basis it becomes diagonal, which is how eigenvalues reads them off), or change coordinates of a quadratic form. It uses \(Q^\top\), not \(Q^{-1}\), so it is a true similarity only for an orthogonal q such as orthogonal or eigenvectors returns.
coordinate-transformsdimensionality-reduction
def similar(a, q) ≔ mul(transpose(q), mul(a, q))
Building blocks
safe(d)— Let \(d\) be a tensor of any shape. Then \(\operatorname{safe}(d) = d + [\,d = 0\,]\) entrywise. Replaces every zero indby one, so it can be used as a divisor; where it matters the term being divided is zero anyway.away(t)— Let \(t\) be a tensor of any shape. Then \(\operatorname{away}(t) = \operatorname{sign}(t) + [\,\operatorname{sign}(t) = 0\,]\) entrywise, which is \(1\) where \(t \ge 0\) and \(-1\) where \(t < 0\). The sign oft, with zero counted as positive (+1).
The inverse: Newton–Schulz
The inverse is found by the Newton–Schulz iteration, which roughly doubles the number of correct digits every turn once it is close:
\[X_{k+1} = X_k (2I - A X_k), \qquad X_0 = \frac{A^\top}{\lVert A \rVert_1 \, \lVert A \rVert_\infty}\]The starting point guarantees convergence for any non-singular matrix, but a badly conditioned matrix spends many turns before the fast phase begins.
inverse
inverse(a) · a non-singular square matrix n×n → matrix n×n
Let \(A \in \mathbb{R}^{n \times n}\) (the argument a) be non-singular, and let
Let \(X_0 = A^{\top} / s\) and \(X_{k+1} = X_k (2 I_n - A X_k)\) for \(k = 0, 1, \dots\). Then
\[\operatorname{inverse}(A) = X_{40} .\]The inverse \(A^{-1}\) after forty Newton–Schulz turns. The matrix need not be symmetric. Forty turns
give machine precision for condition numbers up to about \(10^4\); for a worse-conditioned matrix the
result can be far from the true inverse, and for a symmetric positive-definite system you should prefer
solve to multiplying by an inverse.
Intuition The map that undoes a. Newton–Schulz finds it the way Newton’s method finds \(1/a\) without dividing: guess, measure how far \(AX\) is from \(I\), correct. The scaled transpose as a first guess is small enough to be safe, and once close the number of correct digits roughly doubles every turn.
Use it when you need the matrix itself — a precision matrix from a covariance, \((X^\top X)^{-1}\) for standard errors — or a system that is not symmetric positive-definite. Each of the forty turns costs two matrix products, \(O(n^3)\). Trust it only for condition numbers up to about \(10^4\). To solve a symmetric positive-definite system, solve is cheaper and more accurate than multiplying by the inverse.
linear-systemsuncertainty-quantificationregression
def inverse(a) ≔ nsLoop(40, a, seed(a))
Building blocks
scale(a)— Let \(A \in \mathbb{R}^{n \times p}\). Then \(\operatorname{scale}(A) = \big(\max_j \sum_i \lvert A_{ij} \rvert\big) \cdot \big(\max_i \sum_j \lvert A_{ij} \rvert\big) = \lVert A \rVert_1 \lVert A \rVert_\infty\). The product of the largest absolute column sum and the largest absolute row sum ofa.seed(a)— Let \(A \in \mathbb{R}^{n \times n}\). Then \(\operatorname{seed}(A) = A^{\top} / \operatorname{scale}(A)\). The starting point of the iteration.newton(a, x)— Let \(A, X \in \mathbb{R}^{n \times n}\). Then \(\operatorname{newton}(A, X) = X (2 I_n - A X)\). One Newton–Schulz turn.nsLoop(n, a, x)— Let \(A, X \in \mathbb{R}^{p \times p}\) and \(n \in \mathbb{N}\). Then \(\operatorname{nsLoop}(0, A, X) = X\) and \(\operatorname{nsLoop}(n, A, X) = \operatorname{nsLoop}(n - 1, A, \operatorname{newton}(A, X))\) for \(n \ge 1\).newtonrepeatedntimes, starting fromx; used byinverse.
Cholesky factorisation
For a symmetric positive-definite matrix, \(A = LL^\top\) with \(L\) lower triangular. The module uses the outer-product form, which needs no pivoting on such a matrix: each turn takes the next column of what remains, divides it by the square root of its diagonal entry, writes it into \(L\), and subtracts its outer product from the remainder.
cholesky
cholesky(a) · a symmetric positive-definite matrix n×n → lower-triangular matrix n×n
Let \(A \in \mathbb{R}^{n \times n}\) (the argument a) be symmetric positive-definite. Let \(A^{(0)} = A\)
and \(L^{(0)} = 0\), and for \(s = 1, \dots, n\) (column \(s\), position \(s - 1\) counted from zero) let
Then
\[\operatorname{cholesky}(A) = L^{(n)} , \qquad A = L^{(n)} L^{(n)\top} .\]The Cholesky factor \(L\), with \(A = LL^\top\). The result is exact up to rounding after n turns. A
matrix that is not positive definite produces NaN entries.
Intuition A “square root” of a positive-definite matrix: \(L\) is to a covariance what a standard deviation is to a variance. Multiplying independent standard noise by \(L\) gives noise with covariance \(A\), and the factorisation only exists when \(A\) really is positive definite.
Use it when you draw correlated Gaussian samples (\(\mu + Lz\), with \(z\) supplied by the caller), test whether a matrix is positive definite (NaN entries mean it is not), or compute a determinant (as determinant does). It runs rows(a) turns, each an n×n outer product, so \(O(n^3)\) in total, and it is exact up to rounding, unlike the fixed-turn iterations.
simulationprobabilistic-modellinglinear-systems
def cholesky(a) ≔ chLoop(rows(a), a, zeros(a))
determinant
determinant(a) · a symmetric positive-definite matrix n×n → scalar
Let \(A \in \mathbb{R}^{n \times n}\) (the argument a) be symmetric positive-definite, and let
\(L = \operatorname{cholesky}(A)\). Then
where the last two equalities hold when every \(L_{ii} \neq 0\).
The determinant, computed from the Cholesky factor as the product of its squared diagonal, summed in logarithms. A diagonal entry of \(L\) that is exactly zero contributes a factor of one rather than zero, so a singular matrix does not come back with determinant zero.
It is valid only for symmetric positive-definite matrices, whose determinant is always positive; for
any other matrix it returns NaN or a wrong value.
Intuition The factor by which the matrix scales volume. For a covariance matrix it measures how spread out the whole cloud is (the “generalised variance”); a value near zero means the variables are nearly redundant. Summing logarithms of the Cholesky diagonal avoids multiplying many small or large numbers directly.
Use it when you need the normalising constant of a multivariate Gaussian or compare the total spread of two covariance estimates. For large matrices the determinant itself can overflow or underflow, and you usually want its logarithm: write 2 · einsum("i->", ln(diag(cholesky(a)))) rather than ln(determinant(a)). It does not detect singularity — a zero pivot counts as one — and it is meaningless for non-positive-definite matrices.
probabilistic-modellingdensity-estimationspread
def determinant(a) ≔ exp(2 · einsum("i->", ln(safe(abs(diag(cholesky(a)))))))
Building blocks
at(a, k)— Let \(A \in \mathbb{R}^{n \times p}\) and \(k \in \{1, \dots, n\}\). Then \(\operatorname{at}(A, k) = n - k\), a position counted from zero (column \(s = n - k + 1\)). The column handled by turnk, counting down:rows(a) − k.taken(a, k)— Let \(A \in \mathbb{R}^{n \times n}\), \(k \in \{1, \dots, n\}\) and \(s = n - k + 1\). Then \(\operatorname{taken}(A, k)_i = A_{is} \, [\,i \ge s\,]\) for \(i = 1, \dots, n\). That column of the remainder, with the entries above the diagonal zeroed.pivot(a, k)— Let \(A \in \mathbb{R}^{n \times n}\), \(k \in \{1, \dots, n\}\) and \(s = n - k + 1\). Then \(\operatorname{pivot}(A, k) = \operatorname{taken}(A, k)^{\top} e_s = A_{ss}\). The diagonal entry of that column.factor(a, k)— Let \(A \in \mathbb{R}^{n \times n}\), \(k \in \{1, \dots, n\}\), \(t = \operatorname{taken}(A, k)\) and \(d = \operatorname{pivot}(A, k)\). Then \(\operatorname{factor}(A, k) = t / \sqrt{d + [\,d = 0\,]}\). The column divided by the square root of its pivot: the next column of L.chLoop(k, a, l)— Let \(A, L \in \mathbb{R}^{n \times n}\) and \(k \in \mathbb{N}\), and write \(f = \operatorname{factor}(A, k)\) and \(s = n - k + 1\). Then \(\operatorname{chLoop}(0, A, L) = L\) and \(\operatorname{chLoop}(k, A, L) = \operatorname{chLoop}(k - 1, A - f f^{\top}, L + f e_s^{\top})\) for \(k \ge 1\). One column per turn: subtracts the factor’s outer square fromaand writes the factor intol.
QR by Householder reflections
\(A = QR\) with \(Q\) orthogonal and \(R\) upper triangular. Each turn builds a reflection \(H = I - 2uu^\top\) that takes one column to a multiple of a unit vector, and applies it as a contraction against \(u\) rather than as a full matrix product. The loop carries \(Q\); \(R\) is then computed exactly as \(Q^\top A\).
orthogonal
orthogonal(a) · a square matrix n×n → orthogonal matrix n×n
Let \(A \in \mathbb{R}^{n \times n}\) (the argument a). Let \(Q_0 = I_n\) and \(R_0 = A\), and for
\(s = 1, \dots, n\) (column \(s\), position \(s - 1\) counted from zero) let
where \(\operatorname{away}(t) = 1\) for \(t \ge 0\) and \(-1\) otherwise, and \(\operatorname{safe}(d) = d + [\,d = 0\,]\). Then
\[\operatorname{orthogonal}(A) = Q_n = H_1 H_2 \cdots H_n .\]The factor \(Q\) of the QR decomposition of a.
Intuition A set of perpendicular unit axes that span the same columns as a, found by reflecting each column in turn onto a coordinate axis. Mirrors preserve lengths and angles exactly, which makes this far more stable than subtracting projections one column at a time (Gram–Schmidt).
Use it when you need an orthonormal basis for a set of directions, a well-conditioned change of coordinates, or a random rotation (QR of a matrix of Gaussian noise that you supply). It is also the step inside eigenvectors. It takes square matrices, runs rows(a) turns and costs \(O(n^3)\). Column signs are whatever the reflections choose.
coordinate-transformsnumerical-stabilitygeometry
def orthogonal(a) ≔ qrLoop(rows(a), identity(a), a)
triangular
triangular(a) · a square matrix n×n → upper-triangular matrix n×n
Let \(A \in \mathbb{R}^{n \times n}\) (the argument a) and \(Q = \operatorname{orthogonal}(A)\). Then
The factor \(R = Q^\top A\) of the QR decomposition. The entries below the diagonal are zero up to rounding. Diagonal entries may be negative: the decomposition is unique only up to the signs of the columns of \(Q\) and rows of \(R\).
Intuition What is left of a once its orthogonal part is removed: row \(j\) of \(R\) says how much of each column lies along the first \(j\) axes. A diagonal entry near zero means that column adds almost no new direction beyond the columns before it.
Use it when you want to find nearly collinear columns (small \(\lvert R_{jj} \rvert\)), or the absolute determinant of any square matrix, symmetric or not (\(\prod_j \lvert R_{jj} \rvert\)). Signs on the diagonal may be negative. It recomputes orthogonal(a), so asking for both factors runs the QR loop twice.
linear-systemsmodel-selectionnumerical-stability
def triangular(a) ≔ mul(transpose(orthogonal(a)), a)
Building blocks
column(r, k)— Let \(R \in \mathbb{R}^{n \times n}\) and let \(k \in \{0, \dots, n-1\}\) be a position counted from zero. Then \(\operatorname{column}(R, k)_i = R_{i,k+1} \, [\,i \ge k + 1\,]\) for \(i = 1, \dots, n\). Columnkofrwith the entries above rowkzeroed.alpha(r, k)— Let \(R \in \mathbb{R}^{n \times n}\), \(k\) a position counted from zero, and \(c = \operatorname{column}(R, k)\). Then \(\operatorname{alpha}(R, k) = 0 - \operatorname{away}(c_{k+1}) \, \lVert c \rVert = -\operatorname{away}(R_{k+1,k+1}) \, \lVert c \rVert\). The target of the reflection: minus the sign of the diagonal entry (zero counted as positive) times the column’s length.mirror(r, k)— Let \(R \in \mathbb{R}^{n \times n}\) and \(k\) a position counted from zero. Then \(\operatorname{mirror}(R, k) = \operatorname{column}(R, k) - \operatorname{alpha}(R, k) \, e_{k+1}\). The unnormalised Householder vector: the column lessalphaat positionk.reflector(r, k)— Let \(R \in \mathbb{R}^{n \times n}\), \(k\) a position counted from zero, and \(m = \operatorname{mirror}(R, k)\). Then \(\operatorname{reflector}(R, k) = m / (\lVert m \rVert + [\,\lVert m \rVert = 0\,])\). The Householder vector scaled to unit length (left as zero when it is zero).reflected(r, k)— Let \(R \in \mathbb{R}^{n \times n}\), \(k\) a position counted from zero, and \(u = \operatorname{reflector}(R, k)\). Then \(\operatorname{reflected}(R, k) = R - 2 u (R^{\top} u)^{\top} = (I_n - 2 u u^{\top}) R\).rwith the reflection applied from the left.turned(q, r, k)— Let \(Q, R \in \mathbb{R}^{n \times n}\), \(k\) a position counted from zero, and \(u = \operatorname{reflector}(R, k)\). Then \(\operatorname{turned}(Q, R, k) = Q - 2 (Q u) u^{\top} = Q (I_n - 2 u u^{\top})\).qwith the reflection applied from the right.qrLoop(k, q, r)— Let \(Q, R \in \mathbb{R}^{n \times n}\), \(k \in \mathbb{N}\) and \(c = \operatorname{at}(R, k) = n - k\). Then \(\operatorname{qrLoop}(0, Q, R) = Q\) and \(\operatorname{qrLoop}(k, Q, R) = \operatorname{qrLoop}(k - 1, \operatorname{turned}(Q, R, c), \operatorname{reflected}(R, c))\) for \(k \ge 1\). One column per turn, carryingqandr; returnsq.
The symmetric eigenproblem
The eigenvalues and eigenvectors of a symmetric matrix are found by the QR algorithm: repeatedly replacing \(A\) with \(Q^\top A Q\), an orthogonal similarity that keeps the eigenvalues and drives the off-diagonal entries to zero. Each turn shrinks an off-diagonal entry by roughly \(\lvert \lambda_{j+1} / \lambda_j \rvert\), so a spectrum whose neighbouring eigenvalues differ by a factor of two converges to machine precision well within the forty turns used. When two eigenvalues are nearly equal in magnitude, convergence is slow and those two come back less accurately than a shifted LAPACK iteration would give them.
For a general symmetric matrix the eigenvalues emerge in descending order of magnitude, the order principal components are wanted in. A matrix that is already diagonal is left as it is, so its eigenvalues come back in their original order.
eigenvectors
eigenvectors(a) · a symmetric matrix n×n → orthogonal matrix n×n
Let \(A \in \mathbb{R}^{n \times n}\) (the argument a) be symmetric. Let \(A_0 = A\) and \(V_0 = I_n\), and
for \(t = 0, 1, \dots\) let
Then
\[\operatorname{eigenvectors}(A) = V_{40} = Q_0 Q_1 \cdots Q_{39} .\]The eigenvectors of a, one per column, in the same order as eigenvalues(a). Each vector has unit
length; its sign is arbitrary.
Intuition The directions a symmetric matrix only stretches and never turns. For a covariance matrix they are the axes of the data’s ellipsoid — the principal components — mutually perpendicular. Repeated QR similarities rotate the matrix until it is diagonal, and the accumulated rotations are these axes.
Use it when you do PCA on a covariance or correlation matrix, find principal axes of a shape, or compute spectral embeddings from a symmetric similarity matrix. Each of the forty turns is a full QR, so cost is about \(40 \cdot O(n^3)\). Columns come out by descending magnitude of eigenvalue (a large negative eigenvalue comes early), signs are arbitrary, and nearly equal eigenvalues are resolved poorly. For a non-symmetric matrix the result is not meaningful.
dimensionality-reductionclusteringgeometry
def eigenvectors(a) ≔ eigLoop(40, a, identity(a))
eigenvalues
eigenvalues(a) · a symmetric matrix n×n → vector n
Let \(A \in \mathbb{R}^{n \times n}\) (the argument a) be symmetric and \(V = \operatorname{eigenvectors}(A)\). Then
where \(v_j\) is column \(j\) of \(V\).
The eigenvalues of a, read from the diagonal of \(V^\top A V\) where \(V\) is eigenvectors(a).
Intuition How much the matrix stretches along each of its eigenvector axes. For a covariance matrix, each one is the variance of the data along a principal component, and together they add up to the total variance.
Use it when you choose how many principal components to keep (share of variance explained), check positive definiteness (all eigenvalues above zero), or estimate conditioning (largest over smallest magnitude). It runs eigenvectors internally, with the same cost and the same order by descending magnitude; the loss of accuracy for nearly equal eigenvalues shows up here too.
dimensionality-reductionspreadnumerical-stability
def eigenvalues(a) ≔ diag(similar(a, eigenvectors(a)))
Building blocks
rotate(a, v)— Let \(A, V \in \mathbb{R}^{n \times n}\). Then \(\operatorname{rotate}(A, V) = V \operatorname{orthogonal}(A)\). Accumulates one turn’s orthogonal factor into the eigenvector matrixv.eigLoop(k, a, v)— Let \(A, V \in \mathbb{R}^{n \times n}\), \(k \in \mathbb{N}\) and \(Q = \operatorname{orthogonal}(A)\). Then \(\operatorname{eigLoop}(0, A, V) = V\) and \(\operatorname{eigLoop}(k, A, V) = \operatorname{eigLoop}(k - 1, Q^{\top} A Q, V Q)\) for \(k \ge 1\).kturns of the QR algorithm, carrying the transformed matrix and the accumulated vectors; returns the vectors.
The singular value decomposition
\(A = U \Sigma V^\top\). Here \(V\) and \(\Sigma\) come from the eigenproblem of \(A^\top A\), and
\(U = AV\Sigma^{-1}\). The matrix a may be rectangular, n×p. Forming \(A^\top A\) squares the
condition number, so for a badly conditioned matrix the smallest singular values are less accurate than
a one-sided Jacobi or LAPACK SVD would give. The three factors follow the order of eigenvalues
described above.
rightVectors
rightVectors(a) · a matrix n×p → orthogonal matrix p×p
Let \(A \in \mathbb{R}^{n \times p}\) (the argument a). Then
The right singular vectors \(V\), one per column.
Intuition For a data matrix with one row per observation, the orthogonal directions in feature space ordered by how much of the data’s energy lies along them. Centre the columns first and they are exactly the principal directions of PCA.
Use it when you work from a rectangular data matrix rather than a precomputed covariance: PCA loadings, low-rank approximation, latent factors in a ratings table. It builds \(A^\top A\), a p×p matrix, so it is cheap when there are few columns, but the squared condition number costs accuracy on the weakest directions. Given a covariance matrix already, use eigenvectors.
dimensionality-reductionfeature-extraction
def rightVectors(a) ≔ eigenvectors(inner(a))
singular
singular(a) · a matrix n×p → vector p
Let \(A \in \mathbb{R}^{n \times p}\) (the argument a) and \(\lambda = \operatorname{eigenvalues}(A^{\top} A) \in \mathbb{R}^p\). Then
The singular values, the square roots of the eigenvalues of \(A^\top A\). Eigenvalues that rounding makes slightly negative are floored at zero.
Intuition The lengths of the axes of the ellipse that a makes from the unit sphere. The largest is the most the matrix can stretch any vector; one near zero means some direction is almost flattened out, so the matrix has close to lower rank.
Use it when you estimate the numerical rank (count values above a tolerance), the condition number (largest over smallest), the spectral norm (the largest), or how much energy a rank-\(k\) approximation keeps. Values are non-negative and follow the order of eigenvalues. Small singular values of a badly conditioned matrix are the least accurate, because \(A^\top A\) is formed first.
dimensionality-reductionnumerical-stabilitymodel-selection
def singular(a) ≔ sqrt(max(eigenvalues(inner(a)), 0))
leftVectors
leftVectors(a) · a matrix n×p → matrix n×p
Let \(A \in \mathbb{R}^{n \times p}\) (the argument a), \(V = \operatorname{rightVectors}(A)\) and
\(\sigma = \operatorname{singular}(A)\). Then
The left singular vectors \(U = AV\Sigma^{-1}\), one per column (the thin form, n×p). A column whose singular value is zero comes back as zeros rather than a division by zero.
Intuition Where each observation lands along each principal direction, rescaled so every column has unit length. Multiply back by the singular values and you have the PCA scores \(AV\): the coordinates of each row in the new axes.
Use it when you project observations into a low-dimensional space for plotting or clustering, or rebuild a rank-\(k\) approximation \(U_k \Sigma_k V_k^\top\). A zero singular value yields a zero column, not an error. It computes both rightVectors and singular, repeating the eigen-iteration, so when you only need scores, mul(a, rightVectors(a)) is cheaper.
dimensionality-reductionclusteringfeature-extraction
def leftVectors(a) ≔ mul(a, rightVectors(a)) / safe(singular(a))
Solving a positive-definite system
solve finds \(x\) with \(Ax = b\) by conjugate gradient. In exact arithmetic conjugate gradient reaches
the solution of a p×p symmetric positive-definite system in at most p steps, so solve runs exactly
rows(a) turns.
In floating point a badly scaled matrix can stop it far short of the answer, so the system is preconditioned by its own diagonal (Jacobi preconditioning): each residual is divided by \(\operatorname{diag}(A)\). For a Gram matrix \(X^\top X\) this amounts to working with the correlation matrix instead of the covariance matrix, which removes the effect of columns measured on very different scales. The answer is unchanged in exact arithmetic, and a well-scaled matrix pays almost nothing for it.
solve
solve(a, b) · a symmetric positive-definite matrix p×p · b right-hand side vector p → vector p
Let \(A \in \mathbb{R}^{p \times p}\) (the argument a) be symmetric positive-definite and \(b \in \mathbb{R}^p\).
Let \(M = \operatorname{Diag}(m)\) with \(m_j = A_{jj} + [\,A_{jj} = 0\,]\), and write
\(\operatorname{safe}(d) = d + [\,d = 0\,]\). Let \(x_0 = 0\), \(r_0 = b\) and \(p_0 = M^{-1} b\), and for
\(k = 0, 1, \dots, p - 1\) let
Then
\[\operatorname{solve}(A, b) = x_p ,\]which in exact arithmetic is \(A^{-1} b\).
The solution \(x = A^{-1} b\) by Jacobi-preconditioned conjugate gradient, run for exactly \(p\) turns. A denominator that is exactly zero (a residual already at zero) is replaced by one, so the loop keeps the estimate it has instead of dividing by zero.
The matrix must be symmetric positive-definite. On a non-symmetric or indefinite matrix the loop still runs and returns a vector, but it is not the solution.
Intuition Solving \(Ax = b\) for a positive-definite \(A\) is the same as walking downhill on the bowl \(\tfrac12 x^\top A x - b^\top x\). Conjugate gradient picks each new direction so it never undoes progress along the earlier ones, which is why \(p\) steps suffice in exact arithmetic; dividing by the diagonal evens out columns on very different scales.
Use it when you fit least squares from the normal equations (\(X^\top X \beta = X^\top y\)), ridge regression (\(X^\top X + \lambda I\)), kernel ridge or Gaussian-process weights, or a Newton step with a positive-definite Hessian. It runs \(p\) turns of matrix–vector products, \(O(p^3)\) overall. It is silently wrong on non-symmetric or indefinite matrices; use inverse there.
linear-systemsregressionregularisation
def solve(a, b) ≔ cgLoop(rows(a), a, zeros(b), b, eased(a, b))
Building blocks
In the items below \(A \in \mathbb{R}^{p \times p}\), \(x, r, q \in \mathbb{R}^p\) (the search direction p is written \(q\) here, to keep \(p\) for the order), and \(M = \operatorname{Diag}(\operatorname{pivots}(A))\).
pivots(a)— Let \(A \in \mathbb{R}^{p \times p}\). Then \(\operatorname{pivots}(A)_j = A_{jj} + [\,A_{jj} = 0\,]\) for \(j = 1, \dots, p\). The diagonal ofa, with zeros replaced by ones: the preconditioner.eased(a, r)— Let \(A \in \mathbb{R}^{p \times p}\) and \(r \in \mathbb{R}^p\). Then \(\operatorname{eased}(A, r)_j = r_j / \operatorname{pivots}(A)_j\), that is \(\operatorname{eased}(A, r) = M^{-1} r\). A residual divided elementwise by the preconditioner.step(a, r, p)— Let \(A\), \(r\) and \(q\) be as above. Then \(\operatorname{step}(A, r, q) = \dfrac{r^{\top} M^{-1} r}{\operatorname{safe}(q^{\top} A q)}\) with \(\operatorname{safe}(d) = d + [\,d = 0\,]\). The step length along the search directionp.moved(a, x, r, p)— Let \(A\), \(x\), \(r\) and \(q\) be as above. Then \(\operatorname{moved}(A, x, r, q) = x + \operatorname{step}(A, r, q) \, q\). The next estimate of the solution.left(a, r, p)— Let \(A\), \(r\) and \(q\) be as above. Then \(\operatorname{left}(A, r, q) = r - \operatorname{step}(A, r, q) \, A q\). The next residual.gain(a, r, p)— Let \(A\), \(r\) and \(q\) be as above, and \(r' = \operatorname{left}(A, r, q)\). Then \(\operatorname{gain}(A, r, q) = \dfrac{r'^{\top} M^{-1} r'}{\operatorname{safe}(r^{\top} M^{-1} r)}\). The ratio that mixes the old search direction into the new one.next(a, r, p)— Let \(A\), \(r\) and \(q\) be as above, and \(r' = \operatorname{left}(A, r, q)\). Then \(\operatorname{next}(A, r, q) = M^{-1} r' + \operatorname{gain}(A, r, q) \, q\). The next search direction.cgLoop(k, a, x, r, p)— Let \(A\), \(x\), \(r\) and \(q\) be as above and \(k \in \mathbb{N}\). Then \(\operatorname{cgLoop}(0, A, x, r, q) = x\) and \(\operatorname{cgLoop}(k, A, x, r, q) = \operatorname{cgLoop}\big(k - 1, A, \operatorname{moved}(A, x, r, q), \operatorname{left}(A, r, q), \operatorname{next}(A, r, q)\big)\) for \(k \ge 1\).kconjugate-gradient turns carrying estimate, residual and direction; returns the estimate.
Examples
Matrices are printed row by row.
Invert a matrix
open Prelude.LinearAlgebra
fn inverted(a) ↦ inverse(a)
fn product(a, b) ↦ mul(a, b)
using Examples;
var m = new double[,] { { 4, 1 }, { 1, 3 } };
var inv = Algebra.Inverted(m);
var check = Algebra.Product(m, inv);
Output (verified, rounded to 6 significant digits):
inv = [[0.272727, -0.090909], [-0.090909, 0.363636]]
check = [[1, 0], [0, 1]]
The exact inverse is \(\tfrac{1}{11}\begin{pmatrix} 3 & -1 \\ -1 & 4 \end{pmatrix}\). The same doorway
inverts a non-symmetric matrix: for { { 1, 2 }, { 3, 4 } } it returns [[-2, 1], [1.5, -0.5]].
Cholesky factor and determinant
open Prelude.LinearAlgebra
fn lowerFactor(a) ↦ cholesky(a)
fn volume(a) ↦ determinant(a)
using Examples;
using MLambda.Hilbert.Runtime;
var s = new double[,] { { 4, 2, 0 }, { 2, 5, 3 }, { 0, 3, 10 } };
var l = Algebra.LowerFactor(s);
var det = Algebra.Volume((Tensor)s);
Output (verified, rounded to 6 significant digits):
l = [[2, 0, 0], [1, 2, 0], [0, 1.5, 2.78388]]
det = 124
Multiplying l by its transpose gives back s, and \(124 = (2 \cdot 2 \cdot 2.78388)^2\).
Solve a positive-definite system
open Prelude.LinearAlgebra
fn solved(a, b) ↦ solve(a, b)
using Examples;
using MLambda.Hilbert.Runtime;
var s = new double[,] { { 4, 2, 0 }, { 2, 5, 3 }, { 0, 3, 10 } };
var x = Algebra.Solved((Tensor)s, (Tensor)new[] { 6d, 10d, 13d });
Output (verified):
[1, 1, 1]
Each entry of the right-hand side is the corresponding row sum of s, so the answer is a vector of
ones. The requirement of symmetry is real: for the non-symmetric matrix { { 1, 2 }, { 3, 4 } } and
right-hand side { 5, 11 }, whose solution is [1, 2], solve returns [0.707661, 1.83921]. Use
inverse for such a matrix.
Eigenvalues and eigenvectors
open Prelude.LinearAlgebra
fn spectrum(a) ↦ eigenvalues(a)
fn axes(a) ↦ eigenvectors(a)
using Examples;
using MLambda.Hilbert.Runtime;
var a = new double[,] { { 2, 1 }, { 1, 2 } };
var values = Algebra.Spectrum((Tensor)a);
var vectors = Algebra.Axes(a);
Output (verified, rounded to 6 significant digits):
values = [3, 1]
vectors = [[0.707107, -0.707107], [0.707107, 0.707107]]
The first column, \((1, 1)/\sqrt{2}\), belongs to the eigenvalue 3 and the second to 1. With nearly
equal eigenvalues the fixed turn count shows: for { { 1, 0.001 }, { 0.001, 1 } }, whose eigenvalues
are 1.001 and 0.999, spectrum returns [1.00008, 0.99992].
QR decomposition
open Prelude.LinearAlgebra
fn qFactor(a) ↦ orthogonal(a)
fn rFactor(a) ↦ triangular(a)
using Examples;
var r = new double[,] { { 3, 0 }, { 4, 5 } };
var q = Algebra.QFactor(r);
var upper = Algebra.RFactor(r);
Output (verified):
q = [[-0.6, 0.8], [-0.8, -0.6]]
upper = [[-5, -4], [-0, -3]]
q · upper reproduces the input; the signs are those Householder reflections choose.
Singular value decomposition
open Prelude.LinearAlgebra
fn singularValues(a) ↦ singular(a)
fn leftSingular(a) ↦ leftVectors(a)
fn rightSingular(a) ↦ rightVectors(a)
using Examples;
using MLambda.Hilbert.Runtime;
var a = new double[,] { { 3, 0 }, { 4, 5 } };
var sigma = Algebra.SingularValues((Tensor)a);
var u = Algebra.LeftSingular(a);
var v = Algebra.RightSingular(a);
Output (verified, rounded to 6 significant digits):
sigma = [6.7082, 2.23607]
u = [[0.316228, -0.948683], [0.948683, 0.316228]]
v = [[0.707107, -0.707107], [0.707107, 0.707107]]
The singular values are \(\sqrt{45}\) and \(\sqrt{5}\), the square roots of the eigenvalues of \(A^\top A = \begin{pmatrix} 25 & 20 \\ 20 & 25 \end{pmatrix}\).