Dot Products And Norms🔗

Both algorithms repeatedly compute inner products. Gram-Schmidt uses them to subtract projections onto earlier columns; Cholesky subtracts inner products of partial rows before computing a diagonal square root. Start with the dot product and norm that implement these steps.

The vector operations live on the plain function representation Fin p → α, not on tensors:

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

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

Context α is TorchLean's scalar backend class: it supplies the arithmetic, the comparisons, and MathFunctions.sqrt. Instantiating α := Float gives a program you can run; instantiating α := ℝ gives an object you can prove things about. That single generic definition is why the proofs and the executable path are not two independent transcriptions of the same recurrence.

The body of dotFn fixes the accumulation order with a left fold:

-- Accumulate products in index order; floating addition
-- cannot freely reassociate this fold.
def dotFn {p : Nat} (u v : Fin p → α) : α :=
  (List.finRange p).foldl (fun s i => s + u i * v i) 0

A mathematician would write \sum_i u_iv_i. The generic executable definition instead fixes a summation order. Float has no AddCommMonoid instance supporting the usual Finset.sum: rounded addition need not be associative, so (a+b)+c and a+(b+c) can differ (Goldberg, 1991). The left fold specifies which result to compute. Over the reals, associativity allows a theorem to identify that fold with the mathematical sum:

@dotFn_eq_sum : ∀ {p : ℕ} (u v : Fin p → ℝ), Spec.dotFn u v = ∑ i, u i * v i

The theorem is specific to ℝ; it does not remove the significance of summation order in a floating-point run. To evaluate the fold, read the indexed coordinates of two length-three tensors:

32.000000

That is 1\cdot4+2\cdot5+3\cdot6=4+10+18=32, and the norm of uVec is \sqrt{14}:

3.741657

At the tensor level the same product has a shaped API, which is the one ordinary programs call. Because the shape is in the type, the two arguments cannot have different lengths:

32.000000

@Tensor.dotSpec : {α : Type} → [inst : Storage α] → [Context α] → {s : Shape} → Tensor α s → Tensor α s → α

The shape parameter is an arbitrary Shape, not [n], because the definition is sumSpec (mulSpec a b): it contracts every scalar leaf. On [n] that is the vector dot product; on [m, n] it is the Frobenius inner product of two matrices.

The last α in the signature is the scalar result type, so no shape remains on the returned value. The preceding two occurrences of Tensor α s share one implicit s; Lean infers it from the first tensor and checks the second against it. Two differently shaped tensors cannot be substituted merely because they contain the same number of scalars. An explicit reshape would be needed to describe which coordinates are to be paired.

Two more bridge lemmas connect the real operations to Euclidean inner products and norms:

@dotFn_eq_inner : ∀ {p : ℕ} (u v : Fin p → ℝ), Spec.dotFn u v = inner ℝ (WithLp.toLp 2 u) (WithLp.toLp 2 v)

@normFn_eq_norm : ∀ {p : ℕ} (v : Fin p → ℝ), Spec.normFn v = ‖WithLp.toLp 2 v‖

In classical Gram-Schmidt, subtracting projections makes each new column orthogonal to its predecessors. To reuse Mathlib's proof of this fact, dotFn_eq_inner and normFn_eq_norm identify the real fold operations with the inner product and norm of EuclideanSpace ℝ (Fin m). Then the orthonormality file proves that the j-th executable column is Mathlib's gramSchmidtNormed ℝ of the column map. Orthonormality then comes straight from Mathlib's gramSchmidtNormed_orthonormal' (Community, 2020). Proving that column identity connects the implementation's lists and folds to the existing analytic theorem.

WithLp.toLp 2 keeps the same coordinate values while viewing them in a space equipped with the Euclidean norm. This matters because a raw function type can carry a different norm convention; a sum-of-squares norm should not be inferred from the presence of real coordinates alone. The bridge establishes the intended geometry explicitly. After that identification, a projection coefficient computed by dotFn is the inner product used by Mathlib's Gram-Schmidt theorem, and a divisor computed by normFn is its Euclidean length. Both identifications are needed: matching only the projection formula would leave open whether the resulting columns were normalized.

Cholesky Recurrence🔗

For a symmetric matrix, Cholesky computes a lower-triangular L using

L_{jj} =\sqrt{A_{jj}-\sum_{k<j}L_{jk}^{\,2}},

and, below the diagonal,

L_{ij} =\frac{A_{ij}-\sum_{k<j}L_{ik}L_{jk}}{L_{jj}} \qquad(i>j).

