Skip to content

Commit 2455035

Browse files
Add exact matrix factorizations: Cholesky and QR
The exact, finite half of a verified matrix-factorization foundation for classical / scientific-ML models (Gaussian processes, kernel-ridge regression, PCA, least squares). Spec (NN/Spec/Core/Tensor/Factorizations.lean): choleskySpec (A = L·Lᵀ, lower-triangular L, via a column fold with a strict array @[implemented_by]), the triangular and Cholesky linear solves, and qrSpec (A = Q·R by classical Gram–Schmidt), over the readable Fin n → Fin n → α representation with toMatFn/ofMatFn tensor bridges. Proofs over ℝ via Mathlib (NN/Proofs/Tensor/Basic/Factorizations*.lean): - predicates IsCholesky / IsQR; - Cholesky is lower-triangular (from the column fold) and reconstructs A = L·Lᵀ under SPD pivots; - QR reconstructs A = Q·R with Qᵀ·Q = 1 under full column rank, bridged to Mathlib's gramSchmidt. These are exact finite identities — no iteration, no asymptotic caveat — the bedrock the iterative eigensolver/SVD build on separately. Examples (NN/Examples/Factorization/{Cholesky,QR}): #eval witnesses, each a positive reconstruction / orthonormality check plus a negative control (indefinite A fails Cholesky; rank-deficient A breaks QR orthonormality), wrapped in `#guard_msgs (drop info)` so a regression still fails the build without polluting the log. Verification scope: triSolveLower/Upper, cholSolve, and solveRidge are executable APIs only — this change does not prove their correctness. The verified contribution is the factorizations themselves. The strict-array @[implemented_by] substitution is a trusted boundary (examples are evidence, not an equivalence proof). Blueprint: a focused Ch.4 chapter "Matrix Factorizations: Cholesky and QR". sorry/admit/unsafe/opaque-free. Builds green on Lean 4.31.0: lake build NN.Proofs.Tensor.Basic NN.Examples.Factorization
1 parent 21198a1 commit 2455035

11 files changed

Lines changed: 1783 additions & 0 deletions

File tree

