Matrix factorizations (spec layer)

This file provides real, shape-indexed reference implementations of the two exact, finite matrix factorizations that classical / scientific ML models (Gaussian processes, kernel ridge regression, PCA, least squares) depend on, and which were previously missing from the spec layer:

• choleskySpec — Cholesky factorization $A=LL^\mathsf{T}$ (lower-triangular $L$), proved for matrices with positive executable Cholesky pivots.
• qrSpec — QR factorization $A=QR$ via classical Gram–Schmidt; under positive executable $R$ pivots, $Q$ has orthonormal columns and $R$ is upper-triangular.

It also provides the linear solves that ride on the Cholesky factor:

• triSolveLowerFn / triSolveUpperFn — forward / back triangular substitution;
• cholSolveFn — solve $Ax=b$ from a Cholesky factor of $A$;
• solveRidgeSpec — the Tikhonov / kernel-ridge solve $(K+\gamma I)x=b$.

Verification scope

The verified contribution is the factorizations: choleskySpec / qrSpec come with reconstruction and structural theorems (IsCholesky / IsQR, lower- and upper-triangularity, orthonormality) in NN.Proofs.Tensor.Basic.Factorizations*, under their stated positive-pivot success hypotheses. The triangular- and ridge-solve helpers above (triSolveLowerFn, triSolveUpperFn, cholSolveFn, solveRidgeSpec) are executable APIs only: this PR does not yet prove their correctness (no triangular-solve or solveRidge correctness theorem has landed). They follow the standard substitution formulas over the readable function representation and are exercised by #eval examples, but should not be read as carrying a verified-correctness guarantee.

Intent / tradeoffs

Like the rest of the spec layer (determinantSpec, inverseSpec?, matMulSpec), these prioritize mathematical clarity and shape safety over performance, and are intended for small/medium matrices and proof-oriented reference code. For large-scale numerics, use array-backed runtime kernels.

Internally the algorithms are written over the plain function representation Fin n → Fin n → α (matrices) and Fin n → α (vectors), then wrapped back into Spec.Tensor at the boundary. This keeps the numerical formulas readable and keeps later correctness proofs working on ordinary functions rather than on nested Tensor matches.

Boundary conversions between Spec.Tensor and plain functions

def Spec.toMatFn {α : Type} {m n : ℕ} (A : Tensor α (Shape.dim m (Shape.dim n Shape.scalar))) : Fin m → Fin n → α

View a matrix tensor as a function Fin m → Fin n → α.

def Spec.ofMatFn {α : Type} {m n : ℕ} (f : Fin m → Fin n → α) : Tensor α (Shape.dim m (Shape.dim n Shape.scalar))

Build a matrix tensor from a function Fin m → Fin n → α.

def Spec.toVecFn {α : Type} {n : ℕ} (v : Tensor α (Shape.dim n Shape.scalar)) : Fin n → α

View a vector tensor as a function Fin n → α.

def Spec.ofVecFn {α : Type} {n : ℕ} (f : Fin n → α) : Tensor α (Shape.dim n Shape.scalar)

Build a vector tensor from a function Fin n → α.

Small numeric helpers on the function representation

def Spec.dotFn {α : Type} [Context α] {p : ℕ} (u v : Fin p → α) : α

Dot product of two length-p vectors.

def Spec.normFn {α : Type} [Context α] {p : ℕ} (v : Fin p → α) : α

Euclidean norm of a length-p vector.

Cholesky factorization

For an input whose executable Cholesky pivots are positive, compute the lower-triangular L with $A=LL^\mathsf{T}$. Symmetric positive-definiteness is the standard sufficient condition, but the theorem in this file family is stated against the executable positive-pivot success condition.

The columns are computed left to right. Column j uses only columns 0 .. j-1:

• diagonal: $L_{jj}=\sqrt{A_{jj}-\sum_{k<j}L_{jk}^2}$
• below: $L_{ij}=(A_{ij}-\sum_{k<j}L_{ik}L_{jk})/L_{jj}$ for $i>j$
• above: $L_{ij}=0$ for $i<j$

Trust boundary: the @[implemented_by] performance hooks

Several defs here (choleskyColsFn, cholSolveFn, solveRidgeFn) carry an @[implemented_by …Impl] attribute. The clean closure form is what the correctness proofs reason about; the …Impl companion is a strict, array-backed rewrite that the compiler runs instead, so #eval stays fast (the closure form re-evaluates prefixes exponentially in the interpreter).

This substitution is a trusted runtime boundary. Compiled #eval/runtime code executes the …Impl body while the proofs constrain the clean closure body. The two transcribe the same recurrence, and the numeric examples in NN/Examples/Factorization exercise the compiled path, but a future equivalence theorem should discharge this boundary explicitly. Anything proved about choleskyFn/solveRidgeFn therefore transfers to #eval output only modulo this unverified hook.

def Spec.choleskyColsImpl {α : Type} [Context α] {n : ℕ} (A : Fin n → Fin n → α) : List (Fin n → α)

