Pseudoinverse stability of the sampled leading subspace

import RNLAPhaseB.B3Base
import RNLAPhaseB.Complexity

open MeasureTheory
open scoped Matrix
open scoped Matrix.Norms.L2Operator
open scoped BigOperators
open scoped ENNReal

namespace RNLAPhaseB.DM2017.SampleLem40pfn1

/-- The normalized Hadamard transform matrix `H` of order `n` (Sylvester construction,
meaningful for `n = 2^p`): `H i j = (-1)^(popcount (i AND j)) / sqrt n`. -/
noncomputable def normalizedHadamard (n : ℕ) : Matrix (Fin n) (Fin n) ℝ :=
  fun i j => (-1 : ℝ) ^ ((i.val &&& j.val).bits.count true) / Real.sqrt n

/-- The sampling-and-rescaling matrix `S ∈ ℝ^{n × c}` of RandLowRank: its `t`-th column is
`sqrt(n/c) e_{idx t}`, where `idx t` is the integer selected in trial `t`. -/
noncomputable def samplingRescalingMatrix (n c : ℕ) (idx : Fin c → Fin n) :
    Matrix (Fin n) (Fin c) ℝ :=
  fun i t => if idx t = i then Real.sqrt ((n : ℝ) / c) else 0

/-- The diagonal entries of `D` are independent, each `+1` or `-1` with probability `1/2`. -/
def RademacherDiagonal {Ω : Type*} [MeasurableSpace Ω] (P : Measure Ω) {n : ℕ}
    (d : Ω → Fin n → ℝ) : Prop :=
  (∀ i, Measurable (fun ω => d ω i)) ∧
  (∀ i, P {ω | d ω i = 1} = 1 / 2 ∧ P {ω | d ω i = -1} = 1 / 2) ∧
  ProbabilityTheory.iIndepFun (fun i ω => d ω i) P

/-- `c` i.i.d. trials, with replacement, each selecting an integer of `{1,…,n}` uniformly. -/
def UniformIndexSampling {Ω : Type*} [MeasurableSpace Ω] (P : Measure Ω) {n c : ℕ}
    (idx : Ω → Fin c → Fin n) : Prop :=
  (∀ t, Measurable (fun ω => idx ω t)) ∧
  (∀ t j, P {ω | idx ω t = j} = 1 / (n : ℝ≥0∞)) ∧
  ProbabilityTheory.iIndepFun (fun t ω => idx ω t) P

/-- Eqn. (eqn:lem:HU_eqn2new): `‖(H D V_k)_{i*}‖_2^2 ≤ 2 k ln(40 n k) / n` for all `i`. -/
def HUEqn2newHolds {n k : ℕ} (H : Matrix (Fin n) (Fin n) ℝ) (dvec : Fin n → ℝ)
    (V : Matrix (Fin n) (Fin k) ℝ) : Prop :=
  ∀ i : Fin n, l2Norm ((H * Matrix.diagonal dvec * V) i) ^ 2 ≤
    2 * (k : ℝ) * Real.log (40 * (n : ℝ) * (k : ℝ)) / (n : ℝ)

/-- The Moore-Penrose conditions: `X` is the pseudoinverse `M^†` of `M`. -/
def IsMoorePenroseInverse {k c : ℕ} (M : Matrix (Fin k) (Fin c) ℝ)
    (X : Matrix (Fin c) (Fin k) ℝ) : Prop :=
  M * X * M = M ∧ X * M * X = X ∧ (M * X)ᵀ = M * X ∧ (X * M)ᵀ = X * M

/-- The matrix `V_k^T D H S`. -/
noncomputable def sketchedMatrix {n k c : ℕ} (V : Matrix (Fin n) (Fin k) ℝ)
    (dvec : Fin n → ℝ) (idx : Fin c → Fin n) : Matrix (Fin k) (Fin c) ℝ :=
  Vᵀ * Matrix.diagonal dvec * normalizedHadamard n * samplingRescalingMatrix n c idx

/-- Lemma lem:sample_lem40pfn1: conditioned on Eqn. (eqn:lem:HU_eqn2new) holding, if
`c ≥ 192 k ln(40nk)/ε^2 · ln(192 sqrt 20 k ln(40nk)/ε^2)`, then with probability at least
.95, `‖(V_k^T D H S)^† - (V_k^T D H S)^T‖_2^2 ≤ 2 ε^2`. -/
theorem sample_lem40pfn1
    {Ω : Type*} [MeasurableSpace Ω] (P : Measure Ω) [IsProbabilityMeasure P]
    (p n k c : ℕ) (hn : n = 2 ^ p)
    (V : Matrix (Fin n) (Fin k) ℝ) (hV : OrthonormalColumns V)
    (ε : ℝ) (hε0 : 0 < ε) (hε1 : ε < 1 / 2)
    (d : Ω → Fin n → ℝ) (idx : Ω → Fin c → Fin n)
    (hd : RademacherDiagonal P d) (hS : UniformIndexSampling P idx)
    (hind : ProbabilityTheory.IndepFun d idx P)
    (hc : (c : ℝ) ≥ 192 * (k : ℝ) * Real.log (40 * (n : ℝ) * (k : ℝ)) / ε ^ 2 *
      Real.log (192 * Real.sqrt 20 * (k : ℝ) * Real.log (40 * (n : ℝ) * (k : ℝ)) / ε ^ 2)) :
    (95 / 100 : ℝ≥0∞) ≤ ProbabilityTheory.cond P
      {ω | HUEqn2newHolds (normalizedHadamard n) (d ω) V}
      {ω | ∀ X : Matrix (Fin c) (Fin k) ℝ,
        IsMoorePenroseInverse (sketchedMatrix V (d ω) (idx ω)) X →
          ‖X - (sketchedMatrix V (d ω) (idx ω))ᵀ‖ ^ 2 ≤ 2 * ε ^ 2} := by
  sorry

end RNLAPhaseB.DM2017.SampleLem40pfn1