Full singular value decomposition

import RNLAPhaseB.B3Base
import RNLAPhaseB.Complexity

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

namespace RNLAPhaseB.DM2017.FullSVD

/-- The multiset of eigenvalues of the Gram matrix `Aᵀ A` (over `ℝ`, `Aᴴ = Aᵀ`).
The singular values of `A` are the nonnegative square roots of these eigenvalues.
This multiset is computed from `A` itself. -/
noncomputable def gramEigenvalues {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ) : Multiset ℝ :=
  Multiset.map (Matrix.isHermitian_conjTranspose_mul_self A).eigenvalues Finset.univ.val

/-- Squares of the diagonal entries `S_{jj}` of a rectangular `m × n` matrix, indexed by
the `n` columns. A column `j` with no diagonal entry (`j ≥ m`) contributes `0`. -/
def paddedDiagonalSquares {m n : ℕ} (S : Matrix (Fin m) (Fin n) ℝ) : Fin n → ℝ :=
  fun j => if h : j.val < m then (S ⟨j.val, h⟩ j) ^ 2 else 0

/-- Row index in `Fin m` of the `k`-th diagonal position, for `k < min m n`. -/
def diagRow {m n : ℕ} (k : Fin (min m n)) : Fin m := Fin.castLE (Nat.min_le_left m n) k

/-- Column index in `Fin n` of the `k`-th diagonal position, for `k < min m n`. -/
def diagCol {m n : ℕ} (k : Fin (min m n)) : Fin n := Fin.castLE (Nat.min_le_right m n) k

/-- dm2017 `def:full-svd`. Given `A ∈ ℝ^{m×n}`, a full SVD is
`A = U Σ Vᵀ = ∑_{i=1}^{min(m,n)} σ_i u_i v_iᵀ`. Here:
* `U ∈ ℝ^{m×m}` and `V ∈ ℝ^{n×n}` are orthogonal. Their columns `u_i` and `v_i` are the
  left and right singular vectors.
* `Σ` (written `S`) is an `m × n` rectangular diagonal matrix. Its diagonal carries the
  singular values of `A` in decreasing order.

"Singular values of `A`" is realized from `A` itself: the diagonal entries are nonnegative,
and their squares, padded with zeros, form the eigenvalue multiset of `Aᵀ A`. -/
structure IsFullSVD {m n : ℕ} (A : Matrix (Fin m) (Fin n) ℝ)
    (U : Matrix (Fin m) (Fin m) ℝ) (S : Matrix (Fin m) (Fin n) ℝ)
    (V : Matrix (Fin n) (Fin n) ℝ) : Prop where
  /-- `U` is an orthogonal matrix. -/
  U_orthogonal : Uᵀ * U = 1 ∧ U * Uᵀ = 1
  /-- `V` is an orthogonal matrix. -/
  V_orthogonal : Vᵀ * V = 1 ∧ V * Vᵀ = 1
  /-- `Σ` is (rectangular) diagonal. -/
  S_diagonal : ∀ (i : Fin m) (j : Fin n), (i : ℕ) ≠ (j : ℕ) → S i j = 0
  /-- The diagonal entries are nonnegative, as singular values are. -/
  S_diag_nonneg : ∀ k : Fin (min m n), 0 ≤ S (diagRow k) (diagCol k)
  /-- The diagonal entries are in decreasing order. -/
  S_diag_decreasing : ∀ k l : Fin (min m n), k ≤ l →
    S (diagRow l) (diagCol l) ≤ S (diagRow k) (diagCol k)
  /-- The diagonal entries are the singular values of `A`: their squares, padded with
  zeros, are exactly the eigenvalues of `Aᵀ A`, counted with multiplicity. -/
  S_diag_singularValues :
    gramEigenvalues A = Multiset.map (paddedDiagonalSquares S) Finset.univ.val
  /-- `A = U Σ Vᵀ`. -/
  factorization : A = U * S * Vᵀ
  /-- `A = ∑_{i=1}^{min(m,n)} σ_i u_i v_iᵀ`, where `u_i` and `v_i` are the columns of
  `U` and `V`. -/
  rank_one_expansion : A = ∑ k : Fin (min m n),
    S (diagRow k) (diagCol k) •
      Matrix.vecMulVec (fun r => U r (diagRow k)) (fun c => V c (diagCol k))

end RNLAPhaseB.DM2017.FullSVD