Entries with i<j are set to zero. In the executable specification, choleskyColsFn is a left fold that appends one column at a time; choleskyFn reads the resulting columns as a matrix; choleskySpec wraps that function as a shaped tensor; and Tensor.cholesky is the public executable entry point.

The fold representation matters in the proof. Once column j is appended, later iterations never change it. Generic “fold that appends” lemmas make that invariant reusable instead of reproving list indexing at every matrix entry.

The two propositions are ordinary definitions, and their types say where they live:

@IsCholesky : {n : ℕ} → Matrix (Fin n) (Fin n) ℝ → Matrix (Fin n) (Fin n) ℝ → Prop

@IsQR : {m k : ℕ} → Matrix (Fin m) (Fin k) ℝ → Matrix (Fin m) (Fin k) ℝ → Matrix (Fin k) (Fin k) ℝ → Prop

Both take real-valued Mathlib Matrix values. Their equalities are exact; a floating-point residual tolerance is a separate property.

The factor itself is the generic function, so it can be read at either scalar type:

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

The first structural theorem needs only an index comparison:

@choleskyFn_lower_triangular : ∀ {n : ℕ} (A : Fin n → Fin n → ℝ) {i j : Fin n}, ↑i < ↑j → Spec.choleskyFn A i j = 0

It applies to every real input matrix, including nonsymmetric matrices and matrices with zero pivots. Applying it takes one line, with norm_num discharging the index comparison 0<2:

The column builder emits 0 whenever i < j, independently of the input values. This proves triangularity but says nothing about whether the remaining entries reconstruct the input.

In the printed triangularity type, {i j : Fin n} are inferred indices, and ↑i < ↑j compares their natural-number values. The result fixes one entry at zero. The example specializes n to three and the indices to zero and two; norm_num proves that index comparison without looking at any matrix entries. This is a useful kind of unconditional theorem because the recurrence explicitly writes zeros above the diagonal. Reconstruction is different: its claim depends on the arithmetic values written below the diagonal, so it cannot follow just from the branch that chooses which entries are zero.

Reconstruction needs two hypotheses:

@isCholesky_of_pos : ∀ {n : ℕ} (A : Fin n → Fin n → ℝ),
  (∀ (i j : Fin n), A i j = A j i) →
    (∀ (j : Fin n), 0 < Spec.choleskyFn A j j) → IsCholesky (Matrix.of A) (Matrix.of (Spec.choleskyFn A))

The positive-pivot condition is the algorithm's success condition. It permits the division by L_{jj} and identifies the positive square root. A symmetric positive-definite matrix is the standard sufficient condition, but the current theorem does not prove

\operatorname{PosDef}(A) \Longrightarrow \forall j,\;0<L_{jj}.

It assumes positivity of the pivots computed by the real specification directly. Deriving those premises from positive definiteness is a remaining proof obligation before the theorem can be applied using only an SPD hypothesis.

The diagonal and off-diagonal equations use positivity in two related ways. A positive pivot is nonzero, permitting cancellation after the division in the off-diagonal recurrence. On the diagonal, the square root must recover the radicand when squared. Over the reals, square root is total and returns zero on a negative argument; assuming its value is positive excludes that case. Symmetry then lets the row identities cover both halves of the input matrix. These are facts about the computed real pivots, which explains why the theorem exposes that particular success condition instead of merely asking for a square input tensor.

At tensor level, choleskySpec_reconstruction states the entrywise identity for L=\operatorname{choleskySpec}(A):

@choleskySpec_reconstruction : ∀ {n : ℕ} (A : Tensor ℝ [n, n]),
  (∀ (i j : Fin n), Spec.get2 A i j = Spec.get2 A j i) →
    (∀ (j : Fin n), 0 < Spec.get2 (Spec.choleskySpec A) j j) →
      ∀ (i j : Fin n), Spec.get2 A i j = ∑ k, Spec.get2 (Spec.choleskySpec A) i k * Spec.get2 (Spec.choleskySpec A) j k

This is A_{ij}=\sum_k L_{ik}L_{jk}, which is A=LL^\top written out one entry at a time. The tensor form exists because callers hold tensors; the hypotheses are the same two as before, restated through Spec.get2.

Cholesky Execution🔗

The same A can be passed to the public executable wrapper Tensor.cholesky:

@Tensor.cholesky : {α : Type} → [inst : Storage α] → [Context α] → {n : ℕ} → Tensor α [n, n] → Tensor α [n, n]

#[2.000000, 0.000000, 0.000000, 1.000000, 2.000000, 0.000000, 1.000000, 1.000000, 2.000000]

Tensor.to ... (Array Float) flattens row-major, so read that as three rows of three: 2 0 0, then 1 2 0, then 1 1 2. The displayed entries agree with the candidate L above.

