Explicit spectral-error probabilities for Gaussian rangefinding

import RNLAPhaseB.B3Base
import RNLAPhaseB.Complexity

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

namespace RNLAPhaseB.HMT2011.CorTailSpecErrorGauss

open Classical in
/-- The `j`-th largest element (1-based, with multiplicity) of a finite multiset of reals,
and `0` when `j = 0` or `j` exceeds its cardinality. For a multiset sorted as
`λ₁ ≥ λ₂ ≥ ⋯`, the set of `t` with at least `j` elements `≥ t` is `(-∞, λⱼ]`,
so its supremum is `λⱼ`. -/
noncomputable def multisetKthLargest (s : Multiset ℝ) (j : ℕ) : ℝ :=
  if 1 ≤ j ∧ j ≤ Multiset.card s then
    sSup {t : ℝ | j ≤ Multiset.card (s.filter (fun x => t ≤ x))}
  else 0

/-- The singular value `σⱼ` (1-based) of a real `m × n` matrix `A`: the square root of the
`j`-th largest eigenvalue (with multiplicity) of `Aᵀ A`, i.e. of the `j`-th largest root of
its characteristic polynomial (which splits over `ℝ` since `Aᵀ A` is symmetric positive
semidefinite). Singular values beyond the matrix size are `0` (zero padding), so that
`σ₁ ≥ σ₂ ≥ σ₃ ≥ ⋯` as in the source. -/
noncomputable def singularValue {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (j : ℕ) : ℝ :=
  Real.sqrt (multisetKthLargest (Aᵀ * A).charpoly.roots j)

/-- The tail `∑_{j>k} σⱼ²` of squared singular values of `A` (all `σⱼ` with `j > n` vanish). -/
noncomputable def singularTailSq {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (k : ℕ) : ℝ :=
  ∑ j ∈ Finset.Ioc k n, singularValue A j ^ 2

/-- HMT 2011, Corollary (cor:tail-spec-error-gauss). Hypotheses of Theorem
(thm:avg-frob-error-gauss): `A` a real `m × n` matrix with singular values
`σ₁ ≥ σ₂ ≥ ⋯`, target rank `k ≥ 2`, oversampling `p ≥ 2` with `k + p ≤ min{m,n}`,
`Ω` an `n × (k+p)` standard Gaussian matrix, `Y = AΩ`; further `p ≥ 4`.
Conclusion: the spectral-norm bound
`‖(I - P_Y)A‖ ≤ (1 + 17√(1 + k/p)) σ_{k+1} + (8√(k+p)/(p+1)) (∑_{j>k} σⱼ²)^{1/2}`
fails with probability at most `6 e^{-p}`, and
`‖(I - P_Y)A‖ ≤ (1 + 8√((k+p) p log p)) σ_{k+1} + 3√(k+p) (∑_{j>k} σⱼ²)^{1/2}`
fails with probability at most `6 p^{-p}`. -/
theorem cor_tail_spec_error_gauss
    {Ω : Type*} [MeasurableSpace Ω] (P : Measure Ω)
    {m n k p : ℕ} (A : Matrix (Fin m) (Fin n) ℝ)
    (hk : 2 ≤ k) (hp2 : 2 ≤ p) (hp4 : 4 ≤ p) (hkp : k + p ≤ min m n)
    (Omega : Ω → Matrix (Fin n) (Fin (k + p)) ℝ)
    (hOmega : StandardGaussianTestMatrix P Omega) :
    P {ω | ¬ (‖(1 - orthogonalProjector (A * Omega ω)) * A‖ ≤
        (1 + 17 * Real.sqrt (1 + (k : ℝ) / (p : ℝ))) * singularValue A (k + 1) +
          8 * Real.sqrt ((k : ℝ) + (p : ℝ)) / ((p : ℝ) + 1) *
            Real.sqrt (singularTailSq A k))}
      ≤ ENNReal.ofReal (6 * Real.exp (-(p : ℝ))) ∧
    P {ω | ¬ (‖(1 - orthogonalProjector (A * Omega ω)) * A‖ ≤
        (1 + 8 * Real.sqrt (((k : ℝ) + (p : ℝ)) * (p : ℝ) * Real.log (p : ℝ))) *
            singularValue A (k + 1) +
          3 * Real.sqrt ((k : ℝ) + (p : ℝ)) * Real.sqrt (singularTailSq A k))}
      ≤ ENNReal.ofReal (6 * (p : ℝ) ^ (-(p : ℝ))) := by
  sorry

end RNLAPhaseB.HMT2011.CorTailSpecErrorGauss