Singular Value Decomposition SVD

Stream: general

[MohanadAhmed]

Hello everyone,

So I just finished an attempt at proving that there exists a Singular Value Decomposition (SVD) for complex matrices. The code is on Github. It is still very inefficient (takes about 1 minute for lean to go through). I am hoping to take another run at it in a more structured manner, and then hopefully port it to lean4/mathlib4 (once I learn enough lean4/mathlib4).

I thought I would share it with the community first, see any thoughts on how to structure it better and make it more useful. In its current form the main statement is basically:

lemma svd_decompose{m n r: ℕ} (A: matrix (fin m) (fin n) ℂ) :
r = A.rank → ∃
  (U: matrix (fin m) (fin r ⊕ fin (m - r)) ℂ)
  (Q: matrix (fin r ⊕ fin (m - r)) (fin r ⊕ fin (n - r)) ℝ)
  (V: matrix (fin n) (fin r ⊕ fin (n - r)) ℂ),
    A = U⬝(Q.map RηC)⬝Vᴴ ∧ /- [1] SVD Statement -/
    V⬝Vᴴ = 1 ∧ /- [2] SVD Right Eigenvectors are orthogonoal-/
    U⬝Uᴴ = 1 ∧/- [3] SVD Left Eigenvectors are orthogonoal-/
    Uᴴ⬝U = 1 ∧ /- [4] Orthogonality [3] the other way around -/
    Vᴴ⬝V = 1 ∧ /- [5] Orthogonality [2] the other way around -/
Q.to_blocks₁₂ = 0 ∧ Q.to_blocks₂₁ = 0 ∧  Q.to_blocks₂₂ = 0 /- [6] Off diagonal blocks are zero -/
     ∧ (∀ i j, i ≠ j → Q.to_blocks₁₁ i j = 0) ∧ (∀ i, 0 < Q.to_blocks₁₁ i i) /- [7] Singular values (Non-zero) are positive and sitting in a diagonal matrix -/
      := /- rest of proof here -/

1. My main thought was, SVD will be used as a block in another proof. The above lemma says: "give me the SVD decomposition with its properties". Any suggestions for a better API for SVD?

2. The proof is also currently for complex matrices but will probably not be very difficult to modify to include the reals. The eigendecomposition (on which SVD rests )is based on the file linear_algebra.matrix.spectrum which works for is_R_or_C. Are there any other cases that are worth paying attention to?

Any other thoughts or suggestions are also welcome!

[Patrick Stevens]