The recurrence explains this agreement. Every intermediate value in this run is a small integer or a dyadic rational: \sqrt4=2, then 2/2=1, then \sqrt{5-1}=2, then (3-1)/2=1, then \sqrt{6-1-1}=2. Each is exactly representable in binary64 and each operation is exact, so the float run lands on the real answer. Changing A 0 0 from 4 to 3 makes the first pivot \sqrt3; the floating factor then requires a rounding comparison with the real one. The candidate-factor proof is over \mathbb R, while this Float output is a measurement on one input. A formal connection would also need to account for the executable implementation described later in the chapter.

From the repository root, the shipped example checks the residual instead of the entries:

# Run the shipped concrete factorization checks from the
# repository root.
lake exe torchlean factorizations

The Cholesky part of the output is:

Cholesky A = L·Lᵀ: OK (error = 0.000000)
Cholesky on indefinite A correctly fails (no SPD ⇒ no factor): OK (failure detected)

The example itself calls the public executable API:

-- Application code calls the public tensor wrapper on its
-- chosen input matrix.
def lowerFactor := Tensor.cholesky positiveDefiniteMatrix

The second line comes from the symmetric but indefinite matrix

\begin{pmatrix}1&2\\2&1\end{pmatrix},

whose eigenvalues are 3 and -1. Watch where the failure lands:

#[1.000000, 0.000000, 2.000000, NaN]

The first three entries are finite: L_{00}=\sqrt1=1, L_{01}=0 by triangularity, L_{10}=2/1=2. Only the second pivot fails, because it asks for \sqrt{1-2^2}=\sqrt{-3}. The NaN appears at index (1,1), which is exactly the entry isCholesky_of_pos requires to be positive in the real specification. The real square root and the floating square root handle a negative argument differently, so this run illustrates a failed pivot without constituting a direct application of the theorem.

The example intentionally uses a sum of squared entrywise errors for this negative control because IEEE max can ignore a NaN operand. That detail is part of the test's meaning: even a diagnostic norm must choose NaN behavior deliberately.

Cholesky In PyTorch🔗

torch.linalg.cholesky returns the same displayed factor on the positive-definite input (Paszke et al., 2019):

>>> # Compute the factor for the same small
>>> # positive-definite matrix.
>>> import torch
>>> A = torch.tensor([[4.,2.,2.],[2.,5.,3.],[2.,3.,6.]])
>>> torch.linalg.cholesky(A)
tensor([[2., 0., 0.],
        [1., 2., 0.],
        [1., 1., 2.]])

The factor agrees with the Lean result. On the indefinite matrix, PyTorch reports the failed pivot:

>>> # Inspect both the exception and the leading-minor
>>> # diagnostic for the indefinite matrix.
>>> B = torch.tensor([[1.,2.],[2.,1.]])
>>> torch.linalg.cholesky(B)
torch._C._LinAlgError: linalg.cholesky: The factorization could not be
completed because the input is not positive-definite (the leading minor
of order 2 is not positive-definite).
>>> torch.linalg.cholesky_ex(B).info
tensor(2, dtype=torch.int32)

The failure channels differ. PyTorch raises; TorchLean returns a tensor containing NaN. Tensor.cholesky is a total Lean function, so it must return something of type Tensor α [n, n] for every input, and the scalar backend determines the result of invalid arithmetic. The output type carries the dimensions, not a success flag or a proof of finite entries. Returning a tensor therefore says only that the computation produced an inhabitant of that type. It does not discharge the positive-pivot premise. The indefinite example makes this visible without a large matrix: its first column is usable, but the second diagonal entry cannot serve as a real positive square-root pivot. Passing that tensor directly to a later solve would propagate the failed arithmetic beyond the point at which it first became diagnosable. With the Float backend, callers must check for a failed factorization before using its output.

cholesky_ex returns info = 2. That is the one-based index of the first non-positive leading minor, pointing at the same pivot where the Lean run put its NaN. This is runtime diagnostic data. TorchLean's \forall j,\;0<L_{jj} is instead a hypothesis on the real specification's pivots for a chosen matrix. It is not an equivalence theorem between LAPACK's info result and those real pivots; rounding and the different implementation would need their own analysis.

A successful numerical run is commonly described by A\approx LL^\top, but applying such a statement requires an error measure and bound. The Lean theorem here is exact over \mathbb R; the Float examples only measure residuals. This chapter supplies no proved backward-error bound for either displayed floating result. The FP32 soundness chapter describes rounding machinery that a proof for TorchLean's factorization would need.

Classical Gram-Schmidt As QR🔗