NN/Examples/Factorization.lean

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
/-
2+
Copyright (c) 2026 TorchLean
3+
Released under MIT license as described in the file LICENSE.
4+
Authors: TorchLean Team
5+
-/
6+
7+
module
8+
9+
public import NN.Examples.Factorization.Common
10+
public import NN.Examples.Factorization.Cholesky
11+
public import NN.Examples.Factorization.QR
12+
13+
/-!
14+
# Matrix-factorization examples (Cholesky and QR)
15+
16+
Executable `#eval` witnesses for the exact finite factorizations: Cholesky `A = L·Lᵀ` and QR `A = Q·R`
17+
(with `Qᵀ·Q = I`). Each pairs a positive reconstruction/orthonormality check with a negative control,
18+
over `Float`, sorry/admit-free.
19+
-/
Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
/-
2+
Copyright (c) 2026 TorchLean
3+
Released under MIT license as described in the file LICENSE.
4+
Authors: TorchLean Team
5+
-/
6+
7+
module
8+
9+
public import NN.Examples.Factorization.Common
10+
meta import NN.Examples.Factorization.Common
11+
12+
/-!
13+
# Example: Cholesky factorization
14+
15+
`choleskySpec A` returns the lower-triangular `L` with `A = L · Lᵀ` for a symmetric
16+
positive-definite `A`. Here we factor a 3×3 SPD matrix and check the reconstruction error.
17+
-/
18+
19+
@[expose] public section
20+
21+
22+
namespace NN.Examples.Factorization.Cholesky
23+
24+
/-- A symmetric positive-definite test matrix. -/
25+
def A : Spec.Tensor Float (.dim 3 (.dim 3 .scalar)) :=
26+
mkMat [[4, 2, 2],
27+
[2, 5, 3],
28+
[2, 3, 6]]
29+
30+
/-- The Cholesky factor `L` (lower-triangular). -/
31+
def L : Spec.Tensor Float (.dim 3 (.dim 3 .scalar)) := Spec.choleskySpec A
32+
33+
/-- Reconstruction error `‖A - L·Lᵀ‖_max`. -/
34+
def reconErr : Float := maxMatErr A (mm L (tr L))
35+
36+
-- Inspect the diagonal of the factor.
37+
#guard_msgs (drop info) in
38+
#eval vecToList (Spec.ofVecFn (fun i : Fin 3 => Spec.get2 L i i))
39+
40+
-- Compiled assertion: the factorization reconstructs A (fails the build otherwise).
41+
#guard_msgs (drop info) in
42+
#eval assertLt "Cholesky A = L·Lᵀ" reconErr
43+
44+
/-! ## Negative control: the SPD hypothesis is necessary
45+
46+
`isCholesky_of_pos` requires the executable pivots `L[j,j]` to be positive (`0 < choleskyFn A j j`),
47+
which is exactly the success condition over the reals. The matrix below is symmetric but *not*
48+
positive-definite (eigenvalues `3` and `-1`), so the diagonal step takes `√(negative)` and the
49+
reconstruction is `NaN` — never a small error. This documents that the hypothesis genuinely bites. -/
50+
51+
/-- A symmetric but **indefinite** matrix (eigenvalues `{3, -1}`), outside Cholesky's domain. -/
52+
def Abad : Spec.Tensor Float (.dim 2 (.dim 2 .scalar)) :=
53+
mkMat [[1, 2],
54+
[2, 1]]
55+
56+
def Lbad : Spec.Tensor Float (.dim 2 (.dim 2 .scalar)) := Spec.choleskySpec Abad
57+
-- Use the *summed* Frobenius error here, not `maxMatErr`: IEEE `max` ignores `NaN`, whereas the sum
58+
-- propagates the `NaN` produced by `√(negative)`, faithfully reporting that no factor exists.
59+
def reconErrBad : Float := frobSqErr Abad (mm Lbad (tr Lbad))
60+
61+
#guard_msgs (drop info) in
62+
#eval assertReconFails "Cholesky on indefinite A correctly fails (no SPD ⇒ no factor)" reconErrBad
63+
64+
end NN.Examples.Factorization.Cholesky
Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
/-
2+
Copyright (c) 2026 TorchLean
3+
Released under MIT license as described in the file LICENSE.
4+
Authors: TorchLean Team
5+
-/
6+
7+
module
8+
9+
public import NN.Spec.Core.Tensor.Factorizations
10+
11+
/-!
12+
# Factorization examples — shared helpers
13+
14+
Small `Float`-valued helpers used by the matrix-factorization examples (`Cholesky`, `QR`). These
15+
examples are *executable sanity checks*: each one reconstructs the original matrix from its factors and
16+
asserts (via `#eval`) that the maximum entrywise reconstruction error is below a tolerance, so the build
17+
fails if a factorization is wrong.
18+
19+
These run over `Float` (the executable 64-bit runtime scalar), which is the precision the
20+
factorizations target for Gaussian-process / kernel-method use.
21+
-/
22+
23+
@[expose] public section
24+
25+
26+
namespace NN.Examples.Factorization
27+
28+
/-- Build an `m × n` `Float` matrix tensor from a row-major nested list. Missing entries are `0`. -/
29+
def mkMat {m n : Nat} (rows : List (List Float)) : Spec.Tensor Float (.dim m (.dim n .scalar)) :=
30+
Spec.ofMatFn (fun i j => (rows.getD i.val []).getD j.val 0.0)
31+
32+
/-- Maximum entrywise absolute difference between two `m × n` matrices. -/
33+
def maxMatErr {m n : Nat} (A B : Spec.Tensor Float (.dim m (.dim n .scalar))) : Float :=
34+
(List.finRange m).foldl (fun acc i =>
35+
(List.finRange n).foldl
36+
(fun a j => max a (Float.abs (Spec.get2 A i j - Spec.get2 B i j))) acc) 0.0
37+
38+
/-- Matrix product `A · B` (thin wrapper over `matMulSpec`). -/
39+
def mm {m n p : Nat} (A : Spec.Tensor Float (.dim m (.dim n .scalar)))
40+
(B : Spec.Tensor Float (.dim n (.dim p .scalar))) : Spec.Tensor Float (.dim m (.dim p .scalar)) :=
41+
Spec.matMulSpec A B
42+
43+
/-- Matrix transpose. -/
44+
def tr {m n : Nat} (A : Spec.Tensor Float (.dim m (.dim n .scalar))) :
45+
Spec.Tensor Float (.dim n (.dim m .scalar)) :=
46+
Spec.Tensor.matrixTransposeSpec A
47+
48+
/-- Read a vector tensor back out as a `List Float` (for display). -/
49+
def vecToList {n : Nat} (v : Spec.Tensor Float (.dim n .scalar)) : List Float :=
50+
(List.finRange n).map (fun i => Spec.Tensor.toScalar (Spec.get v i))
51+
52+
/-- Squared Frobenius distance `Σ_{i,j} (A_ij - B_ij)²` between two `m × n` matrices. -/
53+
def frobSqErr {m n : Nat} (A B : Spec.Tensor Float (.dim m (.dim n .scalar))) : Float :=
54+
(List.finRange m).foldl (fun acc i =>
55+
(List.finRange n).foldl
56+
(fun a j => let d := Spec.get2 A i j - Spec.get2 B i j; a + d * d) acc) 0.0
57+
58+
/-- Shared tolerance for reconstruction-error assertions. -/
59+
def tol : Float := 1e-6
60+
61+
/--
62+
Compiled **positive** assertion: print `name: OK (err)` when `err < tol`, otherwise raise an
63+
`IO` error so the build/`#eval` fails. Running this through `#eval` evaluates with the compiler
64+
(fast), unlike `#guard`, which forces slow kernel reduction of the whole factorization.
65+
-/
66+
def assertLt (name : String) (err : Float) (tolerance : Float := tol) : IO Unit :=
67+
if err < tolerance then
68+
IO.println s!"{name}: OK (err = {err})"
69+
else
70+
throw (IO.userError s!"{name}: FAIL (err = {err} ≥ tol = {tolerance})")
71+
72+
/--
73+
Compiled **negative-control** assertion: succeeds only when `err ≥ threshold`, i.e. when a property
74+
that *should not* hold is correctly detected as violated. Gives the metric teeth — a reviewer can see
75+
the same `maxMatErr`/residual that reports `0` on a valid factorization reports a large value on an
76+
invalid one, so the positive checks are not vacuous.
77+
-/
78+
def assertGe (name : String) (err : Float) (threshold : Float := 0.5) : IO Unit :=
79+
if err ≥ threshold then
80+
IO.println s!"{name}: OK (correctly rejected, err = {err}{threshold})"
81+
else
82+
throw (IO.userError s!"{name}: FAIL (err = {err} < {threshold}; expected the property to fail)")
83+
84+
/--
85+
Compiled **negative-control** assertion that a reconstruction *fails*: succeeds when the error is not
86+
below `tol` — including the `NaN` produced when a hypothesis is violated (e.g. Cholesky of a
87+
non-positive-definite matrix takes `√(negative)`). Documents that the success hypotheses (SPD pivots,
88+
full column rank) are genuinely necessary.
89+
-/
90+
def assertReconFails (name : String) (err : Float) (tolerance : Float := tol) : IO Unit :=
91+
if err < tolerance then
92+
throw (IO.userError s!"{name}: FAIL (unexpectedly reconstructed, err = {err} < {tolerance})")
93+
else
94+
IO.println s!"{name}: OK (correctly failed, err = {err})"
95+
96+
end NN.Examples.Factorization

