Posterior geometry.

  1 import Mathlib.LinearAlgebra.Matrix.Rank  2 import Mathlib.LinearAlgebra.Matrix.MvPolynomial  3 import Mathlib.Topology.Instances.Matrix  4 import Mathlib.Topology.Algebra.MvPolynomial  5 import Mathlib.Algebra.MvPolynomial.Equiv  6 import Mathlib.Algebra.MvPolynomial.Funext  7 import Mathlib.MeasureTheory.Constructions.Pi  8 import Mathlib.MeasureTheory.Measure.MutuallySingular  9 import Mathlib.MeasureTheory.Measure.Lebesgue.Basic 10 import Mathlib.MeasureTheory.Measure.Haar.Unique 11 import Mathlib.MeasureTheory.Measure.Lebesgue.EqHaar 12 import Mathlib.Tactic 13  14 /-! 15 Posterior geometry for Singular Bayesian Neural Networks 16 by Mame Diarra Toure and David A. Stephens. 17  18 Constrained support, the measure-zero low-rank set, and singularity. 19 The factor law is arbitrary: no independence or Gaussian assumption is needed. 20 The matrix measure below is the product of Lebesgue measures on its entries. 21 -/ 22  23 open MeasureTheory 24  25 namespace SingularBNN 26  27 abbrev Weight (m n : Nat) := Matrix (Fin m) (Fin n) Real 28  29 noncomputable instance weightMeasureSpace (m n : Nat) : MeasureSpace (Weight m n) := 30   inferInstanceAs (MeasureSpace (Fin m -> Fin n -> Real)) 31  32 instance weightBorelSpace (m n : Nat) : BorelSpace (Weight m n) := 33   inferInstanceAs (BorelSpace (Fin m -> Fin n -> Real)) 34  35 instance weightSecondCountable (m n : Nat) : SecondCountableTopology (Weight m n) := 36   inferInstanceAs (SecondCountableTopology (Fin m -> Fin n -> Real)) 37  38 instance weightLocallyCompact (m n : Nat) : LocallyCompactSpace (Weight m n) := 39   inferInstanceAs (LocallyCompactSpace (Fin m -> Fin n -> Real)) 40  41 instance weightSigmaFinite (m n : Nat) : SigmaFinite (volume : Measure (Weight m n)) := 42   inferInstanceAs (SigmaFinite (volume : Measure (Fin m -> Fin n -> Real))) 43  44 instance weightHaar (m n : Nat) : Measure.IsAddHaarMeasure (volume : Measure (Weight m n)) := by 45   letI := isAddHaarMeasure_volume_pi (Fin n) 46   exact Measure.pi.isAddHaarMeasure (fun _ : Fin m => (volume : Measure (Fin n -> Real))) 47  48 abbrev Factors (m n r : Nat) := Prod (Weight m r) (Weight n r) 49  50 def factorMap {m n r : Nat} (p : Factors m n r) : Weight m n := 51   p.1 * p.2.transpose 52  53 def LowRank (m n r : Nat) : Set (Weight m n) := 54   {W | W.rank <= r} 55  56 noncomputable def inducedPosterior {m n r : Nat} 57     (q : Measure (Factors m n r)) : Measure (Weight m n) := 58   q.map factorMap 59  60 theorem rank_factorMap_le {m n r : Nat} (p : Factors m n r) : 61     (factorMap p).rank <= r := 62   (Matrix.rank_mul_le_left _ _).trans (Matrix.rank_le_width p.1) 63  64 theorem continuous_factorMap {m n r : Nat} : 65     Continuous (factorMap : Factors m n r -> Weight m n) := by 66   unfold factorMap 67   fun_prop 68  69 theorem posterior_lowRank_mass_one {m n r : Nat} 70     (q : Measure (Factors m n r)) [IsProbabilityMeasure q] : 71     inducedPosterior q (LowRank m n r) = 1 := by 72   have hpre : Set.preimage (factorMap (m := m) (n := n) (r := r)) 73       (LowRank m n r) = Set.univ := by 74     ext p 75     simp only [Set.mem_preimage, LowRank, Set.mem_setOf_eq, Set.mem_univ, iff_true] 76     exact rank_factorMap_le p 77   have hmap : inducedPosterior q Set.univ = 1 := by 78     rw [inducedPosterior, Measure.map_apply continuous_factorMap.measurable MeasurableSet.univ] 79     simp 80   apply le_antisymm ((measure_mono (Set.subset_univ _)).trans_eq hmap) 81   have h := show q (Set.preimage factorMap (LowRank m n r)) <= inducedPosterior q (LowRank m n r) from 82     Measure.le_map_apply continuous_factorMap.measurable.aemeasurable (LowRank m n r) 83   simpa only [hpre, measure_univ, inducedPosterior] using h 84  85  86  87  88  89 /-- A nonzero real polynomial in finitely many variables vanishes on a null set. -/ 90 theorem polynomial_zero_volume_fin (k : Nat) (p : MvPolynomial (Fin k) Real) 91     (hp : Not (p = 0)) : 92     volume {x : Fin k -> Real | MvPolynomial.eval x p = 0} = 0 := by 93   classical 94   induction k with 95   | zero => 96     have hn : forall x : Fin 0 -> Real, Not (MvPolynomial.eval x p = 0) := by 97       intro x hx 98       apply hp 99       apply MvPolynomial.funext100       intro y101       have hy : y = x := Subsingleton.elim _ _102       simp [hy, hx]103     have hempty : {x : Fin 0 -> Real | MvPolynomial.eval x p = 0} = {} := by104       ext x105       simp [hn x]106     rw [hempty, measure_empty]107   | succ k ih =>108     let P := MvPolynomial.finSuccEquiv Real k p109     have hP : Not (P = 0) := by110       intro hz111       apply hp112       apply (MvPolynomial.finSuccEquiv Real k).injective113       simpa only [map_zero] using hz114     have hex : exists j : Nat, Not (P.coeff j = 0) := by115       by_contra! h116       exact hP (Polynomial.ext fun i => by simpa using h i)117     let j := Classical.choose hex118     have hj : Not (P.coeff j = 0) := Classical.choose_spec hex119     have htail : Filter.Eventually120         (fun x : Fin k -> Real => Not (MvPolynomial.eval x (P.coeff j) = 0))121         (ae volume) := by122       apply ae_iff.mpr123       simpa only [not_not] using ih (P.coeff j) hj124     have hsec : Filter.Eventually125         (fun x : Fin k -> Real =>126           volume {y : Real | MvPolynomial.eval (Fin.cons y x) p = 0} = 0)127         (ae volume) := by128       filter_upwards [htail] with x hx129       have hspec : Not (Polynomial.map (MvPolynomial.eval x) P = 0) := by130         intro hz131         have hc := congrArg (fun Q : Polynomial Real => Q.coeff j) hz132         simp only [Polynomial.coeff_map, Polynomial.coeff_zero] at hc133         exact hx hc134       simpa only [MvPolynomial.eval_eq_eval_mv_eval', P, Polynomial.IsRoot.def] using135         (Polynomial.finite_setOf_isRoot hspec).measure_zero volume136     have hm : MeasurableSet {z : Prod (Fin k -> Real) Real |137         MvPolynomial.eval (Fin.cons z.2 z.1) p = 0} := by138       exact (isClosed_eq (p.continuous_eval.comp (by fun_prop)) continuous_const).measurableSet139     have hprod : (volume : Measure (Prod (Fin k -> Real) Real))140         {z | MvPolynomial.eval (Fin.cons z.2 z.1) p = 0} = 0 := by141       apply Measure.measure_prod_null_of_ae_null hm142       exact hsec143     have hsplit := volume_preserving_piFinSuccAbove144       (fun _ : Fin (k + 1) => Real) 0145     have hswap : MeasurePreserving Prod.swap146         ((volume : Measure Real).prod (volume : Measure (Fin k -> Real)))147         ((volume : Measure (Fin k -> Real)).prod (volume : Measure Real)) :=148       Measure.measurePreserving_swap149     have h := (hswap.comp hsplit).quasiMeasurePreserving.preimage_null hprod150     convert h using 2151     ext x152     simp only [Set.mem_setOf_eq, Set.mem_preimage, Function.comp_apply,153       MeasurableEquiv.piFinSuccAbove_apply]154     have hcons : Fin.cons ((Fin.insertNthEquiv (fun _ : Fin (k + 1) => Real) 0).symm x).swap.2155         ((Fin.insertNthEquiv (fun _ : Fin (k + 1) => Real) 0).symm x).swap.1 = x := by156       ext i157       refine Fin.cases ?_ (fun j => ?_) i158       next => simp [Fin.insertNthEquiv]159       next => simp [Fin.insertNthEquiv, Fin.tail]160     rw [hcons]161 162 163 /-- Reindexing the variables preserves the null-set conclusion. -/164 theorem polynomial_zero_volume {s : Type} [Fintype s]165     (p : MvPolynomial s Real) (hp : Not (p = 0)) :166     volume {x : s -> Real | MvPolynomial.eval x p = 0} = 0 := by167   classical168   let e := Fintype.equivFin s169   let P := MvPolynomial.renameEquiv Real e p170   have hP : Not (P = 0) := by171     intro hz172     apply hp173     apply (MvPolynomial.renameEquiv Real e).injective174     simpa only [map_zero] using hz175   have hfin := polynomial_zero_volume_fin (Fintype.card s) P hP176   have h := (volume_measurePreserving_piCongrLeft177     (fun _ : Fin (Fintype.card s) => Real) e).quasiMeasurePreserving.preimage_null hfin178   convert h using 2179   ext x180   simp [P, MvPolynomial.renameEquiv_apply, MvPolynomial.eval_rename,181     MeasurableEquiv.piCongrLeft, Equiv.piCongrLeft, Function.comp_def]182 183 /-- Flattening a matrix is a continuous additive equivalence, so it preserves null sets. -/184 theorem matrix_polynomial_zero_volume {m n : Nat}185     (p : MvPolynomial (Prod (Fin m) (Fin n)) Real) (hp : Not (p = 0)) :186     volume {W : Weight m n |187       MvPolynomial.eval (fun ij : Prod (Fin m) (Fin n) => W ij.1 ij.2) p = 0} = 0 := by188   let e : AddEquiv (Weight m n) (Prod (Fin m) (Fin n) -> Real) :=189     { toEquiv := (Equiv.curry (Fin m) (Fin n) Real).symm190       map_add' := by intros; rfl }191   have he : Continuous e := by192     change Continuous (fun W : Weight m n => fun ij : Prod (Fin m) (Fin n) => W ij.1 ij.2)193     fun_prop194   have hei : Continuous e.symm := by195     change Continuous (fun x : Prod (Fin m) (Fin n) -> Real => fun i j => x (i, j))196     fun_prop197   letI : Measure.IsAddHaarMeasure (Measure.map e (volume : Measure (Weight m n))) :=198     e.isAddHaarMeasure_map volume he hei199   have hac := Measure.absolutelyContinuous_isAddHaarMeasure200     (Measure.map e (volume : Measure (Weight m n))) volume201   have hz := hac (polynomial_zero_volume p hp)202   rw [Measure.map_apply he.measurable203     (isClosed_eq p.continuous_eval continuous_const).measurableSet] at hz204   exact hz205 206 noncomputable def minorPolynomial (m n r : Nat) (hm : r < m) (hn : r < n) :207     MvPolynomial (Prod (Fin m) (Fin n)) Real :=208   Matrix.det (Matrix.of (fun i j : Fin (r + 1) => MvPolynomial.X209     (Fin.castLE (Nat.succ_le_of_lt hm) i, Fin.castLE (Nat.succ_le_of_lt hn) j)))210 211 theorem eval_minorPolynomial {m n r : Nat} (hm : r < m) (hn : r < n)212     (W : Weight m n) :213     MvPolynomial.eval (fun ij : Prod (Fin m) (Fin n) => W ij.1 ij.2) (minorPolynomial m n r hm hn) =214       (W.submatrix (Fin.castLE (Nat.succ_le_of_lt hm))215         (Fin.castLE (Nat.succ_le_of_lt hn))).det := by216   unfold minorPolynomial217   rw [(MvPolynomial.eval (fun ij : Prod (Fin m) (Fin n) => W ij.1 ij.2)).map_det]218   congr 1219   ext i j220   simp221 222 theorem minorPolynomial_ne_zero {m n r : Nat} (hm : r < m) (hn : r < n) :223     Not (minorPolynomial m n r hm hn = 0) := by224   classical225   let W : Weight m n := Matrix.of (fun i j => if i.val = j.val then 1 else 0)226   have hsub : W.submatrix (Fin.castLE (Nat.succ_le_of_lt hm))227       (Fin.castLE (Nat.succ_le_of_lt hn)) = (1 : Matrix (Fin (r + 1)) (Fin (r + 1)) Real) := by228     ext i j229     simp [W, Matrix.submatrix_apply, Matrix.one_apply, Fin.ext_iff]230   have he : MvPolynomial.eval (fun ij : Prod (Fin m) (Fin n) => W ij.1 ij.2)231       (minorPolynomial m n r hm hn) = 1 := by232     rw [eval_minorPolynomial, hsub, Matrix.det_one]233   intro hz234   rw [hz, map_zero] at he235   exact zero_ne_one he236 237 theorem lowRank_minor_zero {m n r : Nat} (hm : r < m) (hn : r < n)238     (W : Weight m n) (hW : Membership.mem (LowRank m n r) W) :239     MvPolynomial.eval (fun ij : Prod (Fin m) (Fin n) => W ij.1 ij.2) (minorPolynomial m n r hm hn) = 0 := by240   rw [eval_minorPolynomial]241   by_contra hdet242   let S := W.submatrix (Fin.castLE (Nat.succ_le_of_lt hm))243     (Fin.castLE (Nat.succ_le_of_lt hn))244   have hu : IsUnit S := (Matrix.isUnit_iff_isUnit_det S).mpr245     (isUnit_iff_ne_zero.mpr hdet)246   have hS : S.rank = r + 1 := by247     simpa using Matrix.rank_of_isUnit S hu248   have hle : S.rank <= r :=249     (Matrix.rank_submatrix_le W _ _).trans hW250   omega251 252 /-- The set of real m-by-n matrices of rank at most r has zero Lebesgue measure. -/253 theorem lowRank_volume_zero {m n r : Nat} (hr : r < min m n) :254     volume (LowRank m n r) = 0 := by255   have hm : r < m := lt_of_lt_of_le hr (min_le_left _ _)256   have hn : r < n := lt_of_lt_of_le hr (min_le_right _ _)257   apply measure_mono_null (fun W hW => lowRank_minor_zero hm hn W hW)258   exact matrix_polynomial_zero_volume _ (minorPolynomial_ne_zero hm hn)259 260 /-- Every measure on the factors induces a measure singular to Lebesgue measure. -/261 theorem inducedPosterior_mutuallySingular {m n r : Nat}262     (q : Measure (Factors m n r)) (hr : r < min m n) :263     (inducedPosterior q).MutuallySingular (volume : Measure (Weight m n)) := by264   have hm : r < m := lt_of_lt_of_le hr (min_le_left _ _)265   have hn : r < n := lt_of_lt_of_le hr (min_le_right _ _)266   let Z : Set (Weight m n) := {W |267     MvPolynomial.eval (fun ij : Prod (Fin m) (Fin n) => W ij.1 ij.2) (minorPolynomial m n r hm hn) = 0}268   have hZ : MeasurableSet Z := by269     exact (isClosed_eq270       ((minorPolynomial m n r hm hn).continuous_eval.comp (by fun_prop))271       continuous_const).measurableSet272   have hvol : volume Z = 0 :=273     matrix_polynomial_zero_volume _ (minorPolynomial_ne_zero hm hn)274   have hpre : Set.preimage (factorMap : Factors m n r -> Weight m n) (compl Z) = {} := by275     ext p276     change Not (Membership.mem Z (factorMap p)) <-> False277     have hz : Membership.mem Z (factorMap p) :=278       lowRank_minor_zero hm hn (factorMap p) (rank_factorMap_le p)279     simp only [hz, not_true_eq_false]280   have hq : inducedPosterior q (compl Z) = 0 := by281     change (Measure.map factorMap q) (compl Z) = 0282     have hmap : (Measure.map factorMap q) (compl Z) =283         q (Set.preimage factorMap (compl Z)) :=284       Measure.map_apply continuous_factorMap.measurable hZ.compl285     rw [hmap, hpre, measure_empty]286   refine Exists.intro (compl Z) ?_287   constructor288   next => exact hZ.compl289   constructor290   next => exact hq291   next => simpa only [compl_compl] using hvol292 293 294 end SingularBNN