3.6 Gaussian quadrature and the Jacobi matrix
Gaussian quadrature combines orthogonality with interpolation. The spectral description of its nodes and weights gives the Golub–Welsch construction.
3.6.1 Exactness and positive weights
Integrating the cardinal polynomials gives an interpolatory rule. Choosing the zeros of an orthogonal polynomial as nodes raises its exactness from degree below \(n\) to degree below \(2n\).
For a measure \(\mu \), finite nodes \(v\) indexed by \(s\), and an index \(i\), the quadrature-weight domain requires \(v_i\ne v_j\) for \(j\in s\setminus \{ i\} \) and \(\mu \)-integrability of \(\ell _i(X)=\prod _{j\in s\setminus \{ i\} }(X-v_j)/(v_i-v_j)\). Only this individual polynomial is required to be integrable.
Lean Formalisation
inductive Lagrange.QuadratureInput.{u_1} {ι : Type u_1} [DecidableEq ι] (μ : MeasureTheory.Measure ℝ) (s : Finset ι) (v : ι → ℝ)
(i : ι) : PropLean definition
/-- Local domain for a classical quadrature integral (DLMF §3.5(iv)). -/
structure QuadratureInput (μ : Measure ℝ) (s : Finset ι) (v : ι → ℝ) (i : ι) where
basis : BasisInput s v i
integrable : Integrable (fun x ↦ (classicalBasis s v i basis).eval x) μFor a measure \(\mu \) on \(\mathbb R\) and a finite node set, define \(w_i=\int \ell _i\, d\mu \) on the domain of Definition 3.6.1, where \(\ell _i\) is the Lagrange basis polynomial of the nodes (Theorem 2.4.15). For a polynomial sequence \(P\), its degree-\(n\) Gaussian weights use the finite set \(s_n\) of distinct real roots of \(P_n\), with the identity node map. No ordering of the roots is chosen. The normalized-product integral also defines a weight at an off-root point whenever it is integrable; the identity node map has no vanishing denominator at any real point.
Reference: [DLMF, 3.5(iv)–(v)].
Lean Formalisation
def Lagrange.quadratureWeight.{u_1} {ι : Type u_1} [DecidableEq ι] (μ : MeasureTheory.Measure ℝ) (s : Finset ι) (v : ι → ℝ)
(i : ι) (d : Lagrange.QuadratureInput μ s v i) : ℝLean definition
/-- The interpolatory quadrature weight on its classical domain (DLMF §3.5(iv)). -/
noncomputable def quadratureWeight (μ : Measure ℝ) (s : Finset ι) (v : ι → ℝ)
(i : ι) (d : QuadratureInput μ s v i) : ℝ :=
∫ x, (classicalBasis s v i d.basis).eval x ∂μAssume finite polynomial moments. For distinct nodes indexed by \(s\) and \(\deg p{\lt}|s|\),
Reference: [DLMF, 3.5(iv)].
Lean Formalisation
theorem Lagrange.integral_eq_sum_quadratureWeight_mul_eval.{u_1} {ι : Type u_1} [DecidableEq ι] {μ : MeasureTheory.Measure ℝ}
(hμ : Polynomial.HasFiniteMoments μ) {s : Finset ι} {v : ι → ℝ} (hv : Set.InjOn v ↑s) {p : Polynomial ℝ}
(hp : p.degree < ↑s.card) :
∫ (x : ℝ), Polynomial.eval x p ∂μ = ∑ i, Lagrange.quadratureWeight μ s v ↑i ⋯ * Polynomial.eval (v ↑i) pLean proof
/-- Integrating the Lagrange interpolant gives an exact rule below the node count
(DLMF §3.5(iv)). -/
theorem integral_eq_sum_quadratureWeight_mul_eval {μ : Measure ℝ}
(hμ : HasFiniteMoments μ) {s : Finset ι} {v : ι → ℝ} (hv : Set.InjOn v s)
{p : ℝ[X]} (hp : p.degree < s.card) :
∫ x, p.eval x ∂μ = ∑ i : s,
quadratureWeight μ s v i (QuadratureInput.ofMoments hμ (BasisInput.of_injOn hv i.property)) *
p.eval (v i) := by
classical
have hi (q : ℝ[X]) : Integrable (fun x ↦ q.eval x) μ := by
simpa using hμ.integrable_mul q 1
calc
_ = ∫ x, ∑ i : s, p.eval (v i) *
(classicalBasis s v i (BasisInput.of_injOn hv i.property)).eval x ∂μ := by
apply integral_congr_ae
filter_upwards [] with x
conv_lhs => rw [eq_sum_classicalBasis hv hp]
simp only [eval_finsetSum, eval_mul, eval_C]
_ = _ := by
rw [integral_finsetSum _ (fun i _ ↦ (hi _).const_mul _)]
simp only [integral_const_mul, quadratureWeight, mul_comm]Let \(P\) be an orthogonal polynomial sequence for \(\mu \), under the finite-moment and positive-definiteness hypotheses. Its root set \(s_n\) has cardinality \(n\). For every polynomial \(p\) with \(\deg p{\lt}2n\),
This includes the zero polynomial and requires no boundedness or connectedness of the support.
Lean Formalisation
theorem Polynomial.Sequence.IsOrthogonal.card_gaussNodes {P : Polynomial.Sequence ℝ} {μ : MeasureTheory.Measure ℝ}
(hμ : Polynomial.HasFiniteMoments μ) (hpos : Polynomial.HasPositiveDefiniteIntegral μ) (hP : P.IsOrthogonal μ)
(n : ℕ) : (↑P n).roots.toFinset.card = nLean proof
/-- There are exactly `n` distinct Gaussian nodes (DLMF §3.5(v)). -/
theorem IsOrthogonal.card_gaussNodes (n : ℕ) : (P n).roots.toFinset.card = n := by
rw [Multiset.toFinset_card_of_nodup (hP.roots_nodup hμ hpos n), hP.card_roots hμ hpos]Let \(P\) be an orthogonal polynomial sequence for a measure \(\mu \) with finite polynomial moments and a positive definite polynomial pairing. Let \(s_n\) be the distinct real roots of \(P_n\) and \(w_x\) their Gaussian weights. Then for every \(x\in s_n\),
Reference: [DLMF, 3.5(v)].
Lean Formalisation
theorem Polynomial.Sequence.IsOrthogonal.gaussWeight_eq_integral_sq {P : Polynomial.Sequence ℝ} {μ : MeasureTheory.Measure ℝ}
(hμ : Polynomial.HasFiniteMoments μ) (hpos : Polynomial.HasPositiveDefiniteIntegral μ) (hP : P.IsOrthogonal μ) {n : ℕ}
{x : ℝ} (hx : x ∈ (↑P n).roots.toFinset) :
P.gaussWeight μ n x ⋯ = ∫ (y : ℝ), Polynomial.eval y (Lagrange.classicalBasis (↑P n).roots.toFinset id x ⋯) ^ 2 ∂μLean proof
/-- A Gaussian weight equals the integral of the square of its Lagrange polynomial
(DLMF §3.5(v)). -/
theorem IsOrthogonal.gaussWeight_eq_integral_sq {n : ℕ} {x : ℝ}
(hx : x ∈ (P n).roots.toFinset) :
P.gaussWeight μ n x (Lagrange.QuadratureInput.ofMoments hμ (Lagrange.BasisInput.id _ _)) =
∫ y, (Lagrange.classicalBasis (P n).roots.toFinset id x (Lagrange.BasisInput.id _ _)).eval y ^ 2 ∂μ := by
classical
let l := Lagrange.basis (P n).roots.toFinset id x
have hinj : Set.InjOn (id : ℝ → ℝ) (P n).roots.toFinset := fun _ _ _ _ h ↦ h
have hn : 0 < n := by
rw [← hP.card_gaussNodes hμ hpos n]
exact Finset.card_pos.mpr ⟨x, hx⟩
have hl : l ≠ 0 := Lagrange.basis_ne_zero hinj hx
have hd : (l ^ 2).degree < (2 * n : ℕ) := by
rw [degree_eq_natDegree (pow_ne_zero _ hl), natDegree_pow]
have hld : l.natDegree = n - 1 := by
simpa [l, hP.card_gaussNodes hμ hpos] using Lagrange.natDegree_basis hinj hx
rw [hld]
exact_mod_cast (show 2 * (n - 1) < 2 * n by omega)
have he := hP.integral_eq_sum_gaussWeight_mul_eval hμ hpos hd
simp only [eval_pow] at he
change _ = ∫ y, l.eval y ^ 2 ∂μ
rw [he]
symm
rw [Finset.sum_eq_single_of_mem x hx]
· have hlx : l.eval x = 1 := Lagrange.eval_basis_self hinj hx
simp only [hlx, one_pow, mul_one]
· intro y hy hyx
have hly : l.eval y = 0 := Lagrange.eval_basis_of_ne (v := id) hyx.symm hy
simp only [hly, ne_eq, OfNat.ofNat_ne_zero, not_false_eq_true, zero_pow, mul_zero]Let \(P\) be an orthogonal polynomial sequence for a measure \(\mu \) with finite polynomial moments and a positive definite polynomial pairing. Let \(s_n\) be the distinct real roots of \(P_n\) and \(w_x\) their Gaussian weights. Then for \(n\ge 1\), \(\sum _{x\in s_n}w_x=\mu (\mathbb R)\), with the finite mass viewed as a real number.
Reference: [DLMF, 3.5(v)].
Lean Formalisation
theorem Polynomial.Sequence.IsOrthogonal.sum_gaussWeight {P : Polynomial.Sequence ℝ} {μ : MeasureTheory.Measure ℝ}
(hμ : Polynomial.HasFiniteMoments μ) (hpos : Polynomial.HasPositiveDefiniteIntegral μ) (hP : P.IsOrthogonal μ) {n : ℕ}
(hn : 1 ≤ n) : ∑ x ∈ (↑P n).roots.toFinset, P.gaussWeight μ n x ⋯ = μ.real Set.univLean proof
/-- Gaussian weights sum to the mass of the measure (DLMF §3.5(v)). -/
theorem IsOrthogonal.sum_gaussWeight {n : ℕ} (hn : 1 ≤ n) :
∑ x ∈ (P n).roots.toFinset, P.gaussWeight μ n x (Lagrange.QuadratureInput.ofMoments hμ (Lagrange.BasisInput.id _ _)) = μ.real Set.univ := by
have hd : (1 : ℝ[X]).degree < (2 * n : ℕ) := by
rw [degree_one]
exact_mod_cast (show 0 < 2 * n by omega)
have he := hP.integral_eq_sum_gaussWeight_mul_eval hμ hpos hd
simpa using he.symmQuadrature-weight domains and formulas
On the quadrature-weight domain, \(\ell _i\) is \(\mu \)-integrable and \(\int \ell _i\, d\mu \) is finite. If \(\mu \) has finite polynomial moments, this integrability holds whenever \(v_i-v_j\ne 0\) for all \(j\in s\setminus \{ i\} \), including admissible \(i\notin s\).
Lean Formalisation
theorem Lagrange.QuadratureInput.integrable_eval.{u_1} {ι : Type u_1} [DecidableEq ι] {μ : MeasureTheory.Measure ℝ}
{s : Finset ι} {v : ι → ℝ} {i : ι} (d : Lagrange.QuadratureInput μ s v i) :
MeasureTheory.Integrable (fun x => Polynomial.eval x (Lagrange.classicalBasis s v i ⋯)) μLean proof
/-- The bundled quadrature polynomial has a finite classical integral. -/
theorem QuadratureInput.integrable_eval {μ : Measure ℝ} {s : Finset ι}
{v : ι → ℝ} {i : ι} (d : QuadratureInput μ s v i) :
Integrable (fun x ↦ (classicalBasis s v i d.basis).eval x) μ := d.integrableFor each index \(i\) in the quadrature-weight domain of Definition 3.6.1, the interpolatory weight is \(\int \ell _i\, d\mu \). For a polynomial sequence, its Gaussian weight is the same integral for the identity nodes given by the distinct real roots.
Lean Formalisation
theorem Lagrange.quadratureWeight_apply.{u_1} {ι : Type u_1} [DecidableEq ι] (μ : MeasureTheory.Measure ℝ) (s : Finset ι)
(v : ι → ℝ) (i : ι) (d : Lagrange.QuadratureInput μ s v i) :
Lagrange.quadratureWeight μ s v i d = ∫ (x : ℝ), Polynomial.eval x (Lagrange.classicalBasis s v i ⋯) ∂μLean proof
/-- Classical quadrature-weight integral formula. -/
theorem quadratureWeight_apply (μ : Measure ℝ) (s : Finset ι) (v : ι → ℝ)
(i : ι) (d : QuadratureInput μ s v i) :
quadratureWeight μ s v i d = ∫ x, (classicalBasis s v i d.basis).eval x ∂μ := rfl3.6.2 Gauss–Jacobi quadrature and error
For Jacobi weights, the nodes lie in \((-1,1)\). Positivity transfers a uniform polynomial approximation bound to a quadrature error bound.
Use the roots of \(P_n^{(\alpha ,\beta )}\) as nodes and integrate their Lagrange basis polynomials against the original Jacobi measure, on the domain where the normalized nodal product is integrable, as in Definition 3.6.1. These weights do not renormalise the measure to have total mass one.
Reference: [DLMF, 3.5(v)].
Lean Formalisation
def Polynomial.jacobiGaussWeight (α β : ℝ) (n : ℕ) (x : ℝ)
(d : Lagrange.QuadratureInput (Polynomial.jacobiMeasure α β) (Polynomial.jacobi α β n).roots.toFinset id x) : ℝLean definition
/-- The Gauss--Jacobi weight, in the original Jacobi measure normalisation
(DLMF §3.5(v)). -/
noncomputable def jacobiGaussWeight (α β : ℝ) (n : ℕ) (x : ℝ)
(d : Lagrange.QuadratureInput (jacobiMeasure α β) (jacobi α β n).roots.toFinset id x) : ℝ :=
Lagrange.quadratureWeight (jacobiMeasure α β) (jacobi α β n).roots.toFinset id x dFor \(\alpha ,\beta {\gt}-1\) and \(n\ge 1\), the \(n\) distinct zeros of \(P_n^{(\alpha ,\beta )}\) lie in \((-1,1)\) and carry strictly positive quadrature weights. The resulting rule is exact for all polynomials of degree at most \(2n-1\) (stated in Lean as \(\deg p{\lt}2n\)), and its weights sum to the total Jacobi mass.
Reference: [DLMF, 3.5(v)].
Lean Formalisation
theorem Polynomial.card_jacobiGaussNodes {α β : ℝ} (hα : -1 < α) (hβ : -1 < β) (n : ℕ) :
(Polynomial.jacobi α β n).roots.toFinset.card = nLean proof
/-- The rule has precisely `n` distinct nodes (DLMF §3.5(v)). -/
theorem card_jacobiGaussNodes (n : ℕ) : (jacobi α β n).roots.toFinset.card = n :=
(isOrthogonal_jacobiSequence hα hβ).card_gaussNodes
(hasFiniteMoments_jacobiMeasure hα hβ)
(hasPositiveDefiniteIntegral_jacobiMeasure hα hβ) nUnder the general Gaussian moment, positivity, and orthogonality hypotheses, let \(n\ge 1\) and let \(f\) be integrable. Any real polynomial \(p\) with \(\deg p{\lt}2n\) and \(|f-p|\le \varepsilon \) both almost everywhere and at all Gaussian nodes gives error at most \(2\mu (\mathbb {R})\varepsilon \). No interval or connected-support hypothesis is needed.
Reference: [ATAP, proof of Theorem 19.3].
Lean Formalisation
theorem Polynomial.Sequence.IsOrthogonal.abs_integral_sub_gauss_le_of_competitor {P : Polynomial.Sequence ℝ}
{μ : MeasureTheory.Measure ℝ} (hP : P.IsOrthogonal μ) (hμ : Polynomial.HasFiniteMoments μ)
(hpos : Polynomial.HasPositiveDefiniteIntegral μ) {n : ℕ} (hn : 1 ≤ n) {f : ℝ → ℝ} (hf : MeasureTheory.Integrable f μ)
{p : Polynomial ℝ} (hp : p.degree < ↑(2 * n)) {ε : ℝ} (hae : ∀ᵐ (x : ℝ) ∂μ, |f x - Polynomial.eval x p| ≤ ε)
(hnodes : ∀ x ∈ (↑P n).roots.toFinset, |f x - Polynomial.eval x p| ≤ ε) :
|∫ (x : ℝ), f x ∂μ - ∑ x ∈ (↑P n).roots.toFinset, P.gaussWeight μ n x ⋯ * f x| ≤ 2 * μ.real Set.univ * εLean proof
/-- A polynomial competitor of degree below `2*n` controls the `n`-node Gaussian
error on any admissible measure. Bounds at nodes are separate from the a.e. bound
(ATAP, proof of Theorem 19.3). -/
theorem IsOrthogonal.abs_integral_sub_gauss_le_of_competitor
{P : Sequence ℝ} {μ : Measure ℝ} (hP : P.IsOrthogonal μ)
(hμ : HasFiniteMoments μ) (hpos : HasPositiveDefiniteIntegral μ)
{n : ℕ} (hn : 1 ≤ n) {f : ℝ → ℝ} (hf : Integrable f μ)
{p : ℝ[X]} (hp : p.degree < (2 * n : ℕ)) {ε : ℝ}
(hae : ∀ᵐ x ∂μ, |f x - p.eval x| ≤ ε)
(hnodes : ∀ x ∈ (P n).roots.toFinset, |f x - p.eval x| ≤ ε) :
|(∫ x, f x ∂μ) - ∑ x ∈ (P n).roots.toFinset, P.gaussWeight μ n x (Lagrange.QuadratureInput.ofMoments hμ (Lagrange.BasisInput.id _ _)) * f x| ≤
2 * μ.real univ * ε := by
have := hμ.isFiniteMeasure
apply abs_integral_sub_sum_le_of_nonneg_weights (P n).roots.toFinset id
(fun x ↦ P.gaussWeight μ n x (Lagrange.QuadratureInput.ofMoments hμ (Lagrange.BasisInput.id _ _))) (fun x hx ↦ (hP.gaussWeight_pos hμ hpos hx).le)
(hP.sum_gaussWeight hμ hpos hn) hf (g := fun x ↦ p.eval x)
· simpa using hμ.integrable_mul p 1
· exact hP.integral_eq_sum_gaussWeight_mul_eval hμ hpos hp
· exact hae
· exact hnodesLet \(\alpha ,\beta {\gt}-1\) and \(n\ge 1\). A uniform polynomial competitor on \([-1,1]\) of degree below \(2n\) and error at most \(\varepsilon \) bounds Gauss–Jacobi quadrature error by \(2\mu _{\alpha ,\beta }(\mathbb {R})\varepsilon \), for an integrable \(f\). In particular, if \(\rho {\gt}1\), \(F\) is analytic on the open ellipse \(E_\rho \), and \(|F|\le M\) there, then for \(f(x)=\Re F(x)\),
The original Jacobi mass is retained. This is a baseline geometric estimate, not the sharper Legendre-specific constant of [ATAP, Theorem 19.3].
Reference: [ATAP, proof of Theorem 19.3].
Lean Formalisation
theorem Polynomial.abs_integral_sub_jacobiGauss_le_of_competitor {α β : ℝ} (hα : -1 < α) (hβ : -1 < β) {n : ℕ} (hn : 1 ≤ n)
{f : ℝ → ℝ} (hf : MeasureTheory.Integrable f (Polynomial.jacobiMeasure α β)) {p : Polynomial ℝ}
(hp : p.degree < ↑(2 * n)) {ε : ℝ} (he : ∀ x ∈ Set.Icc (-1) 1, |f x - Polynomial.eval x p| ≤ ε) :
|∫ (x : ℝ), f x ∂Polynomial.jacobiMeasure α β -
∑ x ∈ (Polynomial.jacobi α β n).roots.toFinset, Polynomial.jacobiGaussWeight α β n x ⋯ * f x| ≤
2 * (Polynomial.jacobiMeasure α β).real Set.univ * εLean proof
/-- Gauss--Jacobi error from a real polynomial competitor (ATAP, proof of 19.3). -/
theorem abs_integral_sub_jacobiGauss_le_of_competitor
{α β : ℝ} (hα : -1 < α) (hβ : -1 < β) {n : ℕ} (hn : 1 ≤ n)
{f : ℝ → ℝ} (hf : Integrable f (jacobiMeasure α β))
{p : ℝ[X]} (hp : p.degree < (2 * n : ℕ)) {ε : ℝ}
(he : ∀ x ∈ Icc (-1 : ℝ) 1, |f x - p.eval x| ≤ ε) :
|(∫ x, f x ∂jacobiMeasure α β) -
∑ x ∈ (jacobi α β n).roots.toFinset, jacobiGaussWeight α β n x (Lagrange.QuadratureInput.ofMoments (hasFiniteMoments_jacobiMeasure hα hβ) (Lagrange.BasisInput.id _ _)) * f x| ≤
2 * (jacobiMeasure α β).real univ * ε := by
apply (isOrthogonal_jacobiSequence hα hβ).abs_integral_sub_gauss_le_of_competitor
(hasFiniteMoments_jacobiMeasure hα hβ)
(hasPositiveDefiniteIntegral_jacobiMeasure hα hβ) hn hf hp
· filter_upwards [ae_mem_Ioo_jacobiMeasure α β] with x hx
exact he x (Ioo_subset_Icc_self hx)
· intro x hx
exact he x (Ioo_subset_Icc_self (jacobiGaussNodes_subset_Ioo hα hβ n hx))Gauss–Jacobi weight formulas
Let \(s\) be the distinct real roots of \(P_n^{(\alpha ,\beta )}\) and let \(x\) lie in the quadrature-weight domain of Definition 3.6.1 for the identity nodes on \(s\) and measure \(\mu _{\alpha ,\beta }\). Then
Lean Formalisation
theorem Polynomial.jacobiGaussWeight_apply (α β : ℝ) (n : ℕ) (x : ℝ)
(d : Lagrange.QuadratureInput (Polynomial.jacobiMeasure α β) (Polynomial.jacobi α β n).roots.toFinset id x) :
Polynomial.jacobiGaussWeight α β n x d =
∫ (y : ℝ),
Polynomial.eval y
(Lagrange.classicalBasis (Polynomial.jacobi α β n).roots.toFinset id x ⋯) ∂Polynomial.jacobiMeasure α βLean proof
/-- Classical Gauss--Jacobi weight integral. -/
theorem jacobiGaussWeight_apply (α β : ℝ) (n : ℕ) (x : ℝ)
(d : Lagrange.QuadratureInput (jacobiMeasure α β) (jacobi α β n).roots.toFinset id x) :
jacobiGaussWeight α β n x d =
∫ y, (Lagrange.classicalBasis (jacobi α β n).roots.toFinset id x d.basis).eval y
∂jacobiMeasure α β := rfl3.6.3 The Jacobi matrix and Golub–Welsch
The recurrence matrix gives the nodes as eigenvalues. After normalization it is symmetric; the first components of its unit eigenvectors give the weights.
For a degree-graded sequence, define \(A_n\) by \((A_n)_{ij}=[xP_j]_i\) for \(0\le i,j{\lt}n\). For an orthogonal sequence we call it the finite Jacobi matrix in the original polynomial normalization.
Reference: [DLMF, 3.5(vi)].
Lean Formalisation
def Polynomial.Sequence.jacobiMatrix.{u_1} {K : Type u_1} [Field K] (P : Polynomial.Sequence K) (n : ℕ) :
Matrix (Fin n) (Fin n) KLean definition
/-- The leading section of multiplication by X, with basis images as columns. -/
noncomputable def jacobiMatrix (n : ℕ) : Matrix (Fin n) (Fin n) K :=
fun i j ↦ (P.basis P.isUnit_leadingCoeff).repr (X * P j) iFor any degree-graded sequence over a field,
Thus the spectrum of \(A_n\) consists precisely of the roots of \(P_n\) in the base field. At a root \(x\), the evaluation vector satisfies \(\sum _{j{\lt}n}(A_n)_{ji}P_j(x)=xP_i(x)\).
Lean Formalisation
theorem Polynomial.Sequence.repr_X_mul_mod_next.{u_1} {K : Type u_1} [Field K] (P : Polynomial.Sequence K) (n : ℕ)
(i j : Fin n) :
((P.basis ⋯).repr (Polynomial.X * ↑P ↑j %ₘ (↑P n * Polynomial.C (↑P n).leadingCoeff⁻¹))) ↑i = P.jacobiMatrix n i jLean proof
/-- Taking the remainder modulo the next sequence polynomial preserves the
lower coordinates of multiplication by X. -/
theorem repr_X_mul_mod_next (n : ℕ) (i j : Fin n) :
(P.basis P.isUnit_leadingCoeff).repr
((X * P j) %ₘ (P n * C (P n).leadingCoeff⁻¹)) i = P.jacobiMatrix n i j := by
let g := P n * C (P n).leadingCoeff⁻¹
have hg : g.Monic := monic_mul_leadingCoeff_inv (P.ne_zero n)
have hd : g.natDegree = n := by
simp [g, natDegree_mul, leadingCoeff_ne_zero.mpr (P.ne_zero n)]
have hq : ((X * P j) /ₘ g).natDegree = 0 := by
rw [natDegree_divByMonic _ hg, hd, natDegree_X_mul (P.ne_zero j), P.natDegree_eq]
omega
have he : g * ((X * P j) /ₘ g) =
((P n).leadingCoeff⁻¹ * (((X * P j) /ₘ g).coeff 0)) • P n := by
rw [eq_C_of_natDegree_eq_zero hq]
simp only [g, smul_eq_C_mul, map_mul, coeff_C_zero]
ring
change (P.basis P.isUnit_leadingCoeff).repr ((X * P j) %ₘ g) i = _
rw [modByMonic_eq_sub_mul_div, he, map_sub, map_smul]
have hb : (P.basis P.isUnit_leadingCoeff).repr (P n) = Finsupp.single n 1 := by
rw [← P.basis_eq_self P.isUnit_leadingCoeff n, Module.Basis.repr_self]
simp [hb, ne_of_gt i.isLt, jacobiMatrix]For a real polynomial sequence, a measure and a finite size \(n\), the matrix domain requires integrability of every mixed product \(xP_jP_i\) and every square \(P_i^2\) for \(i,j{\lt}n\), together with \(h_i=\int P_i^2\, d\mu {\gt}0\). On this domain set
Reference: [DLMF, 3.5(vi)].
Lean Formalisation
inductive Polynomial.Sequence.NormalizedJacobiMatrixInput (P : Polynomial.Sequence ℝ) (μ : MeasureTheory.Measure ℝ) (n : ℕ) : PropLean definition
/-- Local existence data for a normalized Jacobi section. -/
structure NormalizedJacobiMatrixInput (P : Sequence ℝ) (μ : Measure ℝ) (n : ℕ) where
mixed : ∀ i j : Fin n, Integrable (fun t => (X * P j).eval t * (P i).eval t) μ
square : ∀ i : Fin n, IntegrableSquare μ (P i)
positive : ∀ i : Fin n, 0 < P.squaredNorm μ i (square i)Assume finite moments, positive-definite polynomial integration and orthogonality. Then \(J_n\) is symmetric tridiagonal with nonzero adjacent entries and
Its spectrum is exactly the Gaussian nodes. Signs of leading coefficients are unrestricted; adjacent entries need not be positive in this convention.
Reference: [DLMF, 3.5(vi)].
Lean Formalisation
theorem Polynomial.Sequence.IsOrthogonal.normalizedJacobiMatrix_eq_zero {μ : MeasureTheory.Measure ℝ}
(hμ : Polynomial.HasFiniteMoments μ) {P : Polynomial.Sequence ℝ} (hP : P.IsOrthogonal μ) {n : ℕ}
(d : P.NormalizedJacobiMatrixInput μ n) (i j : Fin n) (hij : ↑j + 1 < ↑i ∨ ↑i + 1 < ↑j) :
P.normalizedJacobiMatrix μ n d i j = 0Lean proof
/-- Orthogonality makes the symmetric section tridiagonal. -/
theorem IsOrthogonal.normalizedJacobiMatrix_eq_zero {n : ℕ}
(d : NormalizedJacobiMatrixInput P μ n) (i j : Fin n)
(hij : (j : ℕ) + 1 < i ∨ (i : ℕ) + 1 < j) :
P.normalizedJacobiMatrix μ n d i j = 0 := by
simp only [normalizedJacobiMatrix, hP.integralInner_X_mul_eq_zero hμ hij, zero_div]For a point \(x\) and each coordinate \(i{\lt}n\), the Gauss-vector domain requires integrability of the Gaussian weight’s normalized nodal product \(\ell _x\), nonnegativity of \(w_x=\int \ell _x\, d\mu \), and an integrable polynomial square with strictly positive squared norm \(h_i\). No conditions are required when \(n=0\). A point need not be a Gaussian node if these local conditions hold.
Lean Formalisation
inductive Polynomial.Sequence.GaussEigenvectorInput (P : Polynomial.Sequence ℝ) (μ : MeasureTheory.Measure ℝ) (n : ℕ) (x : ℝ) :
PropLean definition
/-- Classical square-root and integral domains, only at actual vector coordinates.
The empty vector requires no weight or norm data. -/
structure GaussEigenvectorInput (P : Sequence ℝ) (μ : Measure ℝ) (n : ℕ) (x : ℝ) where
weight : Fin n → Lagrange.QuadratureInput μ (P n).roots.toFinset id x
nonneg : ∀ i : Fin n, 0 ≤ P.gaussWeight μ n x (weight i)
square : ∀ i : Fin n, IntegrableSquare μ (P i)
positive : ∀ i : Fin n, 0 < P.squaredNorm μ i (square i)On the domain of Definition 3.6.18, define \(u_i(x)=\sqrt{w_x}P_i(x)/\sqrt{h_i}\), \(i{\lt}n\). The square roots have nonnegative arguments and their denominators are positive. Reference: [DLMF, 3.5(vi)].
Lean Formalisation
def Polynomial.Sequence.gaussEigenvector (P : Polynomial.Sequence ℝ) (μ : MeasureTheory.Measure ℝ) (n : ℕ) (x : ℝ)
(d : P.GaussEigenvectorInput μ n x) : Fin n → ℝLean definition
/-- The Gaussian evaluation vector in normalized polynomial coordinates (DLMF §3.5(vi)). -/
noncomputable def gaussEigenvector (P : Sequence ℝ) (μ : Measure ℝ) (n : ℕ) (x : ℝ)
(d : GaussEigenvectorInput P μ n x) : Fin n → ℝ :=
fun i ↦ Real.sqrt (P.gaussWeight μ n x (d.weight i)) * (P i).eval x /
Real.sqrt (P.squaredNorm μ i (d.square i))Under finite moments, positive-definite polynomial integration and orthogonality, the canonical vector \(u(x)\) at each Gaussian node is a unit eigenvector of \(J_n\) with eigenvalue \(x\). For every real unit eigenvector \(v\) for that node,
Thus the formula is independent of eigenvector sign. These statements concern exact matrices and vectors, without a floating-point eigensolver error claim.
Lean Formalisation
theorem Polynomial.Sequence.IsOrthogonal.normalizedJacobiMatrix_mulVec_gaussEigenvector {μ : MeasureTheory.Measure ℝ}
(hμ : Polynomial.HasFiniteMoments μ) (hpos : Polynomial.HasPositiveDefiniteIntegral μ) {P : Polynomial.Sequence ℝ}
(hP : P.IsOrthogonal μ) {n : ℕ} {x : ℝ} (hx : x ∈ (↑P n).roots.toFinset) :
(P.normalizedJacobiMatrix μ n ⋯).mulVec (P.gaussEigenvector μ n x ⋯) = x • P.gaussEigenvector μ n x ⋯Lean proof
/-- The Gaussian evaluation vector is an eigenvector in normalized coordinates. -/
theorem IsOrthogonal.normalizedJacobiMatrix_mulVec_gaussEigenvector {n : ℕ} {x : ℝ}
(hx : x ∈ (P n).roots.toFinset) :
P.normalizedJacobiMatrix μ n *ᵥ P.gaussEigenvector μ n x (hP.gaussEigenvectorInput hμ hpos hx) =
x • P.gaussEigenvector μ n x (hP.gaussEigenvectorInput hμ hpos hx) := by
have hx0 : (P n).eval x = 0 :=
(mem_roots (P.ne_zero n)).mp (Multiset.mem_toFinset.mp hx)
have hs := P.normalizedJacobiMatrix_isSymm μ n
ext i
change (∑ j : Fin n, P.normalizedJacobiMatrix μ n (d := NormalizedJacobiMatrixInput.ofMoments _ _ _) i j * P.gaussEigenvector μ n x (hP.gaussEigenvectorInput hμ hpos hx) j) = _
have he (j : Fin n) :
P.normalizedJacobiMatrix μ n (d := NormalizedJacobiMatrixInput.ofMoments _ _ _) i j * P.gaussEigenvector μ n x (hP.gaussEigenvectorInput hμ hpos hx) j =
(Real.sqrt (P.gaussWeight μ n x (Lagrange.QuadratureInput.ofMoments hμ (Lagrange.BasisInput.id _ _))) / Real.sqrt (P.squaredNorm μ i (Polynomial.Sequence.IntegrableSquare.of _ _))) *
(P.jacobiMatrix n j i * (P j).eval x) := by
rw [← hs.apply i j, hP.normalizedJacobiMatrix_entry hμ hpos]
dsimp only [gaussEigenvector, IsOrthogonal.gaussEigenvectorInput]
have hi := (Real.sqrt_pos.mpr (squaredNorm_pos hpos P i)).ne'
have hj := (Real.sqrt_pos.mpr (squaredNorm_pos hpos P j)).ne'
field_simp
simp_rw [he]
rw [← Finset.mul_sum, P.sum_jacobiMatrix_mul_eval hx0]
simp only [Pi.smul_apply, smul_eq_mul, gaussEigenvector]
ringLet \(P_n\) be either \(P_n^{(\alpha ,\beta )}\) with \(\alpha ,\beta {\gt}-1\), or \(T_n\), and let \(\mu \) be its original orthogonality measure. The eigenvalues of the normalized Jacobi matrix \(J_n\) are exactly the roots of \(P_n\). For \(n\ge 1\) and every real unit eigenvector \(v\) of \(J_n\) with eigenvalue \(x\), the Gaussian weight is
Reference: [DLMF, 3.5(vi)].
Lean Formalisation
theorem Polynomial.mem_spectrum_jacobi_iff {α β : ℝ} (hα : -1 < α) (hβ : -1 < β) (n : ℕ) (x : ℝ) :
x ∈ spectrum ℝ ((Polynomial.jacobiSequence hα hβ).normalizedJacobiMatrix (Polynomial.jacobiMeasure α β) n ⋯) ↔
(Polynomial.jacobi α β n).IsRoot xLean proof
/-- Gauss--Jacobi nodes are the spectrum of the normalized Jacobi section. -/
theorem mem_spectrum_jacobi_iff (n : ℕ) (x : ℝ) :
x ∈ spectrum ℝ ((jacobiSequence hα hβ).normalizedJacobiMatrix (jacobiMeasure α β) n) ↔
(jacobi α β n).IsRoot x :=
(isOrthogonal_jacobiSequence hα hβ).mem_spectrum_normalizedJacobiMatrix_iff
(hasFiniteMoments_jacobiMeasure hα hβ)
(hasPositiveDefiniteIntegral_jacobiMeasure hα hβ) n xExample (Julia).
For the Legendre measure, the total mass is \(2\). Golub–Welsch therefore gives weights equal to twice the squared first eigenvector components. The example compares the computed rule with exact monomial integrals through degree \(2n-1\) and at the next even degree.
using LinearAlgebra
n = 5
b = [k / sqrt(4k^2 - 1) for k in 1:n-1] # orthonormal recurrence
E = eigen(SymTridiagonal(zeros(n), b))
x, w = E.values, 2 .* E.vectors[1, :].^2
maximum(abs(dot(w, x.^k) - (iseven(k) ? 2/(k+1) : 0)) for k in 0:2n-1)
# residual through the exactness degree
dot(w, x.^(2n)) - 2/(2n+1) # error at the next even degree
Normalization, evaluation vectors and eigenvector uniqueness
Every normalizing denominator \(\sqrt{h_i}\sqrt{h_j}\) is positive, and \((J_n)_{ij}\sqrt{h_i}\sqrt{h_j}=\int xP_jP_i\, d\mu \). The resulting matrix is symmetric. Finite polynomial moments and a positive definite polynomial pairing ensure these integrability and positivity conditions for every finite size. The size-zero domain has no integrability or positivity obligations.
Lean Formalisation
theorem Polynomial.Sequence.NormalizedJacobiMatrixInput.denominator_pos {P : Polynomial.Sequence ℝ}
{μ : MeasureTheory.Measure ℝ} {n : ℕ} (d : P.NormalizedJacobiMatrixInput μ n) (i j : Fin n) :
0 < √(P.squaredNorm μ ↑i ⋯) * √(P.squaredNorm μ ↑j ⋯)Lean proof
/-- All normalization denominators are strictly positive. -/
theorem NormalizedJacobiMatrixInput.denominator_pos {P : Sequence ℝ} {μ : Measure ℝ} {n : ℕ}
(d : NormalizedJacobiMatrixInput P μ n) (i j : Fin n) :
0 < Real.sqrt (P.squaredNorm μ i (d.square i)) * Real.sqrt (P.squaredNorm μ j (d.square j)) :=
mul_pos (Real.sqrt_pos.mpr (d.positive i)) (Real.sqrt_pos.mpr (d.positive j))On the Gauss-vector domain, each \(i{\lt}n\) satisfies \(w_x\ge 0\) and \(h_i=\int P_i^2\, d\mu {\gt}0\). For \(n=0\) the domain imposes no conditions. If \(\mu \) has finite polynomial moments and a positive definite polynomial pairing and \(P\) is orthogonal for \(\mu \), every Gaussian node belongs to this domain.
Lean Formalisation
theorem Polynomial.Sequence.GaussEigenvectorInput.valid {P : Polynomial.Sequence ℝ} {μ : MeasureTheory.Measure ℝ} {n : ℕ}
{x : ℝ} (d : P.GaussEigenvectorInput μ n x) (i : Fin n) : 0 ≤ P.gaussWeight μ n x ⋯ ∧ 0 < P.squaredNorm μ ↑i ⋯Lean proof
/-- The coordinate square roots and their denominators lie in their classical domains. -/
theorem GaussEigenvectorInput.valid {P : Sequence ℝ} {μ : Measure ℝ} {n : ℕ} {x : ℝ}
(d : GaussEigenvectorInput P μ n x) (i : Fin n) :
0 ≤ P.gaussWeight μ n x (d.weight i) ∧ 0 < P.squaredNorm μ i (d.square i) :=
⟨d.nonneg i, d.positive i⟩Every coordinate of the vector equals \(\sqrt{w_x}P_i(x)/\sqrt{h_i}\).
Lean Formalisation
theorem Polynomial.Sequence.gaussEigenvector_apply (P : Polynomial.Sequence ℝ) (μ : MeasureTheory.Measure ℝ) (n : ℕ) (x : ℝ)
(d : P.GaussEigenvectorInput μ n x) (i : Fin n) :
P.gaussEigenvector μ n x d i = √(P.gaussWeight μ n x ⋯) * Polynomial.eval x (↑P ↑i) / √(P.squaredNorm μ ↑i ⋯)Lean proof
/-- Classical coordinate formula for the Gaussian evaluation vector. -/
theorem gaussEigenvector_apply (P : Sequence ℝ) (μ : Measure ℝ) (n : ℕ) (x : ℝ)
(d : GaussEigenvectorInput P μ n x) (i : Fin n) :
P.gaussEigenvector μ n x d i =
Real.sqrt (P.gaussWeight μ n x (d.weight i)) * (P i).eval x /
Real.sqrt (P.squaredNorm μ i (d.square i)) := rflLet \(\mu \) be any finite measure and \(Q(f)=\sum _{i\in s}w_i f(v_i)\) a finite rule with \(w_i\ge 0\) and \(\sum _iw_i=\mu (\Omega )\). If \(f,g\) are integrable, \(Q(g)=\int g\, d\mu \), and \(|f-g|\le \varepsilon \) both almost everywhere and at every sampled node, then
The measured space need not be an interval or even a subset of \(\mathbb {R}\).
Reference: the positive-weight argument in the proof of [ATAP, Theorem 19.3].
Lean Formalisation
theorem MeasureTheory.abs_integral_sub_sum_le_of_nonneg_weights.{u_1, u_2} {α : Type u_1} {ι : Type u_2} [MeasurableSpace α]
{μ : MeasureTheory.Measure α} [MeasureTheory.IsFiniteMeasure μ] (s : Finset ι) (v : ι → α) (w : ι → ℝ) {f g : α → ℝ}
{ε : ℝ} (hw : ∀ i ∈ s, 0 ≤ w i) (hmass : ∑ i ∈ s, w i = μ.real Set.univ) (hf : MeasureTheory.Integrable f μ)
(hg : MeasureTheory.Integrable g μ) (hexact : ∫ (x : α), g x ∂μ = ∑ i ∈ s, w i * g (v i))
(hae : ∀ᵐ (x : α) ∂μ, |f x - g x| ≤ ε) (hnodes : ∀ i ∈ s, |f (v i) - g (v i)| ≤ ε) :
|∫ (x : α), f x ∂μ - ∑ i ∈ s, w i * f (v i)| ≤ 2 * μ.real Set.univ * εLean proof
/-- Positive quadrature transfers a uniform competitor bound to an integration
error bound. The measured space and sampled points need not form an interval
(ATAP, proof of Theorem 19.3). -/
theorem abs_integral_sub_sum_le_of_nonneg_weights
{α ι : Type*} [MeasurableSpace α] {μ : Measure α} [IsFiniteMeasure μ]
(s : Finset ι) (v : ι → α) (w : ι → ℝ) {f g : α → ℝ} {ε : ℝ}
(hw : ∀ i ∈ s, 0 ≤ w i) (hmass : ∑ i ∈ s, w i = μ.real univ)
(hf : Integrable f μ) (hg : Integrable g μ)
(hexact : ∫ x, g x ∂μ = ∑ i ∈ s, w i * g (v i))
(hae : ∀ᵐ x ∂μ, |f x - g x| ≤ ε)
(hnodes : ∀ i ∈ s, |f (v i) - g (v i)| ≤ ε) :
|(∫ x, f x ∂μ) - ∑ i ∈ s, w i * f (v i)| ≤ 2 * μ.real univ * ε := by
have hi : |∫ x, (f x - g x) ∂μ| ≤ ε * μ.real univ := by
simpa only [Real.norm_eq_abs] using (norm_integral_le_of_norm_le_const
(f := fun x ↦ f x - g x) (by simpa only [Real.norm_eq_abs] using hae))
have hs : |∑ i ∈ s, w i * (f (v i) - g (v i))| ≤ ε * μ.real univ := by
calc
_ ≤ ∑ i ∈ s, |w i * (f (v i) - g (v i))| := Finset.abs_sum_le_sum_abs _ _
_ ≤ ∑ i ∈ s, w i * ε := by
apply Finset.sum_le_sum
intro i hi
rw [abs_mul, abs_of_nonneg (hw i hi)]
exact mul_le_mul_of_nonneg_left (hnodes i hi) (hw i hi)
_ = _ := by rw [← Finset.sum_mul, hmass, mul_comm]
have he : (∫ x, f x ∂μ) - ∑ i ∈ s, w i * f (v i) =
(∫ x, (f x - g x) ∂μ) - ∑ i ∈ s, w i * (f (v i) - g (v i)) := by
rw [integral_sub hf hg, hexact]
simp only [mul_sub, Finset.sum_sub_distrib]
ring
rw [he]
exact (abs_sub _ _).trans ((add_le_add hi hs).trans_eq (by ring))Under the same moment, positivity and orthogonality hypotheses, let \(x\) be a Gaussian node of order \(n\ge 1\), with weight \(w_x\). Its cardinal polynomial has coordinates \(w_xP_i(x)/h_i\) for \(i{\lt}n\). Consequently
Moreover \(w_x=\mu (\mathbb R)u_0(x)^2\).
Lean Formalisation
theorem Polynomial.Sequence.IsOrthogonal.gauss_lagrange_repr {μ : MeasureTheory.Measure ℝ} (hμ : Polynomial.HasFiniteMoments μ)
(hpos : Polynomial.HasPositiveDefiniteIntegral μ) {P : Polynomial.Sequence ℝ} (hP : P.IsOrthogonal μ) {n : ℕ} {x : ℝ}
(hx : x ∈ (↑P n).roots.toFinset) (i : Fin n) :
((P.basis ⋯).repr (Lagrange.classicalBasis (↑P n).roots.toFinset id x ⋯)) ↑i =
P.gaussWeight μ n x ⋯ * Polynomial.eval x (↑P ↑i) / P.squaredNorm μ ↑i ⋯Lean proof
/-- Coordinates of a Gaussian cardinal polynomial, from quadrature exactness. -/
theorem IsOrthogonal.gauss_lagrange_repr {n : ℕ} {x : ℝ}
(hx : x ∈ (P n).roots.toFinset) (i : Fin n) :
(P.basis P.isUnit_leadingCoeff).repr (Lagrange.classicalBasis (P n).roots.toFinset id x (Lagrange.BasisInput.id _ _)) i =
P.gaussWeight μ n x (Lagrange.QuadratureInput.ofMoments hμ (Lagrange.BasisInput.id _ _)) * (P i).eval x / P.squaredNorm μ i (Polynomial.Sequence.IntegrableSquare.of _ _) := by
classical
let l := Lagrange.basis (P n).roots.toFinset id x
have hl : l ≠ 0 := Lagrange.basis_ne_zero (fun _ _ _ _ h ↦ h) hx
have hd : l.natDegree = n - 1 := by
simpa [l, hP.card_gaussNodes hμ hpos] using
Lagrange.natDegree_basis (v := id) (fun _ _ _ _ h ↦ h) hx
have hdeg : (l * P i).degree < (2 * n : ℕ) := by
rw [degree_eq_natDegree (mul_ne_zero hl (P.ne_zero i)), natDegree_mul hl (P.ne_zero i),
hd, P.natDegree_eq]
exact_mod_cast (show n - 1 + (i : ℕ) < 2 * n by omega)
rw [hP.basis_repr_eq_integralInner_div hμ hpos]
congr 1
change (∫ y, l.eval y * (P i).eval y ∂μ) = _
simp_rw [← eval_mul]
rw [hP.integral_eq_sum_gaussWeight_mul_eval hμ hpos hdeg,
Finset.sum_eq_single_of_mem x hx]
· have hlx : l.eval x = 1 := Lagrange.eval_basis_self (fun _ _ _ _ h ↦ h) hx
simp only [eval_mul, hlx, one_mul]
· intro y hy hyx
have hly : l.eval y = 0 := Lagrange.eval_basis_of_ne (v := id) hyx.symm hy
simp only [eval_mul, hly, zero_mul, mul_zero]Let \(A\) over a field have size \(n+1\), vanish at \(A_{ij}\) when \(j{\gt}i+1\), and have nonzero superdiagonal entries. An eigenvector with first component zero is zero. Two eigenvectors for the same eigenvalue, one with nonzero first component, are proportional with ratio determined by those components. If both squared-coordinate sums equal one, their first-component squares agree. For real vectors these sums are their Euclidean squared norms.
Reference: [GVL13], eigenvalue methods.
Lean Formalisation
theorem Matrix.eq_zero_of_hessenberg_eigenvector.{u_1} {K : Type u_1} [Field K] {n : ℕ}
{A : Matrix (Fin (n + 1)) (Fin (n + 1)) K} (hband : ∀ (i j : Fin (n + 1)), ↑i + 1 < ↑j → A i j = 0)
(hsub : ∀ (i : Fin n), A i.castSucc i.succ ≠ 0) {v : Fin (n + 1) → K} {x : K} (hv : A.mulVec v = x • v)
(hzero : v 0 = 0) : v = 0Lean proof
/-- An eigenvector whose first component vanishes is zero for an unreduced
lower Hessenberg matrix. -/
theorem eq_zero_of_hessenberg_eigenvector
(hband : ∀ i j : Fin (n + 1), (i : ℕ) + 1 < j → A i j = 0)
(hsub : ∀ i : Fin n, A i.castSucc i.succ ≠ 0)
{v : Fin (n + 1) → K} {x : K} (hv : A *ᵥ v = x • v) (hzero : v 0 = 0) : v = 0 := by
have hall (j : ℕ) (hj : j < n + 1) : v ⟨j, hj⟩ = 0 := by
induction j using Nat.strong_induction_on with
| h j ih =>
by_cases hz : j = 0
· subst j; exact hzero
let k : Fin n := ⟨j - 1, by omega⟩
have hk : v k.castSucc = 0 := ih (j - 1) (by omega) (by omega)
have he := congrFun hv k.castSucc
change (∑ l, A k.castSucc l * v l) = x * v k.castSucc at he
rw [hk, mul_zero, Finset.sum_eq_single k.succ] at he
· have hvj := (mul_eq_zero.mp he).resolve_left (hsub k)
have hjk : k.succ = (⟨j, hj⟩ : Fin (n + 1)) := by ext; dsimp [k]; omega
rwa [hjk] at hvj
· intro l _ hl
by_cases hlj : (l : ℕ) < j
· rw [ih l hlj l.isLt, mul_zero]
· rw [hband k.castSucc l (by
have hne : (l : ℕ) ≠ j := by
intro h
apply hl
ext
dsimp [k]
omega
dsimp [k]
omega), zero_mul]
· simp
ext j
exact hall j j.isLt