NN/Examples/Factorization/QR.lean

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
/-
2+
Copyright (c) 2026 TorchLean
3+
Released under MIT license as described in the file LICENSE.
4+
Authors: TorchLean Team
5+
-/
6+
7+
module
8+
9+
public import NN.Examples.Factorization.Common
10+
meta import NN.Examples.Factorization.Common
11+
12+
/-!
13+
# Example: QR factorization
14+
15+
`qrSpec A` returns `(Q, R)` with `A = Q · R`, `Q` having orthonormal columns and `R`
16+
upper-triangular (classical Gram–Schmidt). We check both `A = Q·R` and `Qᵀ·Q = I`.
17+
-/
18+
19+
@[expose] public section
20+
21+
22+
namespace NN.Examples.Factorization.QR
23+
24+
/-- A 3×3 test matrix (the classic Householder/QR example). -/
25+
def A : Spec.Tensor Float (.dim 3 (.dim 3 .scalar)) :=
26+
mkMat [[12, -51, 4],
27+
[6, 167, -68],
28+
[-4, 24, -41]]
29+
30+
/-- Orthonormal `Q` factor. -/
31+
def Q : Spec.Tensor Float (.dim 3 (.dim 3 .scalar)) := Spec.qrQSpec A
32+
/-- Upper-triangular `R` factor. -/
33+
def R : Spec.Tensor Float (.dim 3 (.dim 3 .scalar)) := Spec.qrRSpec A
34+
35+
/-- Reconstruction error `‖A - Q·R‖_max`. -/
36+
def reconErr : Float := maxMatErr A (mm Q R)
37+
/-- Orthonormality error `‖Qᵀ·Q - I‖_max`. -/
38+
def orthoErr : Float := maxMatErr (mm (tr Q) Q) (Spec.identityTensorSpec 3)
39+
40+
-- Compiled assertions (fail the build otherwise).
41+
#guard_msgs (drop info) in
42+
#eval assertLt "QR A = Q·R" reconErr
43+
#guard_msgs (drop info) in
44+
#eval assertLt "QR Qᵀ·Q = I" orthoErr
45+
46+
/-! ## Negative control: full column rank is necessary for orthonormality
47+
48+
`qrSpec_orthonormal` (`Qᵀ Q = 1`) requires full column rank — positive `R`-pivots
49+
(`0 < R[j,j]`). The matrix below has a dependent column (`col₂ = 2·col₁`), so Gram–Schmidt produces a
50+
**zero** `Q` column where the pivot vanishes: `A = Q·R` still holds, but `Qᵀ Q` has a `0` on the
51+
diagonal, so orthonormality fails. This separates the two guarantees and shows the rank hypothesis
52+
genuinely bites. -/
53+
54+
/-- A rank-2 matrix (`col₂ = 2·col₁`): reconstructs, but `Q` cannot be orthonormal. -/
55+
def Adef : Spec.Tensor Float (.dim 3 (.dim 3 .scalar)) :=
56+
mkMat [[1, 2, 0],
57+
[2, 4, 1],
58+
[1, 2, 0]]
59+
60+
def Qdef : Spec.Tensor Float (.dim 3 (.dim 3 .scalar)) := Spec.qrQSpec Adef
61+
def Rdef : Spec.Tensor Float (.dim 3 (.dim 3 .scalar)) := Spec.qrRSpec Adef
62+
63+
/-- Reconstruction still holds even without full rank. -/
64+
def reconErrDef : Float := maxMatErr Adef (mm Qdef Rdef)
65+
/-- Orthonormality fails: `Qᵀ·Q` has a zero diagonal entry, so it is far from `I`. -/
66+
def orthoErrDef : Float := maxMatErr (mm (tr Qdef) Qdef) (Spec.identityTensorSpec 3)
67+
68+
#guard_msgs (drop info) in
69+
#eval assertLt "QR(rank-deficient) A = Q·R still reconstructs" reconErrDef
70+
#guard_msgs (drop info) in
71+
#eval assertGe "QR(rank-deficient) Qᵀ·Q = I correctly fails (needs full column rank)" orthoErrDef
72+
73+
end NN.Examples.Factorization.QR

NN/Proofs/Tensor/Basic.lean

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,9 @@ module
99
public import NN.Proofs.Tensor.Basic.Core
1010
public import NN.Proofs.Tensor.Basic.Folds
1111
public import NN.Proofs.Tensor.Basic.LinearAlgebra
12+
public import NN.Proofs.Tensor.Basic.Factorizations
13+
public import NN.Proofs.Tensor.Basic.FactorizationsReconstruction
14+
public import NN.Proofs.Tensor.Basic.FactorizationsOrthonormal
1215
public import NN.Proofs.Tensor.Basic.BoundsNorms
1316
public import NN.Proofs.Tensor.Basic.Algebra
1417

0 commit comments

Comments
 (0)