Expected eigenvalue error of randomized power iteration

import RNLAPhaseB.B3Base
import RNLAPhaseB.Complexity

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

namespace RNLAPhaseB.MT2020RandPower

/-- The eigenvalues of a real symmetric `n × n` matrix, listed in weakly decreasing order
`λ₁ ≥ λ₂ ≥ … ≥ λₙ` and counted with multiplicity (source convention, book.md:192). -/
noncomputable def sortedEigenvalues {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ}
    (hA : A.IsHermitian) : List ℝ :=
  (List.ofFn hA.eigenvalues).mergeSort (fun a b => decide (b ≤ a))

/-- `λ_j(A)` with the source's 1-based indexing: the `j`-th largest eigenvalue.
An index beyond `n` is padded with `0`. -/
noncomputable def eigenvalueDesc {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ}
    (hA : A.IsHermitian) (j : ℕ) : ℝ :=
  (sortedEigenvalues hA).getD (j - 1) 0

/-- The power-method iterates: `y₀ = ω / ‖ω‖` and `y_q = A y_{q-1} / ‖A y_{q-1}‖` for `q ≥ 1`. -/
noncomputable def powerIterate {n : ℕ} (A : Matrix (Fin n) (Fin n) ℝ) (ω : Fin n → ℝ) :
    ℕ → (Fin n → ℝ)
  | 0 => (l2Norm ω)⁻¹ • ω
  | q + 1 => (l2Norm (A *ᵥ powerIterate A ω q))⁻¹ • (A *ᵥ powerIterate A ω q)

/-- The eigenvalue estimate `ξ_q = y_q^* A y_q` after `q` iterations. -/
noncomputable def eigenvalueEstimate {n : ℕ} (A : Matrix (Fin n) (Fin n) ℝ)
    (ω : Fin n → ℝ) (q : ℕ) : ℝ :=
  powerIterate A ω q ⬝ᵥ (A *ᵥ powerIterate A ω q)

/-- The relative error `e_q = (λ₁ - ξ_q) / λ₁`. -/
noncomputable def relativeError {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ}
    (hA : A.IsHermitian) (ω : Fin n → ℝ) (q : ℕ) : ℝ :=
  (eigenvalueDesc hA 1 - eigenvalueEstimate A ω q) / eigenvalueDesc hA 1

/-- The relative spectral gap `γ = (λ₁ - λ₂) / λ₁`. -/
noncomputable def relativeSpectralGap {n : ℕ} {A : Matrix (Fin n) (Fin n) ℝ}
    (hA : A.IsHermitian) : ℝ :=
  (eigenvalueDesc hA 1 - eigenvalueDesc hA 2) / eigenvalueDesc hA 1

/-- mt2020 thm:rand-power. `A` is a real psd `n × n` matrix; the test vector
`ω ~ N(0, I)` is the single column of an `n × 1` standard Gaussian test matrix.
The expected relative error of the power-method eigenvalue estimate after `q`
iterations satisfies `E e_q ≤ 0.871 log n / (q - 1)` for `q ≥ 2` and
`E e_q ≤ 1.254 √n γ e^{-qγ}` for `q ≥ 1`. -/
theorem rand_power {Ω : Type*} [MeasurableSpace Ω] (P : Measure Ω)
    (n : ℕ) (A : Matrix (Fin n) (Fin n) ℝ) (hA : A.PosSemidef)
    (Omega : Ω → Matrix (Fin n) (Fin 1) ℝ)
    (hOmega : StandardGaussianTestMatrix P Omega) :
    (∀ q : ℕ, 2 ≤ q →
      ∫ w, relativeError hA.1 (fun i => Omega w i 0) q ∂P
        ≤ 0.871 * (Real.log (n : ℝ) / ((q : ℝ) - 1))) ∧
    (∀ q : ℕ, 1 ≤ q →
      ∫ w, relativeError hA.1 (fun i => Omega w i 0) q ∂P
        ≤ 1.254 * Real.sqrt (n : ℝ) * relativeSpectralGap hA.1 *
            Real.exp (-((q : ℝ) * relativeSpectralGap hA.1))) := by
  sorry

end RNLAPhaseB.MT2020RandPower