Matrix factorizations (spec layer)

This file provides shape-indexed scalar-backend recurrences for two matrix factorizations used by Gaussian processes, kernel ridge regression, PCA, and least squares. The factorization theorems apply over ℝ; floating-point execution follows the same operations with the backend's rounding:

• choleskySpec: Cholesky factorization $A=LL^\mathsf{T}$ (lower-triangular $L$), proved for symmetric real 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 and triSolveUpperFn: forward and 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*, over ℝ under positive-pivot success hypotheses and, for Cholesky reconstruction, symmetry. The triangular- and ridge-solve helpers are executable APIs only: their correctness has not been proved. They follow the standard substitution formulas over the function representation and are exercised by #eval examples.

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 TorchLean.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 TorchLean.Tensor and plain functions

def Spec.toMatFn {α : Type} [TorchLean.Storage α] {m n : ℕ} (A : TorchLean.Tensor α [m, n]) : Fin m → Fin n → α

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

Small numeric helpers on the function representation

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

Dot product in increasing index order, with a zero-seeded left fold.

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

Euclidean norm of a length-p vector.

Cholesky factorization

For a symmetric real input whose executable Cholesky pivots are positive, the lower-triangular result satisfies $A=LL^\mathsf{T}$. The reconstruction theorem assumes symmetry and the executable positive-pivot condition. On other scalar backends this definition specifies the recurrence only.

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] attribute. The clean closure form is what the correctness proofs reason about; the internal 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 and runtime code execute the internal array-backed 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.Internal.choleskyColumnsArray {α : Type} [Context α] {n : ℕ} (A : Fin n → Fin n → α) : Array (Array α)

Shared strict Cholesky recurrence used by factorization and ridge solving. 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, which is 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.

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

View the strictly materialized Cholesky columns as finite functions.

@[implemented_by Spec.Internal.choleskyCols]

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 Internal.choleskyCols (strict arrays); the closure form here is the one used by the correctness proofs. The two are intended to compute the same factor. Their equivalence remains part of the trust boundary described 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} [TorchLean.Storage α] [Context α] {n : ℕ} (A : TorchLean.Tensor α [n, n]) : TorchLean.Tensor α [n, n]

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) over ℝ. Execution on every backend still evaluates the full n-term dot, including products with zero. Such terms cannot be skipped for IEEE inputs: ∞ * 0 produces NaN.

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.Internal.solveCholeskyEntries {α : Type} [Context α] (n : ℕ) (entry : ℕ → ℕ → α) (b : Fin n → α) : Fin n → α

Shared forward and backward substitution loops, given a constant-time factor entry reader. Internal.cholSolve supplies a materialized row array; Internal.solveRidge supplies the materialized columns from the shared Cholesky recurrence. 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.

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

Materialize a Cholesky factor once, then use the shared strict substitution loops.

@[implemented_by Spec.Internal.cholSolve]

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 Internal.cholSolve (strict arrays); the closure form here is the logical definition. Neither solve correctness nor equivalence with the runtime implementation has been proved.

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

The regularized matrix $K+\gamma I$ as a function. Over ℝ, a symmetric PSD kernel $K$ and $\gamma>0$ give a symmetric positive-definite matrix. This statement does not cover backend overflow or rounding.

def Spec.Internal.solveRidge {α : 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.Internal.solveRidge]

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 Internal.solveRidge (strict arrays); the logical definition composes choleskyFn with the triangular solves. The factorization theorems do not establish solve correctness or equivalence with this runtime implementation.

def Spec.solveRidgeSpec {α : Type} [TorchLean.Storage α] [Context α] {n : ℕ} (K : TorchLean.Tensor α [n, n]) (γ : α) (b : TorchLean.Tensor α [n]) : TorchLean.Tensor α [n]

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} [TorchLean.Storage α] [Context α] {m n : ℕ} (A : TorchLean.Tensor α [m, n]) : TorchLean.Tensor α [m, n]

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

def Spec.qrRSpec {α : Type} [TorchLean.Storage α] [Context α] {m n : ℕ} (A : TorchLean.Tensor α [m, n]) : TorchLean.Tensor α [n, n]

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

def Spec.qrSpec {α : Type} [TorchLean.Storage α] [Context α] {m n : ℕ} (A : TorchLean.Tensor α [m, n]) : TorchLean.Tensor α [m, n] × TorchLean.Tensor α [n, n]

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).

@[simp]

theorem Spec.qrSpec_fst {α : Type} [TorchLean.Storage α] [Context α] {m n : ℕ} (A : TorchLean.Tensor α [m, n]) : (qrSpec A).1 = qrQSpec A

The first component of the QR pair is the orthonormal factor.

@[simp]

theorem Spec.qrSpec_snd {α : Type} [TorchLean.Storage α] [Context α] {m n : ℕ} (A : TorchLean.Tensor α [m, n]) : (qrSpec A).2 = qrRSpec A

The second component is the upper-triangular factor.

The pair is returned as a plain product rather than a structure, so these two projections are what callers rewrite with instead of destructuring it.