Expected spectral error of randomized SVD

import RNLAPhaseB.B3Base
import RNLAPhaseB.Complexity

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

namespace RNLAPhaseB.HMT2011Main

/-- The `j`th singular value (1-based) of a real `m × n` matrix `A`, written `σ_j` in the source.
It is the square root of the `j`th largest eigenvalue of `Aᴴ A`. Mathlib's `eigenvalues₀` lists
these eigenvalues in decreasing order. For `j > n` the value is padded with `0`. When
`j ≤ min m n` this is the usual singular value. -/
noncomputable def singularValue {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (j : ℕ) : ℝ :=
  if h : 1 ≤ j ∧ j - 1 < Fintype.card (Fin n) then
    Real.sqrt ((Matrix.isHermitian_conjTranspose_mul_self A).eigenvalues₀ ⟨j - 1, h.2⟩)
  else 0

/-- Power-scheme sample matrix `Y = (A Aᵀ)^q A Ω` for an `n × 2k` test matrix `Ω`. -/
noncomputable def powerSample {m n l : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (q : ℕ)
    (G : Matrix (Fin n) (Fin l) ℝ) : Matrix (Fin m) (Fin l) ℝ :=
  (A * Aᵀ) ^ q * A * G

/-- The matrix represented by the rank-`2k` factorization `U Σ Vᵀ` returned by Randomized SVD.
The algorithm takes `Q` with orthonormal columns spanning `range Y` and sets `B = Qᵀ A = Ũ Σ Vᵀ`
and `U = Q Ũ`. Then `U Σ Vᵀ = Q Qᵀ A = P_Y A`, and this does not depend on the choice of `Q` or
of the SVD. -/
noncomputable def randomizedSVDApprox {m n l : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (q : ℕ)
    (G : Matrix (Fin n) (Fin l) ℝ) : Matrix (Fin m) (Fin n) ℝ :=
  orthogonalProjector (powerSample A q G) * A

/-- HMT2011 Theorem (thm:main). Let `A` be a real `m × n` matrix, let `q` be the exponent, and let
`2 ≤ k ≤ 0.5 min{m,n}`. Run Randomized SVD with a standard Gaussian `n × 2k` test matrix. Then the
expected spectral-norm error satisfies
`E ‖A - UΣVᵀ‖ ≤ [1 + 4 √(2 min{m,n}/(k-1))]^{1/(2q+1)} σ_{k+1}`. -/
theorem randomizedSVD_expected_spectral_error
    {Ω : Type*} [MeasurableSpace Ω] (P : Measure Ω)
    {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (q k : ℕ)
    (hk2 : 2 ≤ k) (hkmin : (k : ℝ) ≤ 0.5 * ((min m n : ℕ) : ℝ))
    (G : Ω → Matrix (Fin n) (Fin (2 * k)) ℝ)
    (hG : StandardGaussianTestMatrix P G) :
    ∫ ω, ‖A - randomizedSVDApprox A q (G ω)‖ ∂P ≤
      (1 + 4 * Real.sqrt (2 * ((min m n : ℕ) : ℝ) / ((k : ℝ) - 1))) ^
          ((1 : ℝ) / (2 * (q : ℝ) + 1)) * singularValue A (k + 1) := by
  sorry

end RNLAPhaseB.HMT2011Main