import RNLAPhaseB.B3Base
import RNLAPhaseB.Complexity

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

namespace RNLAPhaseB.MT2020RandKrylov

/-- The eigenvalues of a real symmetric matrix in weakly decreasing order,
`λ₁ ≥ λ₂ ≥ ⋯ ≥ λₙ`, indexed from `0` (so index `0` is `λ₁` and `n-1` is `λₙ`). -/
noncomputable def eigDesc {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} (hA : A.IsHermitian)
    (i : Fin n) : ℝ :=
  hA.eigenvalues₀ (Fin.cast (Fintype.card_fin n).symm i)

/-- The randomized Krylov maximum-eigenvalue estimate after `q` iterations (source
span bd946d348162):
`ξ_q = max_{u ∈ K_{q+1}} u*Au/u*u = max_{deg p ≤ q} ω* A p²(A) ω / ‖p(A) ω‖²`,
the maximum over real polynomials `p` of degree at most `q` with `p(A) ω ≠ 0`, where
`p(A)` is the spectral function induced by `p` (equal to polynomial evaluation at `A`). -/
noncomputable def krylovEstimate {n : ℕ} (A : Matrix (Fin n) (Fin n) ℝ) (w : Fin n → ℝ)
    (q : ℕ) : ℝ :=
  sSup {r : ℝ | ∃ p : Polynomial ℝ, p.natDegree ≤ q ∧
    (Polynomial.aeval A p) *ᵥ w ≠ 0 ∧
    r = (w ⬝ᵥ ((A * (Polynomial.aeval A p) ^ 2) *ᵥ w)) /
      (l2Norm ((Polynomial.aeval A p) *ᵥ w)) ^ 2}

/-- The error relative to the spectral range (source span 452b9e05bbc7):
`f_q = (λ₁ - ξ_q) / (λ₁ - λₙ)`. -/
noncomputable def krylovError {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ} (hA : A.IsHermitian)
    (hn : 0 < n) (w : Fin n → ℝ) (q : ℕ) : ℝ :=
  (eigDesc hA ⟨0, hn⟩ - krylovEstimate A w q) /
    (eigDesc hA ⟨0, hn⟩ - eigDesc hA ⟨n - 1, by omega⟩)

/-- The relative spectral gap `γ = (λ₁ - λ₂) / λ₁` (requires a second eigenvalue). -/
noncomputable def relativeSpectralGap {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ}
    (hA : A.IsHermitian) (h2 : 1 < n) : ℝ :=
  (eigDesc hA ⟨0, by omega⟩ - eigDesc hA ⟨1, h2⟩) / eigDesc hA ⟨0, by omega⟩

/-- Randomized Krylov method for the maximum eigenvalue of a real psd matrix
(mt2020 thm:rand-krylov). The test vector `ω ~ N(0, I)` is the single column of a
standard Gaussian `n × 1` test matrix. Lean conventions: if `λ₁ = λₙ` then `f_q` is
`0/0 = 0`; if `λ₁ = 0` then `γ = 0/0 = 0`; the source does not treat these cases and
no hypothesis excluding them is added. -/
theorem rand_krylov_expected_error
    {Ω : Type*} [MeasurableSpace Ω] (P : Measure Ω)
    (n : ℕ) (hn : 0 < n)
    (A : Matrix (Fin n) (Fin n) ℝ) (hA : A.PosSemidef)
    (G : Ω → Matrix (Fin n) (Fin 1) ℝ) (hG : RNLAPhaseB.StandardGaussianTestMatrix P G) :
    (∀ q : ℕ, 4 ≤ q →
      ∫ ω, krylovError hA.1 hn (fun i => G ω i 0) q ∂P ≤
        2.575 * (Real.log (n : ℝ) / ((q : ℝ) - 1)) ^ 2) ∧
    (∀ h2 : 1 < n, ∀ q : ℕ, 1 ≤ q →
      ∫ ω, krylovError hA.1 hn (fun i => G ω i 0) q ∂P ≤
        2.589 * Real.sqrt (n : ℝ) *
          Real.exp (-2 * ((q : ℝ) - 1) * Real.sqrt (relativeSpectralGap hA.1 h2))) := by
  sorry

end RNLAPhaseB.MT2020RandKrylov