5.1 Polynomial measure modification in one variable
Multiplying a weight by a nonnegative polynomial \(u\) changes the orthogonal polynomials. The Gram matrix records the new inner products in the old basis, and a triangular factor describes the corresponding change of coordinates. This is the Gram/Cholesky viewpoint on measure modification.
The connection from the old basis to the new one is upper triangular with bandwidth \(\deg u\), and multiplication by \(u\) has bandwidths \((\deg u,\deg u)\). Leading sections of the connection are the Cholesky factors of the new Gram matrix, which is itself a section of multiplication by \(u\), with the orthonormal and positive-diagonal conventions specified below. This gives a structured way to build new bases from old ones [Gau70, GOS25].
Theorem
5.1.1Univariate polynomial-modification support
A nonzero nonnegative polynomial density \(u\) preserves finite moments and positive definiteness. The connection from an old degree-\(n\) orthogonal polynomial into the modified sequence has zero coefficient at index \(m\) if \(m{\gt}n\) or \(m+\deg u{\lt}n\). For even original and modified measures, opposite parities also vanish. Equal-degree pairings of positive-leading polynomials with the target orthogonal sequence are positive.
For a real measure with finite moments and positive definite polynomial pairing, multiplication by a polynomial \(u\) in an orthogonal basis has lower and upper bandwidth at most \(\operatorname {natDegree}u\) (zero for the zero polynomial). For nonzero \(u\ge 0\) almost everywhere, conversion from an original orthogonal sequence to one for \(u\, d\mu \) has bandwidths \((0,\deg u)\). If \(C\) denotes conversion and \(J_P,J_Q\) multiplication by \(x\) in source and target coordinates, then \(CJ_P=J_QC\) on finitely supported coefficients. This is not a same-size finite-section similarity assertion.
/-- Multiplication by any polynomial of degree d has lower and upper bandwidth d
in an orthogonal basis, independently of the support of the measure. -/
theorem IsOrthogonal.polynomial_mul_isBanded {μ : Measure ℝ}
(hμ : HasFiniteMoments μ) (hpos : HasPositiveDefiniteIntegral μ)
{P : Sequence ℝ} (hP : P.IsOrthogonal μ)
(hu : ∀ i, IsUnit (P i).leadingCoeff) (u : ℝ[X]) :
((P.basis hu).toMatrix (fun j => u * P j)).IsBanded u.natDegree u.natDegree := by
intro i j hij
rw [Basis.toMatrix_apply, hP.basis_repr_eq_integralInner_div hμ hpos hu]
suffices hz : integralInner μ (IntegrablePair.of μ (u * P j) (P i)) = 0 by rw [hz, zero_div]
have hs : integralInner μ (IntegrablePair.of μ (u * P j) (P i)) = integralInner μ (IntegrablePair.of μ (u * P i) (P j)) := by
unfold integralInner
apply integral_congr_ae
filter_upwards [] with x
simp only [eval_mul]
ring
rcases hij with h | h
· apply hP.integralInner_eq_zero_of_degree_lt hμ
have hd : (u * P j).natDegree < i := natDegree_mul_le.trans_lt (by
simpa only [P.natDegree_eq, Nat.add_comm] using h)
exact degree_le_natDegree.trans_lt (by exact_mod_cast hd)
· rw [hs]
apply hP.integralInner_eq_zero_of_degree_lt hμ
have hd : (u * P i).natDegree < j := natDegree_mul_le.trans_lt (by
simpa only [P.natDegree_eq, Nat.add_comm] using h)
exact degree_le_natDegree.trans_lt (by exact_mod_cast hd)
Let \(Q\) be an orthonormal sequence for a measure \(\nu \) with finite moments and positive definite polynomial pairing, and let \(P\) be any degree-indexed sequence. For \(0\le i,j{\lt}N\), set \(C_{ij}=[Q_i]P_j\). Then
The matrix \(C\) is upper triangular. If both sequences have positive leading coefficients, its diagonal is positive and it is the unique such Gram factor. When \(P\) is orthonormal for \(\mu \) and \(\nu =u\, d\mu \) for a nonzero nonnegative polynomial \(u\), this Gram section equals the leading section of full multiplication by \(u\) in the \(P\) basis. If both sequences have positive leading coefficients, every exact unpivoted LU factorisation with unit lower diagonal has the same upper-factor zeros as \(C\).
theorem Polynomial.Sequence.IsOrthogonal.sum_range_repr {μ : MeasureTheory.Measure ℝ} (hμ : Polynomial.HasFiniteMoments μ)
(hpos : Polynomial.HasPositiveDefiniteIntegral μ) {Q : Polynomial.Sequence ℝ} (hQ : Q.IsOrthogonal μ)
(hu : ∀ (i : ℕ), IsUnit (↑Q i).leadingCoeff) (p : Polynomial ℝ) (N : ℕ) (hp : p.natDegree < N) :
∑ k ∈ Finset.range N, ((Q.basis hu).repr p) k • ↑Q k = p
Lean proof
/-- Every polynomial of degree below N expands in the first N basis vectors. -/
theorem IsOrthogonal.sum_range_repr (p : ℝ[X]) (N : ℕ) (hp : p.natDegree < N) :
∑ k ∈ Finset.range N, (Q.basis hu).repr p k • Q k = p := by
classical
have hs : ((Q.basis hu).repr p).support ⊆ Finset.range N := by
intro k hk
apply Finset.mem_range.mpr
by_contra hn
exact Finsupp.mem_support_iff.mp hk
(hQ.basis_repr_eq_zero_of_natDegree_lt hμ hpos hu (lt_of_lt_of_le hp (by omega)))
have he := (Q.basis hu).linearCombination_repr p
simp only [Finsupp.linearCombination_apply, Finsupp.sum, basis_eq_self] at he
conv_rhs => rw [← he]
symm
apply Finset.sum_subset hs
intro k _ hk
rw [Finsupp.notMem_support_iff.mp hk, zero_smul]
/-- Coordinates of a degree-indexed source form an upper triangular leading section. -/
theorem IsOrthogonal.connection_upperTriangular (P : Sequence ℝ) (N : ℕ) :
(Matrix.of (fun i j : Fin N => (Q.basis hu).repr (P j) i)).IsUpperTriangular := by
intro i j hij
exact hQ.basis_repr_eq_zero_of_natDegree_lt hμ hpos hu (by simpa using hij)
/-- Any exact unit-lower/upper LU of the normalized multiplication section has
exactly the connection support; the upper factor differs by row scaling. -/
theorem IsOrthogonal.modification_lu_support
(hPl : ∀ n, 0 < (P n).leadingCoeff) (hQl : ∀ n, 0 < (Q n).leadingCoeff)
(N : ℕ) (L U : Matrix (Fin N) (Fin N) ℝ)
(hL : L.transpose.IsUpperTriangular) (hd : ∀ i, L i i = 1)
(hU : U.IsUpperTriangular)
(hLU : Matrix.of (fun i j : Fin N => (P.basis hPu).repr (u * P j) i) = L * U) :
∀ i j, U i j = 0 ↔ (Q.basis hQu).repr (P j) i = 0 := by
have hm := hμ.polynomialDensityMeasure u hu
have hp := hpos.polynomialDensityMeasure hμ u hu0 hu
have hg := hP.modification_gram hμ hpos u hu0 hu hQ hPu hQu hPn hQn N
have hr := Matrix.upper_eq_of_diagonal_gram_factor _ _ L U (fun _ => 1)
(fun _ => one_ne_zero) (hQ.connection_upperTriangular hm hp hQu P N)
(fun i => (hQ.connection_diag_pos hm hp hQu P hPl hQl i).ne')
(by simpa using hg.symm) hL hd hU hLU
exact hr.2
Successive nonnegative polynomial measure changes multiply their densities. Leading univariate connection sections compose exactly through an intermediate orthogonal sequence, without padding. Their upper bandwidths add; more generally, composition adds lower and upper bandwidths of coefficient maps separately. For finite matrices with compatible labels in an additive commutative group, if source-minus-target differences lie in finite sets \(S\) and \(T\), those of the product lie in \(S+T\). This applies also to pairs of degree labels. The support bound permits additional cancellations and assumes the same intermediate ordered basis for both factors.
Lean Formalisation
theorem Polynomial.withDensity_polynomial_mul (u v : Polynomial ℝ) (μ : MeasureTheory.Measure { x // 0 ≤ Polynomial.eval x u })
(hv : ∀ᵐ (x : ℝ) ∂u.polynomialDensityMeasure μ, 0 ≤ Polynomial.eval x v) :
v.polynomialDensityMeasure (MeasureTheory.Measure.comap Subtype.val (u.polynomialDensityMeasure μ)) =
(u * v).polynomialDensityMeasure (MeasureTheory.Measure.comap Subtype.val (MeasureTheory.Measure.map Subtype.val μ))
Lean proof
/-- Successive nonnegative polynomial modifications multiply their densities.
The second factor is certified only for the first modified measure. -/
theorem withDensity_polynomial_mul (u v : ℝ[X])
(μ : Measure {x : ℝ // 0 ≤ u.eval x})
(hv : ∀ᵐ x ∂polynomialDensityMeasure u μ, 0 ≤ v.eval x) :
polynomialDensityMeasure v ((polynomialDensityMeasure u μ).comap Subtype.val) =
polynomialDensityMeasure (u * v) ((μ.map Subtype.val).comap Subtype.val) := by
-- These zero extensions are proof devices; each polynomial is evaluated as a
-- density only on its certified nonnegative locus.
let w (p : ℝ[X]) (x : ℝ) : ℝ≥0 :=
if h : 0 ≤ p.eval x then NNReal.mk (p.eval x) h else 0
have hw (p : ℝ[X]) : Measurable (w p) := by
apply measurable_coe_nnreal_real_iff.mp
have he : (fun x => (w p x : ℝ)) =
Set.piecewise {x | 0 ≤ p.eval x} (fun x => p.eval x) 0 := by
funext x
by_cases h : 0 ≤ p.eval x <;> simp [w, Set.piecewise, h]
rw [he]
exact p.continuous.measurable.piecewise
(measurableSet_le measurable_const p.continuous.measurable) measurable_const
have hmap (p : ℝ[X]) (η : Measure {x : ℝ // 0 ≤ p.eval x}) :
polynomialDensityMeasure p η =
(η.map Subtype.val).withDensity (fun x => (w p x : ℝ≥0∞)) := by
ext s hs
rw [polynomialDensityMeasure_apply p η hs, withDensity_apply _ hs,
setLIntegral_map hs (hw p).coe_nnreal_ennreal measurable_subtype_coe]
apply setLIntegral_congr_fun (measurable_subtype_coe hs)
intro x _
simp [w, x.property]
have hcomap (p : ℝ[X]) (η : Measure ℝ) (hp : ∀ᵐ x ∂η, 0 ≤ p.eval x) :
polynomialDensityMeasure p (η.comap Subtype.val) =
η.withDensity (fun x => (w p x : ℝ≥0∞)) := by
rw [hmap, map_comap_polynomial_nonneg p η hp]
have hu : ∀ᵐ x ∂μ.map Subtype.val, 0 ≤ u.eval x := by
rw [ae_map_iff measurable_subtype_coe.aemeasurable
(measurableSet_le measurable_const u.continuous.measurable)]
exact Filter.Eventually.of_forall fun x => x.property
have hv' : ∀ᵐ x ∂μ.map Subtype.val, (w u x : ℝ≥0∞) ≠ 0 → 0 ≤ v.eval x := by
rw [hmap] at hv
exact (ae_withDensity_iff (hw u).coe_nnreal_ennreal).mp hv
have huv : ∀ᵐ x ∂μ.map Subtype.val, 0 ≤ (u * v).eval x := by
filter_upwards [hu, hv'] with x hx hvx
by_cases hz : u.eval x = 0
· simp [eval_mul, hz]
· rw [eval_mul]
apply mul_nonneg hx (hvx ?_)
simp only [w, dite_eq_left hx, ENNReal.coe_ne_zero]
intro he
exact hz (congrArg (fun a : ℝ≥0 => (a : ℝ)) he)
rw [hcomap v _ hv, hmap, hcomap (u * v) _ huv,
← withDensity_mul _ (hw u).coe_nnreal_ennreal (hw v).coe_nnreal_ennreal]
apply withDensity_congr_ae
filter_upwards [hu, hv'] with x hx hvx
by_cases hz : u.eval x = 0
· simp [w, hz, eval_mul]
· have hvx' : 0 ≤ v.eval x := hvx (by
simp only [w, dite_eq_left hx, ENNReal.coe_ne_zero]
intro he
exact hz (congrArg (fun a : ℝ≥0 => (a : ℝ)) he))
change (w u x : ℝ≥0∞) * (w v x : ℝ≥0∞) = (w (u * v) x : ℝ≥0∞)
have hp : 0 ≤ (u * v).eval x := by simpa only [eval_mul] using mul_nonneg hx hvx'
dsimp only [w]
rw [dite_eq_left hx, dite_eq_left hvx', dite_eq_left hp, ← ENNReal.coe_mul]
congr 1
apply NNReal.eq
simp only [NNReal.coe_mul, NNReal.coe_mk, eval_mul]
theorem Module.Basis.toMatrix_comp_isBanded.{u_1, u_2, u_3, u_4} {R : Type u_1} {M : Type u_2} {N : Type u_3} {P : Type u_4}
[CommRing R] [AddCommGroup M] [Module R M] [AddCommGroup N] [Module R N] [AddCommGroup P] [Module R P]
(b : Module.Basis ℕ R M) (c : Module.Basis ℕ R N) (d : Module.Basis ℕ R P) (L : M →ₗ[R] N) (K : N →ₗ[R] P)
{l u l' u' : ℕ} (hL : (c.toMatrix fun j => L (b j)).IsBanded l u)
(hK : (d.toMatrix fun j => K (c j)).IsBanded l' u') : (d.toMatrix fun j => K (L (b j))).IsBanded (l + l') (u + u')
Lean proof
/-- Composing banded polynomial operators adds their lower and upper bandwidths. -/
theorem toMatrix_comp_isBanded (b : Basis ℕ R M) (c : Basis ℕ R N)
(d : Basis ℕ R P) (L : M →ₗ[R] N) (K : N →ₗ[R] P)
{l u l' u' : ℕ}
(hL : (c.toMatrix (fun j => L (b j))).IsBanded l u)
(hK : (d.toMatrix (fun j => K (c j))).IsBanded l' u') :
(d.toMatrix (fun j => K (L (b j)))).IsBanded (l+l') (u+u') := by
classical
intro i j hij
have he := c.coefficientMap_apply d K (c.repr (L (b j))) i
have hc : c.coefficientMap d K (c.repr (L (b j))) = d.repr (K (L (b j))) := by
simp [coefficientMap]
rw [hc] at he
rw [toMatrix_apply, he]
apply Finset.sum_eq_zero
intro k hk
dsimp only
by_cases h : j+l < k ∨ k+u < j
· have hz := hL k j h
change c.repr (L (b j)) k = 0 at hz
rw [hz, mul_zero]
· rw [hK i k (by omega), zero_mul]
/-- Two upper-banded changes of polynomial basis have the sum of the bandwidths.
This applies to compatible successive measure modifications without requiring
the intermediate family to be classical. -/
theorem conversion_isBanded_trans (P Q S : Sequence ℝ)
(hP : ∀ i, IsUnit (P i).leadingCoeff)
(hQ : ∀ i, IsUnit (Q i).leadingCoeff)
(hS : ∀ i, IsUnit (S i).leadingCoeff) {d e : ℕ}
(hPQ : ((Q.basis hQ).toMatrix P).IsBanded 0 d)
(hQS : ((S.basis hS).toMatrix Q).IsBanded 0 e) :
((S.basis hS).toMatrix P).IsBanded 0 (d+e) := by
simpa only [basis_eq_self, LinearMap.id_apply, Nat.zero_add] using
(P.basis hP).toMatrix_comp_isBanded (Q.basis hQ) (S.basis hS)
LinearMap.id LinearMap.id
(by simpa only [basis_eq_self, LinearMap.id_apply] using hPQ)
(by simpa only [basis_eq_self, LinearMap.id_apply] using hQS)
theorem Polynomial.Sequence.IsOrthogonal.connection_comp {μ : MeasureTheory.Measure ℝ} (hμ : Polynomial.HasFiniteMoments μ)
(hpos : Polynomial.HasPositiveDefiniteIntegral μ) {Q : Polynomial.Sequence ℝ} (hQ : Q.IsOrthogonal μ)
(hu : ∀ (i : ℕ), IsUnit (↑Q i).leadingCoeff) (P S : Polynomial.Sequence ℝ)
(hS : ∀ (i : ℕ), IsUnit (↑S i).leadingCoeff) (N : ℕ) :
(Matrix.of fun i j => ((S.basis hS).repr (↑P ↑j)) ↑i) =
(Matrix.of fun i j => ((S.basis hS).repr (↑Q ↑j)) ↑i) * Matrix.of fun i j => ((Q.basis hu).repr (↑P ↑j)) ↑i
Lean proof
/-- Leading connection sections compose exactly, with no additional padding. -/
theorem IsOrthogonal.connection_comp (P S : Sequence ℝ)
(hS : ∀ i, IsUnit (S i).leadingCoeff) (N : ℕ) :
Matrix.of (fun i j : Fin N => (S.basis hS).repr (P j) i) =
Matrix.of (fun i j : Fin N => (S.basis hS).repr (Q j) i) *
Matrix.of (fun i j : Fin N => (Q.basis hu).repr (P j) i) := by
ext i j
simp only [Matrix.mul_apply, Matrix.of_apply]
rw [Fin.sum_univ_eq_sum_range
(fun k => (S.basis hS).repr (Q k) i * (Q.basis hu).repr (P j) k) N]
have he := congrArg (fun p => (S.basis hS).repr p i)
(hQ.sum_range_repr hμ hpos hu (P j) N (by simp))
simpa only [map_sum, map_smul, Finsupp.finsetSum_apply, Finsupp.smul_apply,
smul_eq_mul, mul_comm] using he.symm
def Matrix.HasOffsetSupport.{u_1, u_2, u_4, u_5} {ι : Type u_1} {κ : Type u_2} {R : Type u_4} {G : Type u_5}
[AddCommGroup G] [Zero R] (A : Matrix ι κ R) (row : ι → G) (col : κ → G) (S : Finset G) : Prop
Lean definition
/-- A matrix is supported on specified source-minus-target label differences. -/
def HasOffsetSupport [Zero R] (A : Matrix ι κ R) (row : ι → G) (col : κ → G)
(S : Finset G) : Prop :=
∀ i j, col j-row i ∉ S → A i j = 0
theorem Matrix.HasOffsetSupport.mul.{u_1, u_2, u_3, u_4, u_5} {ι : Type u_1} {κ : Type u_2} {τ : Type u_3} {R : Type u_4}
{G : Type u_5} [AddCommGroup G] [DecidableEq G] [Semiring R] [Fintype κ] {A : Matrix ι κ R} {B : Matrix κ τ R}
{row : ι → G} {mid : κ → G} {col : τ → G} {S T : Finset G} (hA : A.HasOffsetSupport row mid S)
(hB : B.HasOffsetSupport mid col T) : (A * B).HasOffsetSupport row col (S + T)
Lean proof
/-- Compatible finite connections compose with Minkowski-sum support. The
conclusion is containment and permits further exact cancellations. -/
theorem HasOffsetSupport.mul [Semiring R] [Fintype κ]
{A : Matrix ι κ R} {B : Matrix κ τ R}
{row : ι → G} {mid : κ → G} {col : τ → G} {S T : Finset G}
(hA : A.HasOffsetSupport row mid S) (hB : B.HasOffsetSupport mid col T) :
(A * B).HasOffsetSupport row col (S+T) := by
intro i j hij
rw [Matrix.mul_apply]
apply Finset.sum_eq_zero
intro k _
by_cases ha : mid k-row i ∈ S
· have hb : col j-mid k ∉ T := by
intro hb
apply hij
exact Finset.mem_add.mpr ⟨_, ha, _, hb, by abel⟩
rw [hB k j hb, mul_zero]
· rw [hA i k ha, zero_mul]
Suppose \(C\) is upper triangular with nonzero diagonal, \(D\) is invertible diagonal, and \(C^TC=DWD^{-1}\). For any exact unpivoted factorization \(W=LU\) with unit lower triangular \(L\) and upper triangular \(U\),
\[ U_{ij}=D_i^{-1}C_{ii}C_{ij}D_j. \]
Such a factorization exists and is unique, and \(U\) and \(C\) have exactly the same zero entries.
Lean Formalisation
theorem Matrix.exists_lu_of_diagonal_gram_factor {n : ℕ} (W C : Matrix (Fin n) (Fin n) ℝ) (D : Fin n → ℝ)
(hD : ∀ (i : Fin n), D i ≠ 0) (hC : C.IsUpperTriangular) (hc : ∀ (i : Fin n), C i i ≠ 0)
(hg : C.transpose * C = Matrix.diagonal D * W * Matrix.diagonal fun i => (D i)⁻¹) :
∃ L U,
L.transpose.IsUpperTriangular ∧
(∀ (i : Fin n), L i i = 1) ∧
U.IsUpperTriangular ∧
W = L * U ∧ (∀ (i j : Fin n), U i j = (D i)⁻¹ * C i i * C i j * D j) ∧ ∀ (i j : Fin n), U i j = 0 ↔ C i j = 0
Lean proof
/-- Diagonal similarity and Cholesky give an unpivoted LU factorization whose
upper factor has exactly the Cholesky support. -/
theorem exists_lu_of_diagonal_gram_factor (W C : Matrix (Fin n) (Fin n) ℝ)
(D : Fin n → ℝ) (hD : ∀ i, D i ≠ 0)
(hC : C.IsUpperTriangular) (hc : ∀ i, C i i ≠ 0)
(hg : C.transpose * C = diagonal D * W * diagonal (fun i => (D i)⁻¹)) :
∃ L U : Matrix (Fin n) (Fin n) ℝ,
L.transpose.IsUpperTriangular ∧ (∀ i, L i i = 1) ∧ U.IsUpperTriangular ∧
W = L * U ∧ (∀ i j, U i j = (D i)⁻¹ * C i i * C i j * D j) ∧
(∀ i j, U i j = 0 ↔ C i j = 0) := by
let L : Matrix (Fin n) (Fin n) ℝ := fun i j => (D i)⁻¹ * C j i * (C j j)⁻¹ * D j
let U : Matrix (Fin n) (Fin n) ℝ := fun i j => (D i)⁻¹ * C i i * C i j * D j
refine ⟨L, U, ?_, ?_, ?_, ?_, fun _ _ => rfl, ?_⟩
· intro i j hij
change (D j)⁻¹ * C i j * (C i i)⁻¹ * D i = 0
rw [hC hij]
ring
· intro i
dsimp [L]
field_simp [hD, hc]
· intro i j hij
change (D i)⁻¹ * C i i * C i j * D j = 0
rw [hC hij]
ring
· ext i j
have hg' := congrArg (fun M => M i j) hg
simp only [diagonal_mul, mul_diagonal] at hg'
simp only [mul_apply, transpose_apply] at hg'
rw [mul_apply]
have hs : (∑ k, L i k * U k j) = (D i)⁻¹ * (∑ k, C k i * C k j) * D j := by
rw [Finset.mul_sum, Finset.sum_mul]
apply Finset.sum_congr rfl
intro k _
dsimp [L, U]
field_simp [hD, hc]
rw [hs, hg']
field_simp [hD, hc]
· intro i j
dsimp [U]
simp only [mul_eq_zero, inv_eq_zero, hD, hc, false_or, or_false]
theorem Matrix.unitLower_upper_factorization_unique {n : ℕ} (L U L' U' : Matrix (Fin n) (Fin n) ℝ)
(hL : L.transpose.IsUpperTriangular) (hL' : L'.transpose.IsUpperTriangular) (hU : U.IsUpperTriangular)
(hU' : U'.IsUpperTriangular) (hd : ∀ (i : Fin n), L i i = 1) (hd' : ∀ (i : Fin n), L' i i = 1)
(hu : ∀ (i : Fin n), U i i ≠ 0) (he : L * U = L' * U') : L = L' ∧ U = U'
Lean proof
/-- Unpivoted LU with unit lower diagonal and nonzero pivots is unique. -/
theorem unitLower_upper_factorization_unique
(L U L' U' : Matrix (Fin n) (Fin n) ℝ)
(hL : L.transpose.IsUpperTriangular) (hL' : L'.transpose.IsUpperTriangular)
(hU : U.IsUpperTriangular) (hU' : U'.IsUpperTriangular)
(hd : ∀ i, L i i = 1) (hd' : ∀ i, L' i i = 1)
(hu : ∀ i, U i i ≠ 0) (he : L * U = L' * U') : L = L' ∧ U = U' := by
have hall (i : Fin n) : (∀ j, U i j = U' i j) ∧ (∀ j, L j i = L' j i) := by
induction i using WellFoundedLT.induction with | ind i ih
have hrow (j : Fin n) : U i j = U' i j := by
have h := congrArg (fun M => M i j) he
simp only [mul_apply] at h
rw [← Finset.add_sum_erase _ _ (Finset.mem_univ i),
← Finset.add_sum_erase _ _ (Finset.mem_univ i), hd, hd', one_mul, one_mul] at h
have hs : (∑ k ∈ Finset.univ.erase i, L i k * U k j) =
∑ k ∈ Finset.univ.erase i, L' i k * U' k j := by
apply Finset.sum_congr rfl
intro k hk
have hki := (Finset.mem_erase.mp hk).1
rcases lt_or_gt_of_ne hki with hlt | hgt
· rw [(ih k hlt).1 j, (ih k hlt).2 i]
· rw [show L i k = 0 from hL hgt, show L' i k = 0 from hL' hgt, zero_mul, zero_mul]
rw [hs] at h
exact add_right_cancel h
refine ⟨hrow, ?_⟩
intro j
have h := congrArg (fun M => M j i) he
simp only [mul_apply] at h
rw [← Finset.add_sum_erase _ _ (Finset.mem_univ i),
← Finset.add_sum_erase _ _ (Finset.mem_univ i)] at h
have hs : (∑ k ∈ Finset.univ.erase i, L j k * U k i) =
∑ k ∈ Finset.univ.erase i, L' j k * U' k i := by
apply Finset.sum_congr rfl
intro k hk
have hki := (Finset.mem_erase.mp hk).1
rcases lt_or_gt_of_ne hki with hlt | hgt
· rw [(ih k hlt).1 i, (ih k hlt).2 j]
· rw [hU hgt, hU' hgt, mul_zero, mul_zero]
rw [hs, ← hrow i] at h
exact mul_right_cancel₀ (hu i) (add_right_cancel h)
exact ⟨by ext i j; exact (hall j).2 i, by ext i j; exact (hall i).1 j⟩
theorem Matrix.upper_eq_of_diagonal_gram_factor {n : ℕ} (W C L U : Matrix (Fin n) (Fin n) ℝ) (D : Fin n → ℝ)
(hD : ∀ (i : Fin n), D i ≠ 0) (hC : C.IsUpperTriangular) (hc : ∀ (i : Fin n), C i i ≠ 0)
(hg : C.transpose * C = Matrix.diagonal D * W * Matrix.diagonal fun i => (D i)⁻¹) (hL : L.transpose.IsUpperTriangular)
(hd : ∀ (i : Fin n), L i i = 1) (hU : U.IsUpperTriangular) (hLU : W = L * U) :
(∀ (i j : Fin n), U i j = (D i)⁻¹ * C i i * C i j * D j) ∧ ∀ (i j : Fin n), U i j = 0 ↔ C i j = 0
Lean proof
/-- Any exact unpivoted LU of the diagonally normalized Gram matrix has the
same upper-factor zeros as Cholesky, with the precise diagonal scaling. -/
theorem upper_eq_of_diagonal_gram_factor
(W C L U : Matrix (Fin n) (Fin n) ℝ) (D : Fin n → ℝ)
(hD : ∀ i, D i ≠ 0) (hC : C.IsUpperTriangular) (hc : ∀ i, C i i ≠ 0)
(hg : C.transpose * C = diagonal D * W * diagonal (fun i => (D i)⁻¹))
(hL : L.transpose.IsUpperTriangular) (hd : ∀ i, L i i = 1)
(hU : U.IsUpperTriangular) (hLU : W = L * U) :
(∀ i j, U i j = (D i)⁻¹ * C i i * C i j * D j) ∧
(∀ i j, U i j = 0 ↔ C i j = 0) := by
obtain ⟨L', U', hL', hd', hU', hLU', hf, hz⟩ :=
exists_lu_of_diagonal_gram_factor W C D hD hC hc hg
have hu' : ∀ i, U' i i ≠ 0 := fun i => by
intro h
exact hc i ((hz i i).mp h)
have he := unitLower_upper_factorization_unique L' U' L U hL' hL hU' hU hd' hd hu'
(hLU'.symm.trans hLU)
rw [he.2] at hf hz
exact ⟨hf, hz⟩
This makes no assertion about floating-point thresholding.
Example (Julia).
Modify the Legendre weight by \(u=1+x^2\). In the orthonormal Legendre basis, the Gram matrix of the new measure is the leading section of \(u(J)\). To get that section exactly, compute \(J^2\) from a padded \(J\) and then truncate. Its Cholesky factor is the connection matrix: upper triangular with bandwidth \(\deg u=2\) and a positive diagonal.
using LinearAlgebra
N = 8
b = [k / sqrt(4k^2 - 1) for k in 1:N+1] # orthonormal Legendre, padded
J = SymTridiagonal(zeros(N+2), b)
W = (I + J^2)[1:N, 1:N] # Gram of P for (1+x^2) dx
R = cholesky(Symmetric(W)).U
maximum(abs(R[i, j]) for i in 1:N, j in 1:N if j > i + 2) # 0.0
Supporting results: triangular Gram factors
Theorem
5.1.6Gram matrices in orthonormal coordinatesMathlib
For vectors \(P_j\) in a finite-dimensional real inner-product space and an orthonormal basis \(Q\), set \(R_{ij}=\langle Q_i,P_j\rangle \). Then \(G(P)=R^\top R\). This identity is independent of polynomial variables or ordering.
Lean Formalisation
theorem Matrix.gram_eq_conjTranspose_mul.{u_1, u_2, u_4, u_5} {E : Type u_1} {n : Type u_2} {𝕜 : Type u_4} [RCLike 𝕜]
[NormedAddCommGroup E] [InnerProductSpace 𝕜 E] {ι : Type u_5} [Fintype ι] (b : OrthonormalBasis ι 𝕜 E) (v : n → E) :
Matrix.gram 𝕜 v =
(Matrix.of fun i j => (b.repr (v j)).ofLp i).conjTranspose * Matrix.of fun i j => (b.repr (v j)).ofLp i
Every ordered basis of a finite-dimensional real inner-product space admits orthonormalization with upper triangular coordinate matrix \(R\) and strictly positive diagonal. It satisfies \(G(P)=R^\top R\). Among upper triangular matrices with positive diagonal, this factor is unique.
theorem InnerProductSpace.inner_gramSchmidtNormed_self.{u_1} {E : Type u_1} [NormedAddCommGroup E] [InnerProductSpace ℝ E]
{n : ℕ} (v : Fin n → E) (i : Fin n) :
inner ℝ (InnerProductSpace.gramSchmidtNormed ℝ v i) (v i) = ‖InnerProductSpace.gramSchmidt ℝ v i‖
Lean proof
/-- The Gram–Schmidt diagonal coefficient is the norm of the residual vector. -/
theorem inner_gramSchmidtNormed_self (v : Fin n → E) (i : Fin n) :
inner ℝ (gramSchmidtNormed ℝ v i) (v i) = ‖gramSchmidt ℝ v i‖ := by
have hi : inner ℝ (gramSchmidt ℝ v i) (v i) = ‖gramSchmidt ℝ v i‖ ^ 2 := by
rw [gramSchmidt_def'' ℝ v i]
rw [inner_add_right, inner_sum]
have hz : (∑ j ∈ Finset.Iio i,
inner ℝ (gramSchmidt ℝ v i)
((inner ℝ (gramSchmidt ℝ v j) (v i) / ‖gramSchmidt ℝ v j‖ ^ 2) •
gramSchmidt ℝ v j)) = 0 := by
apply Finset.sum_eq_zero
intro j hj
rw [inner_smul_right, gramSchmidt_orthogonal ℝ v
(ne_of_gt (Finset.mem_Iio.mp hj)), mul_zero]
simp only [RCLike.ofReal_real_eq_id, id_eq] at *
rw [hz, add_zero, real_inner_self_eq_norm_sq]
rw [gramSchmidtNormed, real_inner_smul_left, hi]
by_cases h : gramSchmidt ℝ v i = 0
· simp [h]
have hn : ‖gramSchmidt ℝ v i‖ ≠ 0 := norm_ne_zero_iff.mpr h
simp only [RCLike.ofReal_real_eq_id, id_eq]
field_simp
theorem InnerProductSpace.gram_eq_transpose_toMatrix_mul.{u_1, u_2} {E : Type u_1} [NormedAddCommGroup E]
[InnerProductSpace ℝ E] {ι : Type u_2} [Fintype ι] (q : OrthonormalBasis ι ℝ E) (v : ι → E) :
Matrix.gram ℝ v = (q.toBasis.toMatrix v).transpose * q.toBasis.toMatrix v
Lean proof
/-- The Gram matrix is the transpose product of the coordinates in an
orthonormal basis, expressed using the existing basis change matrix. -/
theorem gram_eq_transpose_toMatrix_mul {ι : Type*} [Fintype ι]
(q : OrthonormalBasis ι ℝ E) (v : ι → E) :
gram ℝ v = (q.toBasis.toMatrix v).transpose * q.toBasis.toMatrix v := by
have hR : q.toBasis.toMatrix v = Matrix.of (fun i j => q.repr (v j) i) := by
ext i j
exact q.coe_toBasis_repr_apply (v j) i
rw [hR]
simpa only [conjTranspose_eq_transpose_of_trivial] using gram_eq_conjTranspose_mul q v
theorem Matrix.eq_of_pos_upperTriangular_gram_eq {n : ℕ} {R S : Matrix (Fin n) (Fin n) ℝ} (hR : R.IsUpperTriangular)
(hS : S.IsUpperTriangular) (hRp : ∀ (i : Fin n), 0 < R i i) (hSp : ∀ (i : Fin n), 0 < S i i)
(hgram : R.transpose * R = S.transpose * S) : R = S
Lean proof
/-- A real upper triangular Gram factor with positive diagonal is unique.
This is the finite uniqueness statement used for Cholesky measure conversion. -/
theorem eq_of_pos_upperTriangular_gram_eq {n : ℕ}
{R S : Matrix (Fin n) (Fin n) ℝ}
(hR : R.IsUpperTriangular) (hS : S.IsUpperTriangular)
(hRp : ∀ i, 0 < R i i) (hSp : ∀ i, 0 < S i i)
(hgram : R.transpose * R = S.transpose * S) : R = S := by
have hrow : ∀ i, ∀ j, R i j = S i j := by
intro i
induction i using WellFoundedLT.induction with
| ind i ih =>
have hprod (j : Fin n) : R i i * R i j = S i i * S i j := by
have heq := congrArg (fun M : Matrix (Fin n) (Fin n) ℝ => M i j) hgram
simp only [mul_apply, transpose_apply] at heq
have hsum : (∑ k, (R k i * R k j - S k i * S k j)) =
R i i * R i j - S i i * S i j := by
apply Finset.sum_eq_single i
· intro k _ hki
rcases lt_or_gt_of_ne hki with hk | hk
· rw [ih k hk i, ih k hk j, sub_self]
· rw [hR hk, hS hk, zero_mul, zero_mul, sub_self]
· simp
rw [Finset.sum_sub_distrib, heq, sub_self] at hsum
exact sub_eq_zero.mp hsum.symm
have hd : R i i = S i i := by
have hsquare := hprod i
nlinarith [hRp i, hSp i]
intro j
have hj := hprod j
rw [← hd] at hj
exact mul_left_cancel₀ (ne_of_gt (hRp i)) hj
exact Matrix.ext hrow
Repeated polynomial modifications correspond to compositions of the connection maps. In the next chapter, integer parameter shifts of the semiclassical Jacobi family give examples of these modifications, together with sparse differentiation identities.