A programmer might answer your point 1 as: write a def which builds the SVD, and then separately prove that your def satisfies the properties of the SVD. That way you can obtain the SVD without having to extract some specific part of an existential statement. For the first example of this pattern which GitHub search gave me, see e.g. the definition of Cantor normal form of ordinals (https://github.com/leanprover-community/mathlib4/blob/6a12e9c90cfa8b8026aeb23b4aeb80737ae1f940/Mathlib/SetTheory/Ordinal/CantorNormalForm.lean#L78), which defines the CNF without any properties, and then goes on to prove in many small chunks that this actually is a normal form.

[MohanadAhmed]

So two questions:

1. What is preferable about this pattern compared to other methods? It seems it is quite common in mathlib, for example in docs#linear_algebra.matrix.spectrum the eigenvalues of a hermitian are treated in this way.

2. How do you work with the inference system to make sure lean understands you are always referring to the same matrix all the time. I was trying to mockup the suggestion you made above. I defined the three matrices U, Q, V and then wanted to state the 7 properties of SVD above sorried out. It seems I cannot convince lean to infer the matrix A in the 2nd lemma below. Any thoughts on how to write this one out in the shortest / idiomatic way?

import data.matrix.basic
import data.matrix.notation
import data.complex.basic

open matrix
open_locale matrix big_operators

namespace svd

variables {m n: ℕ}
variables (A: matrix (fin m) (fin n) ℂ)


def RηC := algebra_map ℝ ℂ

def V : matrix (fin n) (fin n) ℂ := sorry
def U : matrix (fin n) (fin n) ℂ := sorry
def Q : matrix (fin m) (fin n) ℝ := sorry

lemma svd_eq: A = svd.U⬝(svd.Q.map RηC)⬝svd.V := sorry /- This one works out -/
lemma right_conj_transpose_eigenvector_mul: (svd.Vᴴ) ⬝ svd.V   = 1 := sorry /- Does not typecheck -/

-- More lemmas about the U, Q, V matrices should be here

end svd

[Mario Carneiro]

I think you need to fake a dependency on A in V U Q, or else use include A, because as written lean will assume V doesn't depend on A

[Mario Carneiro]

this issue is mostly an artifact of filling the definitions by sorry instead of the real thing (which presumably refers to A)

[MohanadAhmed]

Mario Carneiro said:

I think you need to fake a dependency on A in V U Q, or else use include A, because as written lean will assume V doesn't depend on A

So I tried what I understood from your first suggestion (about fake dependency). The A is not recognized inside the V defintion. Any thoughts on what I am missing?

import data.matrix.basic
import data.matrix.notation
import data.complex.basic
import linear_algebra.matrix.pos_def
import linear_algebra.matrix.spectrum

open matrix
open_locale matrix big_operators

namespace svd

variables {m n: ℕ}
variables (A: matrix (fin m) (fin n) ℂ)


def RηC := algebra_map ℝ ℂ
-- include A
noncomputable def V : matrix (fin n) (fin n) ℂ := begin
  let eigs := (is_hermitian_transpose_mul_self A).eigenvalues, /- Unknown identifier A-/
  exact (diagonal eigs),
end
def U : matrix (fin n) (fin n) ℂ := sorry
def Q : matrix (fin m) (fin n) ℝ := sorry

lemma svd_eq: A = svd.U⬝(svd.Q.map RηC)⬝svd.V := sorry
lemma right_conj_transpose_eigenvector_mul: (svd.Vᴴ) ⬝ svd.V = 1 := sorry /- Still does not typecheck-/

-- More lemmas about the U, Q, V matrices should be here

end svd

[Mario Carneiro]

if you do it inside begin ... end it will have already decided that A is not part of the definition

[Mario Carneiro]

so you have to use include or put the let outside begin ... end

[MohanadAhmed]

So you mean it must be part of the signature?

[Mario Carneiro]

yes, that will also work

[Mario Carneiro]

in general variable declarations are only included in definitions when they are used

[MohanadAhmed]

So if there is not begin end but I just a direct term, will that work?

[Mario Carneiro]

yes

[MohanadAhmed]

I guess not! Just tried it and now both svd_eq and right_conj_transpose_eigenvector_mul do not typecheck!!

import data.matrix.basic
import data.matrix.notation
import data.complex.basic
import linear_algebra.matrix.pos_def
import linear_algebra.matrix.spectrum

open matrix
open_locale matrix big_operators

namespace svd

variables {m n: ℕ}
variables (A: matrix (fin m) (fin n) ℂ)


def RηC := algebra_map ℝ ℂ
-- include A
noncomputable def V : matrix (fin n) (fin n) ℂ :=
((matrix.diagonal (is_hermitian_transpose_mul_self A).eigenvalues).map RηC)

noncomputable def U : matrix (fin m) (fin m) ℂ :=
((matrix.diagonal (is_hermitian_mul_conj_transpose_self A).eigenvalues).map RηC)

def Q : matrix (fin m) (fin n) ℝ := (A.map (λ x: ℂ, x.re))

/-Now this one does not typecheck as well-/
lemma svd_eq: A = svd.U⬝(svd.Q.map RηC)⬝svd.V := sorry
/- Still does not typecheck-/
lemma right_conj_transpose_eigenvector_mul: (svd.Vᴴ) ⬝ svd.V = 1 := sorry

-- More lemmas about the U, Q, V matrices should be here

end svd

[Mario Carneiro]

right, that's the other part, svd.V and svd.Q are functions, you need to apply them to A

[Mario Carneiro]

so svd_eq would be something like A = svd.U A ⬝ (svd.Q A).map RηC ⬝ svd.V A

[MohanadAhmed]

OK I see! so I will always be mentioning A.
Is there a way to make svd act like the is_hermitian docs#matrix.is_hermitian.eigenvalues where the hypothesis (hA: is_hermitian A) allows you to write hA.eigenvalues andhA.eigenvector_matrix and so on?

[Mario Carneiro]

if you put it in the namespace of A's type you can write it as A.V and A.Q etc, although the name matrix.V is kind of vague

[Eric Wieser]

Another approach (which more closely matches things like numpy etc) would be to havematrix.svd (A : matrix (fin n) (fin n) ℂ) : matrix.svd_parts where structure svd_parts := (V : _) (U : _) (Q : _)

[cvlvxi]

hello