Expected error of single-view matrix approximation

import RNLAPhaseB.B3Base
import RNLAPhaseB.Complexity

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

namespace RNLAPhaseB.MT2020SingleView

/-- Squared Frobenius norm of a complex matrix: the sum of the squared moduli of its entries. -/
noncomputable def frobSq {a b : ℕ} (M : Matrix (Fin a) (Fin b) ℂ) : ℝ :=
  ∑ i, ∑ j, ‖M i j‖ ^ 2

/-- The `j`-th largest value (1-based) of a finite real family, counted with multiplicity.
It is zero-padded: it equals `0` for `j` larger than the family size. -/
noncomputable def kthLargest {d : ℕ} (f : Fin d → ℝ) (j : ℕ) : ℝ := by
  classical
  exact sSup {x : ℝ | j ≤ (Finset.univ.filter (fun i => x ≤ f i)).card}

/-- `σ_j`: the `j`-th largest singular value of `A`, 1-based and zero-padded. It is computed
from `A` as square roots of the eigenvalues of `Aᴴ A`. -/
noncomputable def singularValue {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℂ) (j : ℕ) : ℝ :=
  kthLargest (fun i => Real.sqrt ((Matrix.isHermitian_conjTranspose_mul_self A).eigenvalues i)) j

/-- The four Moore-Penrose conditions characterizing `X = M^†`. -/
def IsMoorePenroseInverse {a b : ℕ} (M : Matrix (Fin a) (Fin b) ℂ)
    (X : Matrix (Fin b) (Fin a) ℂ) : Prop :=
  M * X * M = M ∧ X * M * X = X ∧ (M * X)ᴴ = M * X ∧ (X * M)ᴴ = X * M

/-- The Moore-Penrose pseudoinverse `M^†`. It exists and is unique, and is selected by choice. -/
noncomputable def pinv {a b : ℕ} (M : Matrix (Fin a) (Fin b) ℂ) : Matrix (Fin b) (Fin a) ℂ := by
  classical
  exact if h : ∃ X, IsMoorePenroseInverse M X then Classical.choose h else 0

/-- `[Q,~] = qr_econ(Y)`: `Q` has `min(a,b)` orthonormal columns. Also `Y = Q R` for some
upper-triangular `R` of size `min(a,b) × b`, as in economy-size QR. -/
def IsEconQR {a b : ℕ} (Y : Matrix (Fin a) (Fin b) ℂ)
    (Q : Matrix (Fin a) (Fin (min a b)) ℂ) : Prop :=
  Qᴴ * Q = 1 ∧ ∃ R : Matrix (Fin (min a b)) (Fin b) ℂ,
    (∀ (i : Fin (min a b)) (j : Fin b), (j : ℕ) < (i : ℕ) → R i j = 0) ∧ Y = Q * R

/-- Coordinates of the four test matrices: an entry index and a flag (false: real part,
true: imaginary part). -/
abbrev TestCoord (m n ℓ s : ℕ) :=
  ((Fin m × Fin ℓ) ⊕ (Fin n × Fin ℓ) ⊕ (Fin m × Fin s) ⊕ (Fin n × Fin s)) × Bool

/-- Real and imaginary parts of all entries of `Υ, Ω, Φ, Ψ` as real random variables. -/
def testEntries {Ω : Type*} {m n ℓ s : ℕ}
    (Ups : Ω → Matrix (Fin m) (Fin ℓ) ℂ) (Om : Ω → Matrix (Fin n) (Fin ℓ) ℂ)
    (Phi : Ω → Matrix (Fin m) (Fin s) ℂ) (Psi : Ω → Matrix (Fin n) (Fin s) ℂ)
    (t : TestCoord m n ℓ s) (ω : Ω) : ℝ :=
  let z : ℂ := match t.1 with
    | Sum.inl p => Ups ω p.1 p.2
    | Sum.inr (Sum.inl p) => Om ω p.1 p.2
    | Sum.inr (Sum.inr (Sum.inl p)) => Phi ω p.1 p.2
    | Sum.inr (Sum.inr (Sum.inr p)) => Psi ω p.1 p.2
  if t.2 then z.im else z.re

/-- The test matrices `Υ, Ω, Φ, Ψ` are drawn from the standard normal distribution over ℂ.
All entries are independent standard circular complex Gaussians. Real and imaginary parts
are independent `N(0, 1/2)`, so each entry has variance one. -/
structure IndependentStandardComplexGaussianTestMatrices {Ω : Type*} [MeasurableSpace Ω]
    {m n ℓ s : ℕ} (μ : Measure Ω)
    (Ups : Ω → Matrix (Fin m) (Fin ℓ) ℂ) (Om : Ω → Matrix (Fin n) (Fin ℓ) ℂ)
    (Phi : Ω → Matrix (Fin m) (Fin s) ℂ) (Psi : Ω → Matrix (Fin n) (Fin s) ℂ) : Prop where
  isProbabilityMeasure : IsProbabilityMeasure μ
  coord_hasLaw : ∀ t, ProbabilityTheory.HasLaw (testEntries Ups Om Phi Psi t)
    (ProbabilityTheory.gaussianReal 0 (1 / 2 : NNReal)) μ
  coords_iIndep : ProbabilityTheory.iIndepFun (testEntries Ups Om Phi Psi) μ

