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.PowerMethodSpecGauss

/-- The `j`th singular value (1-based) of a real `m × n` matrix, `σ_j(A)`.
It is given by the standard rank characterization (Schmidt–Mirsky/Eckart–Young
in the spectral norm): `σ_j(A) = inf { ‖A - X‖ : rank X < j }`, where `‖·‖` is
the ℓ₂ operator norm. This equals the `j`th largest singular value for
`1 ≤ j ≤ min m n` and is `0` (zero padding) for `j > min m n`. -/
noncomputable def singularValue {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (j : ℕ) : ℝ :=
  sInf {r : ℝ | ∃ X : Matrix (Fin m) (Fin n) ℝ, X.rank < j ∧ r = ‖A - X‖}

/-- The source's matrix `B = (A A^*)^q A` (real case, so `A^* = Aᵀ`). -/
noncomputable def powerMatrix {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (q : ℕ) :
    Matrix (Fin m) (Fin n) ℝ :=
  (A * Aᵀ) ^ q * A

/-- The source's sample matrix `Z = B Ω` with `B = (A A^*)^q A`. -/
noncomputable def sampleMatrix {m n ℓ : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (q : ℕ)
    (Omega : Matrix (Fin n) (Fin ℓ) ℝ) : Matrix (Fin m) (Fin ℓ) ℝ :=
  powerMatrix A q * Omega

/-- HMT2011, Corollary cor:power-method-spec-gauss. Hypotheses of
Theorem thm:avg-frob-error-gauss: `A` real `m × n`, target rank `k ≥ 2`,
oversampling `p ≥ 2`, `k + p ≤ min{m,n}`, `Ω` an `n × (k+p)` standard Gaussian
matrix. With `B = (AA^*)^q A` (`q` a nonnegative integer) and `Z = BΩ`,
`E ‖(I - P_Z) A‖ ≤ [(1 + √(k/(p-1))) σ_{k+1}^{2q+1}
   + (e √(k+p)/p) (∑_{j>k} σ_j^{2(2q+1)})^{1/2}]^{1/(2q+1)}`,
with `‖·‖` the spectral (ℓ₂ operator) norm and `σ_j` the singular values of `A`. -/
theorem power_method_spec_gauss
    {Ω : Type*} [MeasurableSpace Ω] (P : Measure Ω)
    {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) (k p q : ℕ)
    (hk : 2 ≤ k) (hp : 2 ≤ p) (hkp : k + p ≤ min m n)
    (G : Ω → Matrix (Fin n) (Fin (k + p)) ℝ)
    (hG : StandardGaussianTestMatrix P G) :
    ∫ ω, ‖(1 - orthogonalProjector (sampleMatrix A q (G ω))) * A‖ ∂P ≤
      ((1 + Real.sqrt ((k : ℝ) / ((p : ℝ) - 1))) * singularValue A (k + 1) ^ (2 * q + 1)
        + Real.exp 1 * Real.sqrt ((k : ℝ) + (p : ℝ)) / (p : ℝ) *
          Real.sqrt (∑ j ∈ Finset.Ioc k (min m n), singularValue A j ^ (2 * (2 * q + 1))))
        ^ ((1 : ℝ) / (2 * (q : ℝ) + 1)) := by
  sorry

end RNLAPhaseB.HMT2011.PowerMethodSpecGauss