EllipticAlgorithms.lean

1 import KomorebiVerification.PaperIBlockCompletion
2 import KomorebiVerification.SeedRankBasics
3 import KomorebiVerification.SingleOperatorGraph
4 
5 /-!
6 # Paper I: elliptic coordinates and finite algorithm specifications
7 
8 This module checks the displayed two-packet model of
9 `k[x,y] / (y^2 - x^3 + x)`.  The normal form is represented by Mathlib's
10 quadratic algebra over `k[x]`: an element is exactly a pair `a + y b`, and
11 the defining multiplication law makes `y^2 = x^3 - x` hold definitionally.
12 
13 It also gives mathematical specifications and correctness theorems for the
14 three finite algorithms stated in Paper I.  The orbit-closure construction is
15 an ascending chain of subspaces and is proved to terminate in finite
16 dimension.  The transfer-block construction is an actual finite table whose
17 sum is proved to reconstruct the input map.  The graph routine is specified as
18 the maximum exact eigenspace dimension over the finite spectrum; evaluation
19 and performance of a concrete symbolic backend are deliberately outside the
20 scope of these declarations.
21 -/
22 
23 set_option autoImplicit false
24 
25 -- Suppress only style/unused lints in the frozen completion proofs; warnings remain fatal.
26 set_option linter.unnecessarySeqFocus false
27 set_option linter.unusedSectionVars false
28 set_option linter.unusedSimpArgs false
29 
30 namespace Komorebi
31 
32 /-! ## The elliptic coordinate algebra -/
33 
34 namespace EllipticCoordinates
35 
36 open Polynomial
37 
38 noncomputable section
39 
40 universe u
41 
42 variable (k : Type u) [Field k]
43 
44 /-- The coefficient `x^3 - x` in the elliptic relation. -/
45 def curveWeight : Polynomial k := X ^ 3 - X
46 
47 /-- The relation polynomial in the outer variable `y`. -/
48 def relationPolynomial : Polynomial (Polynomial k) :=
49   X ^ 2 - C (curveWeight k)
50 
51 /-- The literal quotient presentation `k[x][y] / (y^2 - (x^3 - x))`. -/
52 abbrev QuotientModel := AdjoinRoot (relationPolynomial k)
53 
54 /--
55 The unique-normal-form model.  An element is exactly `a + y b`, stored as the
56 two fields `re = a` and `im = b`; the quadratic generator squares to
57 `curveWeight k`.
58 -/
59 abbrev NormalForm :=
60   QuadraticAlgebra (Polynomial k) (curveWeight k) 0
61 
62 /-- The distinguished normal-form generator `y = 0 + y * 1`. -/
63 def y : NormalForm k := QuadraticAlgebra.omega
64 
65 /-- Coefficient coordinates `(a + y b) \mapsto (a,b)`, explicitly `k`-linear. -/
66 def coefficientCoordinates :
67     NormalForm k ≃ₗ[k] Polynomial k × Polynomial k where
68   toFun z := (z.re, z.im)
69   invFun p := ⟨p.1, p.2⟩
70   left_inv z := by cases z <;> rfl
71   right_inv p := by cases p <;> rfl
72   map_add' _ _ := rfl
73   map_smul' _ _ := rfl
74 
75 @[simp] theorem coefficientCoordinates_apply (z : NormalForm k) :
76     coefficientCoordinates k z = (z.re, z.im) := rfl
77 
78 @[simp] theorem coefficientCoordinates_symm_apply
79     (p : Polynomial k × Polynomial k) :
80     (coefficientCoordinates k).symm p = ⟨p.1, p.2⟩ := rfl
81 
82 /-- The even packet projector `P₀(a + yb) = a`. -/
83 def P0 : Module.End k (NormalForm k) where
84   toFun z := ⟨z.re, 0⟩
85   map_add' _ _ := by apply QuadraticAlgebra.ext <;> simp
86   map_smul' _ _ := by ext <;> simp
87 
88 /-- The odd packet projector `P₁(a + yb) = yb`. -/
89 def P1 : Module.End k (NormalForm k) where
90   toFun z := ⟨0, z.im⟩
91   map_add' _ _ := by apply QuadraticAlgebra.ext <;> simp
92   map_smul' _ _ := by ext <;> simp
93 
94 @[simp] theorem P0_apply (z : NormalForm k) : P0 k z = ⟨z.re, 0⟩ := by
95   simp [P0]
96 
97 @[simp] theorem P1_apply (z : NormalForm k) : P1 k z = ⟨0, z.im⟩ := by
98   simp [P1]
99 
100 theorem P0_idempotent : (P0 k).comp (P0 k) = P0 k := by
101   apply LinearMap.ext
102   intro z
103   apply QuadraticAlgebra.ext <;> simp [P0]
104 
105 theorem P1_idempotent : (P1 k).comp (P1 k) = P1 k := by
106   apply LinearMap.ext
107   intro z
108   apply QuadraticAlgebra.ext <;> simp [P1]
109 
110 theorem P0_P1_orthogonal : (P0 k).comp (P1 k) = 0 := by
111   apply LinearMap.ext
112   intro z
113   apply QuadraticAlgebra.ext <;> simp [P0, P1]
114 
115 theorem P1_P0_orthogonal : (P1 k).comp (P0 k) = 0 := by
116   apply LinearMap.ext
117   intro z
118   apply QuadraticAlgebra.ext <;> simp [P0, P1]
119 
120 theorem P0_add_P1 : P0 k + P1 k = LinearMap.id := by
121   apply LinearMap.ext
122   intro z
123   apply QuadraticAlgebra.ext <;> simp
124 
125 /-- The two coefficient packets as a complete finite projector family. -/
126 def packetProjectors :
127     Foundation.LinearProjectorFamily k (NormalForm k) Bool where
128   proj i := if i then P1 k else P0 k
129   idempotent i z := by
130     cases i <;> simp [P0, P1]
131   orthogonal {i j} hij z := by
132     cases i <;> cases j
133     · exact (hij rfl).elim
134     · apply QuadraticAlgebra.ext <;> simp [P0, P1]
135     · apply QuadraticAlgebra.ext <;> simp [P0, P1]
136     · exact (hij rfl).elim
137   complete z := by
138     apply QuadraticAlgebra.ext <;> simp [P0, P1]
139 
140 /-- Multiplication by the coefficient variable `x`, in normal-form coordinates. -/
141 def Mx : Module.End k (NormalForm k) where
142   toFun z := ⟨X * z.re, X * z.im⟩
143   map_add' z w := by ext <;> simp [mul_add]
144   map_smul' c z := by ext <;> simp [mul_smul_comm]
145 
146 /-- Multiplication by `y`, in normal-form coordinates. -/
147 def My : Module.End k (NormalForm k) where
148   toFun z := ⟨curveWeight k * z.im, z.re⟩
149   map_add' z w := by ext <;> simp [mul_add]
150   map_smul' c z := by ext <;> simp [mul_smul_comm]
151 
152 @[simp] theorem Mx_apply (z : NormalForm k) :
153     Mx k z = ⟨X * z.re, X * z.im⟩ := rfl
154 
155 @[simp] theorem My_apply (z : NormalForm k) :
156     My k z = ⟨curveWeight k * z.im, z.re⟩ := rfl
157 
158 /-- `M_x` really is multiplication by the class of the coefficient variable. -/
159 theorem Mx_eq_multiplication (z : NormalForm k) :
160     Mx k z =
161       (algebraMap (Polynomial k) (NormalForm k) X) * z := by
162   apply QuadraticAlgebra.ext <;> simp [Mx]
163 
164 /-- `M_y` really is multiplication by the quadratic generator. -/
165 theorem My_eq_multiplication (z : NormalForm k) :
166     My k z = y k * z := by
167   apply QuadraticAlgebra.ext <;> simp [My, y]
168 
169 /-- Multiplication by `y` sends the odd coefficient packet into the even one. -/
170 theorem P0_My (z : NormalForm k) :
171     P0 k (My k z) = ⟨curveWeight k * z.im, 0⟩ := by
172   simp
173 
174 /-- Multiplication by `y` sends the even coefficient packet into the odd one. -/
175 theorem P1_My (z : NormalForm k) :
176     P1 k (My k z) = ⟨0, z.re⟩ := by
177   simp
178 
179 /-- A concrete witness that `P₀` is not linear over the full coordinate algebra. -/
180 theorem P0_not_full_algebra_linear :
181     P0 k (y k * 1) ≠ y k * P0 k 1 := by
182   intro h
183   have him := congrArg QuadraticAlgebra.im h
184   simp only [P0, y, LinearMap.coe_mk, AddHom.coe_mk,
185     QuadraticAlgebra.im_mul, QuadraticAlgebra.im_one,
186     QuadraticAlgebra.re_one, QuadraticAlgebra.im_omega,
187     QuadraticAlgebra.re_omega, mul_zero, zero_mul, one_mul, add_zero, zero_add] at him
188   exact zero_ne_one him
189 
190 /-- A concrete witness that `P₁` is not linear over the full coordinate algebra. -/
191 theorem P1_not_full_algebra_linear :
192     P1 k (y k * 1) ≠ y k * P1 k 1 := by
193   intro h
194   have him := congrArg QuadraticAlgebra.im h
195   simp only [P1, y, LinearMap.coe_mk, AddHom.coe_mk,
196     QuadraticAlgebra.im_mul, QuadraticAlgebra.im_one,
197     QuadraticAlgebra.re_one, QuadraticAlgebra.im_omega,
198     QuadraticAlgebra.re_omega, mul_zero, zero_mul, one_mul, add_zero, zero_add] at him
199   exact one_ne_zero him
200 
201 /--
202 The unit is a genuine one-element module seed when the coordinate algebra acts
203 on itself by multiplication.
204 -/
205 theorem unitSeed_generates :
206     SeedRankBasics.Generates (NormalForm k) (NormalForm k)
207       ({1} : Set (NormalForm k)) := by
208   rw [SeedRankBasics.Generates,
209     Submodule.span_singleton_eq_top_iff]
210   intro z
211   exact ⟨z, by simp⟩
212 
213 /-- Iterating `M_x` is exactly multiplication by the corresponding power of `x`. -/
214 theorem Mx_pow_apply (i : ℕ) (z : NormalForm k) :
215     (Mx k ^ i) z =
216       (algebraMap (Polynomial k) (NormalForm k) X) ^ i * z := by
217   induction i generalizing z with
218   | zero => simp
219   | succ i ih =>
220       rw [pow_succ, Module.End.mul_apply, ih, Mx_eq_multiplication]
221       simp only [pow_succ, mul_assoc]
222 
223 /-- Iterating `M_y` is exactly multiplication by the corresponding power of `y`. -/
224 theorem My_pow_apply (j : ℕ) (z : NormalForm k) :
225     (My k ^ j) z = (y k) ^ j * z := by
226   induction j generalizing z with
227   | zero => simp
228   | succ j ih =>
229       rw [pow_succ, Module.End.mul_apply, ih, My_eq_multiplication]
230       simp only [pow_succ, mul_assoc]
231 
232 /--
233 The displayed address formula: applying `M_x^i M_y^j` to the unit seed gives
234 the quotient-ring monomial `x^i y^j`.
235 -/
236 theorem Mx_pow_My_pow_unitSeed (i j : ℕ) :
237     ((Mx k ^ i) * (My k ^ j)) (1 : NormalForm k) =
238       (algebraMap (Polynomial k) (NormalForm k) X) ^ i * (y k) ^ j := by
239   rw [Module.End.mul_apply, Mx_pow_apply, My_pow_apply]
240   simp only [mul_one]
241 
242 /-- The displayed block-diagonal coefficient formula for `M_x`. -/
243 @[simp] theorem coefficientCoordinates_Mx (z : NormalForm k) :
244     coefficientCoordinates k (Mx k z) =
245       (X * z.re, X * z.im) := rfl
246 
247 /-- The displayed off-diagonal coefficient formula for `M_y`. -/
248 @[simp] theorem coefficientCoordinates_My (z : NormalForm k) :
249     coefficientCoordinates k (My k z) =
250       (curveWeight k * z.im, z.re) := rfl
251 
252 /-- The exact two-packet product law from the paper. -/
253 theorem product_normal_form
254     (a0 a1 b0 b1 : Polynomial k) :
255     (⟨a0, a1⟩ : NormalForm k) * ⟨b0, b1⟩ =
256       ⟨a0 * b0 + curveWeight k * a1 * b1,
257         a0 * b1 + a1 * b0⟩ := by
258   apply QuadraticAlgebra.ext <;> simp
259 
260 /-- The normal-form generator obeys `y^2 = x^3 - x`. -/
261 theorem y_sq :
262     y k ^ 2 =
263       (algebraMap (Polynomial k) (NormalForm k)) (curveWeight k) := by
264   simpa [y, Algebra.smul_def] using
265     (QuadraticAlgebra.omega_pow_two_eq_add
266       (R := Polynomial k) (a := curveWeight k) (b := 0))
267 
268 /-- The normal-form generator is a root of the quotient relation polynomial. -/
269 theorem relationPolynomial_eval_y :
270     (relationPolynomial k).eval₂
271       (algebraMap (Polynomial k) (NormalForm k)) (y k) = 0 := by
272   rw [relationPolynomial, Polynomial.eval₂_sub, Polynomial.eval₂_pow,
273     Polynomial.eval₂_X, Polynomial.eval₂_C, y_sq]
274   exact sub_self _
275 
276 /-- The quotient generator has square `x^3 - x`. -/
277 theorem quotient_root_sq :
278     (AdjoinRoot.root (relationPolynomial k)) ^ 2 =
279       (algebraMap (Polynomial k) (QuotientModel k)) (curveWeight k) := by
280   have h := AdjoinRoot.eval₂_root (relationPolynomial k)
281   rw [relationPolynomial, Polynomial.eval₂_sub, Polynomial.eval₂_pow,
282     Polynomial.eval₂_X, Polynomial.eval₂_C] at h
283   exact sub_eq_zero.mp h
284 
285 /-- The quotient map sending its class of `y` to the normal-form generator. -/
286 def quotientToNormalForm :
287     QuotientModel k →ₐ[Polynomial k] NormalForm k :=
288   AdjoinRoot.liftAlgHom (relationPolynomial k)
289     (Algebra.ofId (Polynomial k) (NormalForm k)) (y k)
290     (relationPolynomial_eval_y k)
291 
292 /-- The normal-form map sending its generator to the quotient class of `y`. -/
293 def normalFormToQuotient :
294     NormalForm k →ₐ[Polynomial k] QuotientModel k :=
295   QuadraticAlgebra.lift
296     ⟨AdjoinRoot.root (relationPolynomial k), by
297       rw [← pow_two, quotient_root_sq]
298       simp [Algebra.smul_def]⟩
299 
300 @[simp] theorem quotientToNormalForm_root :
301     quotientToNormalForm k (AdjoinRoot.root (relationPolynomial k)) = y k := by
302   exact AdjoinRoot.liftAlgHom_root _ _ _ _
303 
304 @[simp] theorem normalFormToQuotient_y :
305     normalFormToQuotient k (y k) =
306       AdjoinRoot.root (relationPolynomial k) := by
307   simp [normalFormToQuotient, y]
308 
309 /--
310 The quotient presentation and the unique coefficient normal form are exactly
311 isomorphic as `k[x]`-algebras.
312 -/
313 def quotientNormalFormEquiv :
314     QuotientModel k ≃ₐ[Polynomial k] NormalForm k where
315   __ := quotientToNormalForm k
316   invFun := normalFormToQuotient k
317   left_inv q := by
318     have hcomp :
319         (normalFormToQuotient k).comp (quotientToNormalForm k) =
320           AlgHom.id (Polynomial k) (QuotientModel k) := by
321       apply AdjoinRoot.algHom_ext
322       simp
323     exact DFunLike.congr_fun hcomp q
324   right_inv z := by
325     have hcomp :
326         (quotientToNormalForm k).comp (normalFormToQuotient k) =
327           AlgHom.id (Polynomial k) (NormalForm k) := by
328       apply QuadraticAlgebra.algHom_ext
329       change quotientToNormalForm k (normalFormToQuotient k (y k)) = y k
330       rw [normalFormToQuotient_y, quotientToNormalForm_root]
331     exact DFunLike.congr_fun hcomp z
332 
333 /--
334 The operator relation from the displayed elliptic example.  The characteristic
335 assumption is the paper's elliptic-curve hypothesis; the quotient identity
336 itself is valid more generally.
337 -/
338 theorem My_sq_eq_Mx_cube_sub_Mx
339     (_hchar : (2 : k) ≠ 0) :
340     My k ^ 2 = Mx k ^ 3 - Mx k := by
341   apply LinearMap.ext
342   intro z
343   have hMy : (My k ^ 2) z =
344       (⟨curveWeight k * z.re, curveWeight k * z.im⟩ : NormalForm k) := by
345     rw [pow_two, Module.End.mul_apply]
346     apply QuadraticAlgebra.ext <;> simp [My]
347   have hMx : (Mx k ^ 3) z =
348       (⟨X ^ 3 * z.re, X ^ 3 * z.im⟩ : NormalForm k) := by
349     rw [show (3 : ℕ) = 2 + 1 by omega, pow_succ, pow_two,
350       Module.End.mul_apply, Module.End.mul_apply]
351     apply QuadraticAlgebra.ext <;> simp [Mx] <;> ring
352   rw [hMy, LinearMap.sub_apply, hMx, Mx_apply]
353   apply QuadraticAlgebra.ext
354   · rw [QuadraticAlgebra.re_sub]
355     simp [curveWeight]
356     ring
357   · rw [QuadraticAlgebra.im_sub]
358     simp [curveWeight]
359     ring
360 
361 /-- In characteristic different from two, the displayed discriminant is nonzero. -/
362 theorem discriminant_ne_zero (hchar : (2 : k) ≠ 0) : (64 : k) ≠ 0 := by
363   rw [show (64 : k) = (2 : k) ^ 6 by norm_num]
364   exact pow_ne_zero 6 hchar
365 
366 end
367 
368 end EllipticCoordinates
369 
370 /-! ## Orbit-hull closure algorithm -/
371 
372 namespace FiniteAlgorithms
373 
374 universe u v w x y
375 
376 section OrbitHull
377 
378 variable {k : Type u} {V : Type v} {ι : Type w}
379 variable [Field k] [AddCommGroup V] [Module k V] [Fintype ι]
380 
381 /-- A subspace is invariant under every operator in the finite input family. -/
382 def InvariantUnder (operators : ι → Module.End k V) (W : Submodule k V) : Prop :=
383   ∀ i, ∀ x ∈ W, operators i x ∈ W
384 
385 /-- The least invariant subspace containing the input seed set. -/
386 def orbitHull (operators : ι → Module.End k V) (seeds : Set V) : Submodule k V :=
387   sInf {W : Submodule k V | seeds ⊆ W ∧ InvariantUnder operators W}
388 
389 theorem seeds_subset_orbitHull
390     (operators : ι → Module.End k V) (seeds : Set V) :
391     seeds ⊆ orbitHull operators seeds := by
392   intro x hx
393   apply (Submodule.mem_sInf).2
394   intro W hW
395   exact hW.1 hx
396 
397 theorem orbitHull_invariant
398     (operators : ι → Module.End k V) (seeds : Set V) :
399     InvariantUnder operators (orbitHull operators seeds) := by
400   intro i x hx
401   apply (Submodule.mem_sInf).2
402   intro W hW
403   exact hW.2 i x ((sInf_le hW) hx)
404 
405 theorem orbitHull_le
406     (operators : ι → Module.End k V) (seeds : Set V)
407     {W : Submodule k V} (hSeeds : seeds ⊆ W)
408     (hInvariant : InvariantUnder operators W) :
409     orbitHull operators seeds ≤ W := by
410   exact sInf_le ⟨hSeeds, hInvariant⟩
411 
412 /-- One closure pass: retain the old subspace and add all operator images. -/
413 def orbitClosureStep (operators : ι → Module.End k V)
414     (W : Submodule k V) : Submodule k V :=
415   W ⊔ ⨆ i, W.map (operators i)
416 
417 theorem le_orbitClosureStep
418     (operators : ι → Module.End k V) (W : Submodule k V) :
419     W ≤ orbitClosureStep operators W :=
420   le_sup_left
421 
422 theorem orbitClosureStep_mono
423     (operators : ι → Module.End k V) :
424     Monotone (orbitClosureStep operators) := by
425   intro W Z hWZ
426   apply sup_le_sup hWZ
427   exact iSup_mono fun i => Submodule.map_mono hWZ
428 
429 /-- The successive finite closure passes, starting from the span of the seeds. -/
430 def orbitClosureIter (operators : ι → Module.End k V)
431     (seeds : Set V) : ℕ → Submodule k V
432   | 0 => Submodule.span k seeds
433   | n + 1 => orbitClosureStep operators (orbitClosureIter operators seeds n)
434 
435 theorem orbitClosureIter_le_succ
436     (operators : ι → Module.End k V) (seeds : Set V) (n : ℕ) :
437     orbitClosureIter operators seeds n ≤
438       orbitClosureIter operators seeds (n + 1) := by
439   exact le_orbitClosureStep operators _
440 
441 theorem orbitClosureIter_monotone
442     (operators : ι → Module.End k V) (seeds : Set V) :
443     Monotone (orbitClosureIter operators seeds) := by
444   exact monotone_nat_of_le_succ (orbitClosureIter_le_succ operators seeds)
445 
446 /-- The closure iterates bundled as an ascending chain. -/
447 def orbitClosureChain (operators : ι → Module.End k V)
448     (seeds : Set V) : ℕ →o Submodule k V where
449   toFun := orbitClosureIter operators seeds
450   monotone' := orbitClosureIter_monotone operators seeds
451 
452 theorem seeds_subset_orbitClosureIter
453     (operators : ι → Module.End k V) (seeds : Set V) (n : ℕ) :
454     seeds ⊆ orbitClosureIter operators seeds n := by
455   intro x hx
456   exact (orbitClosureIter_monotone operators seeds (Nat.zero_le n))
457     (Submodule.subset_span hx)
458 
459 theorem orbitClosureIter_le_orbitHull
460     (operators : ι → Module.End k V) (seeds : Set V) (n : ℕ) :
461     orbitClosureIter operators seeds n ≤ orbitHull operators seeds := by
462   induction n with
463   | zero =>
464       exact Submodule.span_le.mpr (seeds_subset_orbitHull operators seeds)
465   | succ n ih =>
466       rw [orbitClosureIter]
467       apply sup_le ih
468       refine iSup_le fun i => ?_
469       rw [Submodule.map_le_iff_le_comap]
470       intro x hx
471       exact orbitHull_invariant operators seeds i x (ih hx)
472 
473 theorem invariant_of_orbitClosureStep_eq
474     (operators : ι → Module.End k V) {W : Submodule k V}
475     (hStop : orbitClosureStep operators W = W) :
476     InvariantUnder operators W := by
477   intro i x hx
478   rw [← hStop]
479   exact ((le_iSup (fun j => W.map (operators j)) i).trans le_sup_right)
480     ⟨x, hx, rfl⟩
481 
482 /-- A stopped closure pass is exactly the operator hull. -/
483 theorem orbitClosureIter_eq_orbitHull_of_stop
484     (operators : ι → Module.End k V) (seeds : Set V) (n : ℕ)
485     (hStop : orbitClosureIter operators seeds (n + 1) =
486       orbitClosureIter operators seeds n) :
487     orbitClosureIter operators seeds n = orbitHull operators seeds := by
488   apply le_antisymm
489   · exact orbitClosureIter_le_orbitHull operators seeds n
490   · apply orbitHull_le operators seeds
491     · exact seeds_subset_orbitClosureIter operators seeds n
492     · apply invariant_of_orbitClosureStep_eq operators
493       simpa [orbitClosureIter] using hStop
494 
495 /--
496 In finite dimension the ascending closure algorithm stops after finitely many
497 passes, and its stopped output is exactly the least invariant seed hull.
498 -/
499 theorem orbitHull_algorithm_terminates_and_correct
500     [FiniteDimensional k V]
501     (operators : ι → Module.End k V) (seeds : Set V) :
502     ∃ n : ℕ,
503       orbitClosureIter operators seeds (n + 1) =
504         orbitClosureIter operators seeds n ∧
505       orbitClosureIter operators seeds n = orbitHull operators seeds := by
506   have hEventually := eventuallyConst_of_isNoetherian
507     (orbitClosureChain operators seeds)
508   rw [Filter.eventuallyConst_atTop] at hEventually
509   obtain ⟨n, hn⟩ := hEventually
510   have hStop : orbitClosureIter operators seeds (n + 1) =
511       orbitClosureIter operators seeds n := by
512     exact hn (n + 1) (Nat.le_succ n)
513   exact ⟨n, hStop,
514     orbitClosureIter_eq_orbitHull_of_stop operators seeds n hStop⟩
515 
516 end OrbitHull
517 
518 /-! ## Transfer-block algorithm -/
519 
520 section TransferBlocks
521 
522 variable {R : Type u} {M : Type v} {N : Type w}
523 variable {ι : Type x} {κ : Type y}
524 variable [Semiring R]
525 variable [AddCommMonoid M] [Module R M]
526 variable [AddCommMonoid N] [Module R N]
527 variable [Fintype ι] [Fintype κ]
528 
529 /-- The finite table produced by the transfer-block algorithm. -/
530 def transferBlockTable
531     (P : Foundation.LinearProjectorFamily R M ι)
532     (Q : Foundation.LinearProjectorFamily R N κ)
533     (T : M →ₗ[R] N) : κ → ι → (M →ₗ[R] N) :=
534   fun j i => Foundation.transferBlock P Q T j i
535 
536 /-- Reassemble a finite transfer-block table by summing every entry. -/
537 def assembleTransferBlockTable
538     (B : κ → ι → (M →ₗ[R] N)) : M →ₗ[R] N :=
539   ∑ j, ∑ i, B j i
540 
541 /-- Correctness of the transfer-block algorithm: its table reconstructs `T`. -/
542 theorem transferBlock_algorithm_correct
543     (P : Foundation.LinearProjectorFamily R M ι)
544     (Q : Foundation.LinearProjectorFamily R N κ)
545     (T : M →ₗ[R] N) :
546     assembleTransferBlockTable (transferBlockTable P Q T) = T := by
547   exact Foundation.transferBlock_reconstruction_linearMap P Q T
548 
549 /-- The computed zero pattern is exactly the transfer graph by construction. -/
550 def transferGraph
551     (P : Foundation.LinearProjectorFamily R M ι)
552     (Q : Foundation.LinearProjectorFamily R N κ)
553     (T : M →ₗ[R] N) : κ → ι → Prop :=
554   fun j i => transferBlockTable P Q T j i ≠ 0
555 
556 theorem transferGraph_iff_block_ne_zero
557     (P : Foundation.LinearProjectorFamily R M ι)
558     (Q : Foundation.LinearProjectorFamily R N κ)
559     (T : M →ₗ[R] N) (j : κ) (i : ι) :
560     transferGraph P Q T j i ↔
561       Foundation.transferBlock P Q T j i ≠ 0 :=
562   Iff.rfl
563 
564 end TransferBlocks
565 
566 /-! ## Exact graph spectral-rank algorithm -/
567 
568 section GraphRank
569 
570 open scoped Matrix
571 open Polynomial
572 
573 variable {K : Type u} {ν : Type v}
574 variable [Field K] [Fintype ν] [DecidableEq ν]
575 
576 /-- The exact finite spectrum of a square matrix. -/
577 noncomputable def graphSpectrum (A : Matrix ν ν K) : Finset K :=
578   (Matrix.finite_spectrum A).toFinset
579 
580 @[simp] theorem mem_graphSpectrum_iff (A : Matrix ν ν K) (μ : K) :
581     μ ∈ graphSpectrum A ↔ μ ∈ spectrum K A := by
582   classical
583   simp [graphSpectrum]
584 
585 /-- Exact geometric multiplicity of an eigenvalue. -/
586 noncomputable def graphEigenMultiplicity
587     (A : Matrix ν ν K) (μ : K) : ℕ :=
588   Module.finrank K (Module.End.eigenspace (Matrix.toLin' A) μ)
589 
590 /--
591 The spectral computation in Paper I: take the maximum exact eigenspace
592 dimension over the finite spectrum, with value zero for the empty graph.
593 -/
594 noncomputable def exactGraphSeedRank (A : Matrix ν ν K) : ℕ :=
595   (graphSpectrum A).sup (graphEigenMultiplicity A)
596 
597 theorem graphEigenMultiplicity_le_exactGraphSeedRank
598     (A : Matrix ν ν K) {μ : K} (hμ : μ ∈ spectrum K A) :
599     graphEigenMultiplicity A μ ≤ exactGraphSeedRank A := by
600   classical
601   exact Finset.le_sup (s := graphSpectrum A)
602     (f := graphEigenMultiplicity A) ((mem_graphSpectrum_iff A μ).2 hμ)
603 
604 /--
605 End-to-end correctness from exact spectral certificates.  The eigenbasis and
606 Lagrange-projector data prove that `R` is the least number of operator seeds;
607 the coverage and dimension certificates prove that the finite-spectrum
608 maximum computed above is exactly that same `R`.
609 -/
610 theorem exactGraphSeedRank_isLeast
611     {Λ : Type w} [Fintype Λ]
612     (A : Matrix ν ν K) (R : ℕ)
613     (multiplicity : Λ → ℕ) (eigenvalue : Λ → K)
614     (eigenbasis :
615       Module.Basis (Σ eigen : Λ, Fin (multiplicity eigen)) K (ν → K))
616     (hEigen : ∀ z : Σ eigen : Λ, Fin (multiplicity eigen),
617       Matrix.toLin' A (eigenbasis z) = eigenvalue z.1 • eigenbasis z)
618     (hbound : ∀ eigen, multiplicity eigen ≤ R)
619     (maxEigen : Λ) (hmax : multiplicity maxEigen = R)
620     (projectorPolynomial : Λ → K[X])
621     (hproject : ∀ (eigen : Λ) (j : Fin (multiplicity eigen)),
622       aeval (Matrix.toLin' A) (projectorPolynomial eigen)
623         (SingleOperatorGraph.paddedBlockSum R multiplicity
624           (fun e j ↦ eigenbasis ⟨e, j⟩)
625           ⟨j.val, lt_of_lt_of_le j.isLt (hbound eigen)⟩) =
626         eigenbasis ⟨eigen, j⟩)
627     (hSpectrumCover : ∀ μ ∈ spectrum K A, ∃ eigen, eigenvalue eigen = μ)
628     (hMultiplicity : ∀ eigen,
629       graphEigenMultiplicity A (eigenvalue eigen) = multiplicity eigen)
630     (hMaxSpectrum : eigenvalue maxEigen ∈ spectrum K A) :
631     IsLeast
632       {r : ℕ | SingleOperatorGraph.SeedNumberLE (Matrix.toLin' A) r}
633       (exactGraphSeedRank A) := by
634   classical
635   have hLeast :
636       IsLeast
637         {r : ℕ | SingleOperatorGraph.SeedNumberLE (Matrix.toLin' A) r} R :=
638     SingleOperatorGraph.graph_seedNumber_eq_max_eigenspace_multiplicity
639       (Matrix.toLin' A) R multiplicity eigenvalue eigenbasis hEigen
640       hbound maxEigen hmax projectorPolynomial hproject
641   have hComputed : exactGraphSeedRank A = R := by
642     apply Nat.le_antisymm
643     · apply Finset.sup_le
644       intro μ hμ
645       have hμSpectrum : μ ∈ spectrum K A :=
646         (mem_graphSpectrum_iff A μ).1 hμ
647       obtain ⟨eigen, heigen⟩ := hSpectrumCover μ hμSpectrum
648       rw [← heigen, hMultiplicity]
649       exact hbound eigen
650     · rw [← hmax, ← hMultiplicity maxEigen]
651       exact Finset.le_sup
652         ((mem_graphSpectrum_iff A (eigenvalue maxEigen)).2 hMaxSpectrum)
653   rw [hComputed]
654   exact hLeast
655 
656 end GraphRank
657 
658 end FiniteAlgorithms
659 
660 end Komorebi
661 
662 #print axioms Komorebi.EllipticCoordinates.My_sq_eq_Mx_cube_sub_Mx
663 #print axioms Komorebi.EllipticCoordinates.unitSeed_generates
664 #print axioms Komorebi.EllipticCoordinates.Mx_pow_My_pow_unitSeed
665 #print axioms Komorebi.FiniteAlgorithms.orbitHull_algorithm_terminates_and_correct
666 #print axioms Komorebi.FiniteAlgorithms.transferBlock_algorithm_correct
667 #print axioms Komorebi.FiniteAlgorithms.exactGraphSeedRank_isLeast

Back to source controls