Error bound for SVD via row extraction

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

/-- The fixed-precision condition [eq:fixed_precision2] (restating [eq:fixed_precision]).
`ε` is a positive error tolerance, `Q` has `k` orthonormal columns, and
`‖A - Q Qᵀ A‖ ≤ ε` holds in the `ℓ₂` operator norm. -/
def FixedPrecisionCondition {m n k : ℕ} (A : Matrix (Fin m) (Fin n) ℝ)
    (Q : Matrix (Fin m) (Fin k) ℝ) (ε : ℝ) : Prop :=
  0 < ε ∧ OrthonormalColumns Q ∧ ‖A - Q * Qᵀ * A‖ ≤ ε

/-- The interpolative decomposition `Q = X Q_(J,:)` as described in the source:
`J` marks `k` distinct rows of `Q`, `X` is `m × k` with entries bounded in
magnitude by two, and `X_(J,:) = I_k`. -/
def IsInterpolativeDecomposition {m k : ℕ} (Q : Matrix (Fin m) (Fin k) ℝ)
    (J : Fin k → Fin m) (X : Matrix (Fin m) (Fin k) ℝ) : Prop :=
  Function.Injective J ∧
  Q = X * Q.submatrix J id ∧
  X.submatrix J id = 1 ∧
  ∀ i j, |X i j| ≤ 2

/-- Output relation of Algorithm [alg:extractrows] (SVD via Row Extraction).
1. ID `Q = X Q_(J,:)`;
2. QR factorization `A_(J,:) = Rᵀ Wᵀ` (`W` orthonormal columns, `R` upper triangular);
3. `Z = X Rᵀ`;
4. SVD `Z = U Σ Ṽᵀ` (`U` orthonormal columns, `Σ` nonnegative diagonal, `Ṽ` orthogonal);
5. `V = W Ṽ`. -/
def RowExtractionSVDOutput {m n k : ℕ} (A : Matrix (Fin m) (Fin n) ℝ)
    (Q : Matrix (Fin m) (Fin k) ℝ) (U : Matrix (Fin m) (Fin k) ℝ)
    (S : Matrix (Fin k) (Fin k) ℝ) (V : Matrix (Fin n) (Fin k) ℝ) : Prop :=
  ∃ (J : Fin k → Fin m) (X : Matrix (Fin m) (Fin k) ℝ)
    (R : Matrix (Fin k) (Fin k) ℝ) (W : Matrix (Fin n) (Fin k) ℝ)
    (Z : Matrix (Fin m) (Fin k) ℝ) (Vtilde : Matrix (Fin k) (Fin k) ℝ),
    IsInterpolativeDecomposition Q J X ∧
    OrthonormalColumns W ∧ R.BlockTriangular id ∧
    A.submatrix J id = Rᵀ * Wᵀ ∧
    Z = X * Rᵀ ∧
    OrthonormalColumns U ∧
    (∃ s : Fin k → ℝ, (∀ i, 0 ≤ s i) ∧ S = Matrix.diagonal s) ∧
    Vtildeᵀ * Vtilde = 1 ∧
    Z = U * S * Vtildeᵀ ∧
    V = W * Vtilde

/-- HMT2011 Lemma [lem:convert_ID]: if `A` (`m × n`) and `Q` (`m × k`) satisfy
[eq:fixed_precision2] and `U, Σ, V` are constructed by Algorithm [alg:extractrows], then
`‖A - U Σ Vᵀ‖ ≤ [1 + sqrt(1 + 4k(n-k))] ε` in the spectral norm.

SOURCE WARNING sw1: the source defines `X` as an `m × k` matrix containing `I_k` as a
submatrix, so its remaining block is `(m-k) × k`; yet the proof speaks of the
"remaining `(n-k) × k` submatrix" and the stated bound [eq:IPA] uses `n-k`. These
dimensions do not match (the natural bound on `‖X‖` from the proof argument would involve
`m-k`). As instructed, the bound is stated exactly as written in the source, with `n - k`
(computed in ℝ), with no added hypothesis relating `m` and `n`. -/
theorem convert_ID {m n k : ℕ} (A : Matrix (Fin m) (Fin n) ℝ)
    (Q : Matrix (Fin m) (Fin k) ℝ) (ε : ℝ)
    (U : Matrix (Fin m) (Fin k) ℝ) (S : Matrix (Fin k) (Fin k) ℝ)
    (V : Matrix (Fin n) (Fin k) ℝ)
    (hfp : FixedPrecisionCondition A Q ε)
    (halg : RowExtractionSVDOutput A Q U S V) :
    ‖A - U * S * Vᵀ‖ ≤ (1 + Real.sqrt (1 + 4 * (k : ℝ) * ((n : ℝ) - (k : ℝ)))) * ε := by
  sorry

end RNLAPhaseB.HMT2011.ConvertID