CirculantFourier.lean

1 module
2 public import Mathlib.Analysis.Fourier.ZMod
3 public import Mathlib.LinearAlgebra.Matrix.Circulant
4 public import Mathlib.Tactic
5 
6 /-!
7 Author: Komorebi P. Paper VIII 4.3.
8 The actual unitary-normalized Fourier matrix diagonalizes the actual circulant.
9 The standard character is explicitly linked to the complex exponential.
10 -/
11 @[expose] public section
12 set_option autoImplicit false
13 namespace KomorebiPaperVIII
14 open Matrix ZMod
15 open scoped BigOperators
16 variable {n : ℕ} [NeZero n]
17 
18 noncomputable def fourierScale (n : ℕ) : ℂ := ((Real.sqrt (n:ℝ))⁻¹ : ℝ)
19 noncomputable def fourierMatrix (n : ℕ) [NeZero n] : Matrix (ZMod n) (ZMod n) ℂ :=
20   fun i j => fourierScale n * ZMod.stdAddChar (-(i*j))
21 
22 theorem fourier_scale_normalization :
23     fourierScale n * fourierScale n * (n:ℂ) = 1 := by
24   have hn : (0:ℝ) < n := Nat.cast_pos.mpr (NeZero.pos n)
25   have hsq : Real.sqrt (n:ℝ) * Real.sqrt (n:ℝ) = n := Real.mul_self_sqrt hn.le
26   have hs : Real.sqrt (n:ℝ) ≠ 0 := (Real.sqrt_pos.mpr hn).ne'
27   have h : (Real.sqrt (n:ℝ))⁻¹ * (Real.sqrt (n:ℝ))⁻¹ * (n:ℝ) = 1 := by
28     field_simp [hs]
29     nlinarith [hsq]
30   unfold fourierScale
31   exact_mod_cast h
32 
33 theorem fourier_scale_sqrt : (Real.sqrt (n:ℝ) : ℂ) * fourierScale n = 1 := by
34   have hs : Real.sqrt (n:ℝ) ≠ 0 := (Real.sqrt_pos.mpr (Nat.cast_pos.mpr (NeZero.pos n))).ne'
35   unfold fourierScale
36   norm_cast
37   exact mul_inv_cancel₀ hs
38 
39 theorem fourier_matrix_exponential (i j : ZMod n) :
40     fourierMatrix n i j = fourierScale n *
41       Complex.exp (-(2 * Real.pi * Complex.I * (i.val:ℂ) * (j.val:ℂ) / (n:ℂ))) := by
42   unfold fourierMatrix
43   have heq : ((-((i.val:ℤ)*(j.val:ℤ)) : ℤ) : ZMod n) = -(i*j) := by simp
44   rw [← heq, ZMod.stdAddChar_coe]
45   congr 2
46   push_cast
47   ring
48 
49 theorem fourier_action (x : ZMod n → ℂ) :
50     fourierMatrix n *ᵥ x = fourierScale n • ZMod.dft x := by
51   ext k
52   simp only [fourierMatrix, Matrix.mulVec, dotProduct, Pi.smul_apply,
53     ZMod.dft_apply, smul_eq_mul, Finset.mul_sum]
54   apply Finset.sum_congr rfl
55   intro i _
56   rw [mul_comm i k]
57   ring
58 
59 theorem fourier_adjoint_entry (i j : ZMod n) :
60     (fourierMatrix n)ᴴ i j = fourierScale n * ZMod.stdAddChar (i*j) := by
61   simp only [Matrix.conjTranspose_apply, fourierMatrix, star_mul, fourierScale,
62     Complex.star_def, Complex.conj_ofReal]
63   rw [← AddChar.map_neg_eq_conj, neg_neg, mul_comm j i]
64   ring
65 
66 theorem fourier_adjoint_action (x : ZMod n → ℂ) :
67     (fourierMatrix n)ᴴ *ᵥ x = (fourierScale n * (n:ℂ)) • ZMod.dft.symm x := by
68   have hn : (n:ℂ) ≠ 0 := by exact_mod_cast NeZero.ne n
69   ext k
70   simp only [Matrix.mulVec, dotProduct, fourier_adjoint_entry, Pi.smul_apply,
71     ZMod.invDFT_apply, smul_eq_mul, Finset.mul_sum]
72   calc
73     _ = fourierScale n * ∑ j, ZMod.stdAddChar (j*k) * x j := by
74       rw [Finset.mul_sum]
75       apply Finset.sum_congr rfl
76       intro j _
77       rw [mul_comm k j]
78       ring
79     _ = _ := by
80       field_simp
81       simp only [Finset.mul_sum, mul_assoc]
82 
83 theorem fourier_unitary :
84     (fourierMatrix n)ᴴ * fourierMatrix n = 1 ∧
85     fourierMatrix n * (fourierMatrix n)ᴴ = 1 := by
86   have hnorm := fourier_scale_normalization (n := n)
87   constructor
88   · apply Matrix.ext_iff_mulVec.mpr
89     intro x
90     rw [← Matrix.mulVec_mulVec, fourier_action, fourier_adjoint_action,
91       _root_.map_smul, LinearEquiv.symm_apply_apply, smul_smul, Matrix.one_mulVec]
92     have hh : (fourierScale n * (n:ℂ)) * fourierScale n = 1 := by
93       calc
94         _ = (fourierScale n * fourierScale n) * (n:ℂ) := by ring
95         _ = 1 := hnorm
96     rw [hh, one_smul]
97   · apply Matrix.ext_iff_mulVec.mpr
98     intro x
99     rw [← Matrix.mulVec_mulVec, fourier_adjoint_action, fourier_action,
100       _root_.map_smul, LinearEquiv.apply_symm_apply, smul_smul, Matrix.one_mulVec]
101     have hh : fourierScale n * (fourierScale n * (n:ℂ)) = 1 := by rw [← mul_assoc, hnorm]
102     rw [hh, one_smul]
103 
104 /-- The convolution theorem uses an explicit bijective change of finite index. -/
105 theorem dft_circulant (c x : ZMod n → ℂ) :
106     ZMod.dft (Matrix.circulant c *ᵥ x) = fun k => ZMod.dft c k * ZMod.dft x k := by
107   funext k
108   simp only [ZMod.dft_apply, Matrix.mulVec, dotProduct, Matrix.circulant_apply,
109     smul_eq_mul, Finset.mul_sum, Finset.sum_mul]
110   rw [Finset.sum_comm]
111   apply Finset.sum_congr rfl
112   intro j _
113   refine Fintype.sum_equiv (Equiv.subRight j) _ _ ?_
114   intro a
115   simp only [Equiv.subRight_apply]
116   have hchar : ZMod.stdAddChar (-((a-j)*k)) * ZMod.stdAddChar (-(j*k)) =
117       ZMod.stdAddChar (-(a*k)) := by
118     rw [← AddChar.map_add_eq_mul]
119     congr 1
120     ring
121   calc
122     _ = (ZMod.stdAddChar (-((a-j)*k)) * ZMod.stdAddChar (-(j*k))) * (c (a-j)*x j) := by
123       rw [hchar]
124     _ = _ := by ring
125 
126 
127 /-- The diagonal intertwining identity precedes reconstruction. -/
128 theorem fourier_circulant_intertwine (c : ZMod n → ℂ) :
129     fourierMatrix n * Matrix.circulant c = Matrix.diagonal (ZMod.dft c) * fourierMatrix n := by
130   apply Matrix.ext_iff_mulVec.mpr
131   intro x
132   rw [← Matrix.mulVec_mulVec, ← Matrix.mulVec_mulVec, fourier_action, dft_circulant,
133     fourier_action]
134   ext k
135   simp [Matrix.mulVec_diagonal, Pi.smul_apply]
136   ring
137 
138 theorem circulant_spectral_dft (c : ZMod n → ℂ) :
139     Matrix.circulant c = (fourierMatrix n)ᴴ * Matrix.diagonal (ZMod.dft c) * fourierMatrix n := by
140   calc
141     _ = ((fourierMatrix n)ᴴ * fourierMatrix n) * Matrix.circulant c := by
142       rw [(fourier_unitary (n := n)).1, Matrix.one_mul]
143     _ = _ := by rw [Matrix.mul_assoc, fourier_circulant_intertwine, ← Matrix.mul_assoc]
144 
145 theorem normalized_fourier_eigenvalues (c : ZMod n → ℂ) :
146     (Real.sqrt (n:ℝ) : ℂ) • (fourierMatrix n *ᵥ c) = ZMod.dft c := by
147   rw [fourier_action, smul_smul, fourier_scale_sqrt, one_smul]
148 
149 /-- Exactly the normalized Fourier-matrix formula in the manuscript. -/
150 theorem circulant_spectral_representation (c : ZMod n → ℂ) :
151     Matrix.circulant c = (fourierMatrix n)ᴴ *
152       Matrix.diagonal ((Real.sqrt (n:ℝ) : ℂ) • (fourierMatrix n *ᵥ c)) * fourierMatrix n := by
153   rw [normalized_fourier_eigenvalues, circulant_spectral_dft]
154 
155 theorem circulant_fourier_application (c x : ZMod n → ℂ) :
156     Matrix.circulant c *ᵥ x = (fourierMatrix n)ᴴ *ᵥ
157       (fun k => ZMod.dft c k * (fourierMatrix n *ᵥ x) k) := by
158   rw [circulant_spectral_dft, ← Matrix.mulVec_mulVec, ← Matrix.mulVec_mulVec]
159   congr 1
160   funext k
161   exact Matrix.mulVec_diagonal _ _ _
162 
163 end KomorebiPaperVIII
164 
165 -- FINAL AXIOM AUDITS
166 #print axioms KomorebiPaperVIII.fourier_scale_normalization
167 #print axioms KomorebiPaperVIII.fourier_scale_sqrt
168 #print axioms KomorebiPaperVIII.fourier_matrix_exponential
169 #print axioms KomorebiPaperVIII.fourier_action
170 #print axioms KomorebiPaperVIII.fourier_adjoint_entry
171 #print axioms KomorebiPaperVIII.fourier_adjoint_action
172 #print axioms KomorebiPaperVIII.fourier_unitary
173 #print axioms KomorebiPaperVIII.dft_circulant
174 #print axioms KomorebiPaperVIII.fourier_circulant_intertwine
175 #print axioms KomorebiPaperVIII.circulant_spectral_dft
176 #print axioms KomorebiPaperVIII.normalized_fourier_eigenvalues
177 #print axioms KomorebiPaperVIII.circulant_spectral_representation
178 #print axioms KomorebiPaperVIII.circulant_fourier_application
179 
180 #print axioms KomorebiPaperVIII.fourierMatrix

Back to source controls