Optimality of truncated SVD

import RNLAPhaseB.B3Base
import RNLAPhaseB.Complexity

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

namespace RNLAPhaseB.DM2017BestRankK

/-- Thin SVD `A = U Σ Vᵀ` of `A ∈ ℝ^{m×n}` with `ρ = rank(A)`.
`U ∈ ℝ^{m×ρ}` and `V ∈ ℝ^{n×ρ}` have orthonormal columns. `Σ = diag(σ₁,…,σ_ρ)` is
`ρ × ρ` with `σ₁ ≥ … ≥ σ_ρ > 0`, and index `i : Fin ρ` stands for `σ_{i+1}`.
The source's Linear Algebra section (book.md:105) takes the basic operations of transposition
and matrix multiplication as known. They are realized here by Mathlib's `ᵀ` and `*`. -/
def IsThinSVD {m n ρ : ℕ} (A : Matrix (Fin m) (Fin n) ℝ)
    (U : Matrix (Fin m) (Fin ρ) ℝ) (σ : Fin ρ → ℝ) (V : Matrix (Fin n) (Fin ρ) ℝ) : Prop :=
  OrthonormalColumns U ∧ OrthonormalColumns V ∧ Antitone σ ∧ (∀ i, 0 < σ i) ∧
    A = U * Matrix.diagonal σ * Vᵀ

/-- `A_k = ∑_{i=1}^k σ_i u_i v_iᵀ = U_k Σ_k V_kᵀ`, where `u_i`, `v_i` are the columns of
`U`, `V`. In 0-based indexing these are the indices `i` with `i.val < k`.
The outer product `u_i v_iᵀ` is one of the basic operations assumed by the source's Linear Algebra
section (book.md:105). It is realized by `Matrix.vecMulVec`, whose `(a,b)` entry is
`u_i a * v_i b`. -/
noncomputable def truncatedSVD {m n ρ : ℕ} (U : Matrix (Fin m) (Fin ρ) ℝ) (σ : Fin ρ → ℝ)
    (V : Matrix (Fin n) (Fin ρ) ℝ) (k : ℕ) : Matrix (Fin m) (Fin n) ℝ :=
  ∑ i ∈ Finset.univ.filter (fun i : Fin ρ => i.val < k),
    σ i • Matrix.vecMulVec (fun a => U a i) (fun b => V b i)

/-- The squared Frobenius norm `‖M‖_F^2 = ∑_{i,j} M_{ij}^2`. -/
def frobeniusNormSq {m n : ℕ} (M : Matrix (Fin m) (Fin n) ℝ) : ℝ :=
  ∑ i, ∑ j, M i j ^ 2

/-- Best rank-k approximation (Eckart–Young–Mirsky), as stated in dm2017.
Let `A = UΣVᵀ` be the thin SVD of `A`, with `k < rank(A) = ρ` an integer.
The first part is `σ_{k+1} = min_{rank(B)=k} ‖A - B‖₂ = ‖A - A_k‖₂`.
The second part is `∑_{j=k+1}^ρ σ_j² = min_{rank(B)=k} ‖A - B‖_F² = ‖A - A_k‖_F²`.
Here `‖·‖` is the spectral (ℓ₂ operator) norm. `σ_{k+1}` is `σ ⟨k, hk⟩`, and the tail
sum runs over the 0-based indices `j` with `k ≤ j.val`. -/
theorem best_rank_k_approximation {m n ρ : ℕ} (A : Matrix (Fin m) (Fin n) ℝ)
    (U : Matrix (Fin m) (Fin ρ) ℝ) (σ : Fin ρ → ℝ) (V : Matrix (Fin n) (Fin ρ) ℝ)
    (hSVD : IsThinSVD A U σ V) (hρ : A.rank = ρ) (k : ℕ) (hk : k < ρ) :
    (IsLeast {x : ℝ | ∃ B : Matrix (Fin m) (Fin n) ℝ, B.rank = k ∧ x = ‖A - B‖}
        (σ ⟨k, hk⟩) ∧
      ‖A - truncatedSVD U σ V k‖ = σ ⟨k, hk⟩) ∧
    (IsLeast {x : ℝ | ∃ B : Matrix (Fin m) (Fin n) ℝ, B.rank = k ∧ x = frobeniusNormSq (A - B)}
        (∑ j ∈ Finset.univ.filter (fun j : Fin ρ => k ≤ j.val), σ j ^ 2) ∧
      frobeniusNormSq (A - truncatedSVD U σ V k) =
        ∑ j ∈ Finset.univ.filter (fun j : Fin ρ => k ≤ j.val), σ j ^ 2) := by
  sorry

end RNLAPhaseB.DM2017BestRankK