/-- Line 4: `[P,~] = qr_econ(X)` with `X = Aᴴ Υ`. -/
noncomputable def singleViewP {m n ℓ : ℕ}
    (qrP : Matrix (Fin n) (Fin ℓ) ℂ → Matrix (Fin n) (Fin (min n ℓ)) ℂ)
    (A : Matrix (Fin m) (Fin n) ℂ) (Ups : Matrix (Fin m) (Fin ℓ) ℂ) :
    Matrix (Fin n) (Fin (min n ℓ)) ℂ :=
  qrP (Aᴴ * Ups)

/-- Line 4: `[Q,~] = qr_econ(Y)` with `Y = A Ω`. -/
noncomputable def singleViewQ {m n ℓ : ℕ}
    (qrQ : Matrix (Fin m) (Fin ℓ) ℂ → Matrix (Fin m) (Fin (min m ℓ)) ℂ)
    (A : Matrix (Fin m) (Fin n) ℂ) (Om : Matrix (Fin n) (Fin ℓ) ℂ) :
    Matrix (Fin m) (Fin (min m ℓ)) ℂ :=
  qrQ (A * Om)

/-- Line 5: `C = (Φᴴ Q)^† Z (Pᴴ Ψ)^†` with `Z = Φᴴ A Ψ`. -/
noncomputable def singleViewC {m n ℓ s : ℕ}
    (A : Matrix (Fin m) (Fin n) ℂ) (Phi : Matrix (Fin m) (Fin s) ℂ)
    (Psi : Matrix (Fin n) (Fin s) ℂ)
    (P : Matrix (Fin n) (Fin (min n ℓ)) ℂ) (Q : Matrix (Fin m) (Fin (min m ℓ)) ℂ) :
    Matrix (Fin (min m ℓ)) (Fin (min n ℓ)) ℂ :=
  pinv (Phiᴴ * Q) * (Phiᴴ * A * Psi) * pinv (Pᴴ * Psi)

/-- mt2020 thm:singleview (complex field). Algorithm alg:SingleViewSVD is run on
`A ∈ ℂ^{m×n}` with sampling sizes `ℓ, s`, where `s ≥ 2ℓ`. The test matrices `Υ, Ω, Φ, Ψ` are
independent standard normal complex matrices. `qr_econ` is any economy-QR routine, a fixed
function of its argument. Then
`E ‖A - Q C Pᴴ‖_F^2 ≤ s/(s-ℓ) · min_{k<ℓ} (ℓ+k)/(ℓ-k) ∑_{j=k+1}^{min(m,n)} σ_j^2`.
The expectation is a lower Lebesgue integral. The hypothesis `0 < ℓ` is the domain on which
`min_{k<ℓ}` is defined. The informal remark about the real field is not formalized. -/
theorem singleView_expected_frobenius_error
    {Ω : Type*} [MeasurableSpace Ω] (μ : Measure Ω)
    (m n ℓ s : ℕ) (A : Matrix (Fin m) (Fin n) ℂ)
    (hℓ : 0 < ℓ) (hs : 2 * ℓ ≤ s)
    (Ups : Ω → Matrix (Fin m) (Fin ℓ) ℂ) (Om : Ω → Matrix (Fin n) (Fin ℓ) ℂ)
    (Phi : Ω → Matrix (Fin m) (Fin s) ℂ) (Psi : Ω → Matrix (Fin n) (Fin s) ℂ)
    (hG : IndependentStandardComplexGaussianTestMatrices μ Ups Om Phi Psi)
    (qrQ : Matrix (Fin m) (Fin ℓ) ℂ → Matrix (Fin m) (Fin (min m ℓ)) ℂ)
    (hqrQ : ∀ Y, IsEconQR Y (qrQ Y))
    (qrP : Matrix (Fin n) (Fin ℓ) ℂ → Matrix (Fin n) (Fin (min n ℓ)) ℂ)
    (hqrP : ∀ X, IsEconQR X (qrP X)) :
    ∫⁻ ω, ENNReal.ofReal (frobSq (A -
        singleViewQ qrQ A (Om ω) *
          singleViewC A (Phi ω) (Psi ω) (singleViewP qrP A (Ups ω)) (singleViewQ qrQ A (Om ω)) *
          (singleViewP qrP A (Ups ω))ᴴ)) ∂μ
      ≤ ENNReal.ofReal ((s : ℝ) / ((s : ℝ) - (ℓ : ℝ)) *
        ⨅ k : Fin ℓ, (((ℓ : ℝ) + ((k : ℕ) : ℝ)) / ((ℓ : ℝ) - ((k : ℕ) : ℝ)) *
          ∑ j ∈ Finset.Icc ((k : ℕ) + 1) (min m n), singularValue A j ^ 2)) := by
  sorry

end RNLAPhaseB.MT2020SingleView