For columns a_0,\ldots,a_{k-1}, classical Gram-Schmidt computes

\begin{aligned} v_j &= a_j-\sum_{i<j}\langle q_i,a_j\rangle q_i,\\ r_{jj} &= \|v_j\|,\\ q_j &= v_j/r_{jj},\\ r_{ij} &= \langle q_i,a_j\rangle\quad(i<j). \end{aligned}

Here \langle\cdot,\cdot\rangle and \|\cdot\| correspond to the real instances of dotFn and normFn. Subtracting the projections leaves the part of the new column outside the span of the earlier columns; normalizing that residual gives the next column of Q.

The coefficient uses the original column a_j throughout the projection sum. In exact arithmetic, its component along each earlier orthonormal column can be removed independently, leaving a residual orthogonal to all of them. That is the classical Gram-Schmidt recurrence in this source. If the new column is already in their span, the residual is zero and there is no direction to normalize. A positive r_jj both rules out this degeneracy and selects the division branch of the implementation. Thus the pivot hypothesis describes the step where a new independent feature direction is actually added to the factor.

The TorchLean specification uses the same column-building pattern as Cholesky. gramSchmidtFn threads lists of Q and R columns, while qrQSpec, qrRSpec, and qrSpec expose tensor-shaped results. qrSpec runs Gram-Schmidt once to materialize both factors.

The public operation returns a small structure rather than a pair, so the two factors have names:

TorchLean.Tensor.QRFactors (α : Type) [Storage α] (m n : ℕ) : Type

@Tensor.qr : {α : Type} → [inst : Storage α] → [Context α] → {m n : ℕ} → Tensor α [m, n] → Tensor.QRFactors α m n

The fields carry the reduced shapes:

-- The shared inner dimension allows the two reduced factors
-- to multiply back to the input shape.
result.q : Tensor α [m, min m n]
result.r : Tensor α [min m n, n]

When n ≤ m, those shapes reduce to [m, n] and [n, n], and the public operation uses the theorem-backed qrSpec computation. For a wide matrix, it builds at most m independent basis columns and does not let an early dependent source column consume a basis slot. That distinction is necessary: merely truncating the old square-shaped factors can lose a later independent direction and fail to reconstruct the input. The current IsQR theorem is about the tall/square specification and does not silently certify the separate wide algorithm.

Three separate theorems correspond to the three conjuncts of IsQR. Upper-triangularity, like its Cholesky counterpart, needs no hypothesis:

@Rmat_upper_triangular : ∀ {n m : ℕ} (A : Fin m → Fin n → ℝ) {k j : Fin n}, ↑j < ↑k → Rmat A k j = 0

Rmat and Qmat are the proof layer's readable views of the two factors, defined as plain functions of the row and column index so that induction over columns does not have to travel through a tensor. Orthonormality and reconstruction are where the rank hypothesis enters:

@QT_mul_Q_eq_one : ∀ {m n : ℕ} (A : Fin m → Fin n → ℝ),
  (∀ (j : Fin n), 0 < Rmat A j j) → ((Matrix.of fun i k => Qmat A i k).transpose * Matrix.of fun i k => Qmat A i k) = 1

@isQR_of_pos : ∀ {m n : ℕ} (A : Fin m → Fin n → ℝ),
  (∀ (j : Fin n), 0 < Rmat A j j) →
    IsQR (Matrix.of A) (Matrix.of fun i k => Qmat A i k) (Matrix.of fun k j => Rmat A k j)

\forall j,\;0<R_{jj} is, for classical Gram-Schmidt, the executable form of full column rank: every new column has a nonzero component orthogonal to its predecessors. The API again states the pivot condition directly rather than deriving it from a separately formalized rank predicate, the same unfinished bridge as on the Cholesky side.

At the tensor boundary the two results are restated through Spec.get2:

@qrSpec_reconstruction : ∀ {n m : ℕ} (A : Tensor ℝ [m, n]),
  (∀ (j : Fin n), 0 < Spec.get2 (Spec.qrRSpec A) j j) →
    ∀ (i : Fin m) (j : Fin n), Spec.get2 A i j = ∑ k, Spec.get2 (Spec.qrQSpec A) i k * Spec.get2 (Spec.qrRSpec A) k j

@qrSpec_orthonormal : ∀ {m n : ℕ} (A : Tensor ℝ [m, n]),
  (∀ (j : Fin n), 0 < Spec.get2 (Spec.qrRSpec A) j j) →
    ∀ (a b : Fin n), ∑ i, Spec.get2 (Spec.qrQSpec A) i a * Spec.get2 (Spec.qrQSpec A) i b = if a = b then 1 else 0