Strict, array-backed runtime implementation of choleskyColsFn (registered via @[implemented_by]). Each column is materialized into an Array α, so a back-reference L[i,k] is an O(1) lookup rather than a closure that re-evaluates the whole prefix. The closure form below is mathematically clean (and is what the proofs reason about), but reading the full factor L from it re-evaluates columns exponentially — ruinous in the interpreter (#eval). It is intended to compute the same factor strictly; this equivalence is trusted, not proved (see the trust-boundary note above), with the numeric examples ($A=LL^\mathsf{T}$ and ridge-solve residual $\approx0$) as evidence rather than a proof.

@[implemented_by Spec.choleskyColsImpl]

def Spec.choleskyColsFn {α : Type} [Context α] {n : ℕ} (A : Fin n → Fin n → α) : List (Fin n → α)

The list of columns of the Cholesky factor L, as length-n vectors, computed left to right. Element j of the result is column j of L. Built by a left fold so that when column j is formed, cols already holds columns 0 .. j-1.

The runtime implementation is choleskyColsImpl (strict arrays); the closure form here is the one the correctness proofs reason about. The two are intended to compute the same factor — trusted, not proved; see the trust-boundary note above.

def Spec.choleskyFn {α : Type} [Context α] {n : ℕ} (A : Fin n → Fin n → α) : Fin n → Fin n → α

Cholesky factor as a function: $L_{ij}=(\mathtt{choleskyColsFn}(A))_j(i)$.

def Spec.choleskySpec {α : Type} [Context α] {n : ℕ} (A : Tensor α (Shape.dim n (Shape.dim n Shape.scalar))) : Tensor α (Shape.dim n (Shape.dim n Shape.scalar))

Cholesky factorization candidate for A, returning a lower-triangular factor. Over ℝ, the proved reconstruction theorem assumes symmetry and positive executable Cholesky pivots.

PyTorch analogue: torch.linalg.cholesky(A).

Triangular solves and the kernel-ridge (Tikhonov) linear solve

Once $A$ is factored as $A=LL^\mathsf{T}$ (Cholesky), the linear system $Ax=b$ is solved by two triangular substitutions: forward-solve $Lz=b$, then back-solve $L^\mathsf{T}x=z$. Each substitution visits the unknowns in an order such that, when row i is reached, every unknown it depends on has already been computed; the accumulator acc holds those values and 0 everywhere else, so the dot dotFn (row i) acc is exactly the required partial sum (the not-yet-solved and structurally-zero terms drop out).

def Spec.triSolveLowerFn {α : Type} [Context α] {n : ℕ} (L : Fin n → Fin n → α) (b : Fin n → α) : Fin n → α

Forward substitution: solve $Ly=b$ for a lower-triangular $L$ with nonzero diagonal. Unknowns are visited $0,1,\ldots,n-1$; when row $i$ is reached, acc holds $y_0,\ldots,y_{i-1}$ (and $0$ elsewhere), so $\mathtt{dotFn}(L_i,\mathtt{acc})=\sum_{k<i}L_{ik}y_k$ by lower-triangularity.

def Spec.triSolveUpperFn {α : Type} [Context α] {n : ℕ} (U : Fin n → Fin n → α) (y : Fin n → α) : Fin n → α

Back substitution: solve $Ux=y$ for an upper-triangular $U$ with nonzero diagonal. Unknowns are visited $n-1,\ldots,1,0$; when row $i$ is reached, acc holds $x_{i+1},\ldots,x_{n-1}$ (and $0$ elsewhere), so $\mathtt{dotFn}(U_i,\mathtt{acc})=\sum_{k>i}U_{ik}x_k$ by upper-triangularity.

def Spec.cholSolveImpl {α : Type} [Context α] {n : ℕ} (L : Fin n → Fin n → α) (b : Fin n → α) : Fin n → α

Strict, array-backed runtime implementation of cholSolveFn (registered via @[implemented_by]). It materializes L into a strict Array (Array α) once, then runs both triangular substitutions over Arrays, so a back-reference is an O(1) lookup. The closure form below (triSolveUpperFn over triSolveLowerFn) is mathematically clean — and is what the correctness proofs reason about — but reads the Function.update accumulator chain on every step, which is ruinous in the interpreter (#eval) when L is itself an unmaterialized closure (e.g. choleskyFn of a kernel matrix). It is intended to compute the same solution strictly; this equivalence is trusted, not proved (see the trust-boundary note above), with the numeric examples (the ridge residual $\approx0$) as evidence rather than a proof.

@[implemented_by Spec.cholSolveImpl]

def Spec.cholSolveFn {α : Type} [Context α] {n : ℕ} (L : Fin n → Fin n → α) (b : Fin n → α) : Fin n → α

Solve $Ax=b$ given a Cholesky factor $L$ of $A$ (so $A=LL^\mathsf{T}$): forward-solve $Lz=b$, then back-solve $L^\mathsf{T}x=z$.

The runtime implementation is cholSolveImpl (strict arrays); the closure form here is what the correctness proofs reason about. The two are intended to compute the same solution — trusted, not proved; see the trust-boundary note above.

def Spec.addScaledIdFn {α : Type} [Context α] {n : ℕ} (K : Fin n → Fin n → α) (γ : α) : Fin n → Fin n → α

The regularized matrix $K+\gamma I$ as a function. For a symmetric PSD kernel $K$ and $\gamma>0$, this is symmetric positive-definite, so its Cholesky factorization succeeds.

def Spec.solveRidgeImpl {α : Type} [Context α] {n : ℕ} (K : Fin n → Fin n → α) (γ : α) (b : Fin n → α) : Fin n → α

Strict, array-backed runtime implementation of solveRidgeFn (registered via @[implemented_by]). It factors $K+\gamma I=LL^\mathsf{T}$ and runs both triangular substitutions entirely over Arrays, so no step materializes the deep Fin n → α closures the functional definition builds — those re-evaluate columns / the substitution accumulator exponentially, which is ruinous in the interpreter (#eval). Intended to be the same linear solve; this equivalence is trusted, not proved (see the trust-boundary note above), with the numeric examples (residual $(K+\gamma I)x-b\approx0$) as evidence rather than a proof.

@[implemented_by Spec.solveRidgeImpl]

def Spec.solveRidgeFn {α : Type} [Context α] {n : ℕ} (K : Fin n → Fin n → α) (γ : α) (b : Fin n → α) : Fin n → α

The Tikhonov-regularized (kernel-ridge) solve $(K+\gamma I)x=b$, via the Cholesky factorization of $K+\gamma I$.

The runtime implementation is solveRidgeImpl (strict arrays); the closure form here, built from the choleskyFn / triSolve* pieces the correctness proofs reason about. The two are intended to compute the same solution — trusted, not proved; see the trust-boundary note above.

def Spec.solveRidgeSpec {α : Type} [Context α] {n : ℕ} (K : Tensor α (Shape.dim n (Shape.dim n Shape.scalar))) (γ : α) (b : Tensor α (Shape.dim n Shape.scalar)) : Tensor α (Shape.dim n Shape.scalar)

Tensor-level kernel-ridge solve: $(K+\gamma I)x=b$.

PyTorch analogue: torch.linalg.solve(K + gamma * I, b) (specialized to the SPD Cholesky path).

QR factorization (classical Gram–Schmidt)

For $A\in\mathbb{R}^{m\times n}$, compute classical Gram–Schmidt factors. Under positive executable $R$ pivots, the proved real theorem gives $A=QR$, with $Q$ having orthonormal columns and $R$ upper-triangular. This uses classical Gram–Schmidt: each $r_{kj}=q_k^\mathsf{T}a_j$ is the inner product against the original column $a_j$, and all projections are subtracted in a single pass (modified Gram–Schmidt would instead dot each $q_k$ against the running residual). In exact real arithmetic the two coincide; the classical form is what the recurrence below implements and what the reconstruction proof matches.

structure Spec.GSState (m n : ℕ) (α : Type) : Type

Internal state for the Gram–Schmidt fold: computed Q columns and R columns so far.

• qs : List (Fin m → α)

  Orthonormal Q columns produced so far (each of length m).

• rcols : List (Fin n → α)

  R columns produced so far (each of length n, upper-triangular).

def Spec.gramSchmidtFn {α : Type} [Context α] {m n : ℕ} (A : Fin m → Fin n → α) : GSState m n α

Run classical Gram–Schmidt over the columns of A, returning the Q columns and R columns. Column j is orthogonalized against the previously produced Q columns.

def Spec.qrQSpec {α : Type} [Context α] {m n : ℕ} (A : Tensor α (Shape.dim m (Shape.dim n Shape.scalar))) : Tensor α (Shape.dim m (Shape.dim n Shape.scalar))

The Q factor candidate of the QR factorization of A. Its columns are proved orthonormal under positive executable R pivots.

def Spec.qrRSpec {α : Type} [Context α] {m n : ℕ} (A : Tensor α (Shape.dim m (Shape.dim n Shape.scalar))) : Tensor α (Shape.dim n (Shape.dim n Shape.scalar))

The R factor (upper-triangular) of the QR factorization of A.

def Spec.qrSpec {α : Type} [Context α] {m n : ℕ} (A : Tensor α (Shape.dim m (Shape.dim n Shape.scalar))) : Tensor α (Shape.dim m (Shape.dim n Shape.scalar)) × Tensor α (Shape.dim n (Shape.dim n Shape.scalar))

QR factorization candidate of $A\in\mathbb{R}^{m\times n}$ via classical Gram–Schmidt. Over ℝ, the full $A=QR$, orthonormal-column, and upper-triangular specification is proved under positive executable $R$ pivots.

PyTorch analogue: torch.linalg.qr(A).