FiniteRankFreeModuleMorphism: Add method SVD (singular value decomposition)

[mkoeppe]

This operation is currently available for matrices over RDF.

The version for FiniteRankFreeModuleMorphism would define a new basis on each of the domain and codomain (with an orthogonal change of basis) so that the matrix of the morphism in the new bases is the diagonal of the singular values.

Example:

sage: M = FiniteRankFreeModule(ZZ, 3, name='M') 
sage: N = FiniteRankFreeModule(ZZ, 2, name='N') 
sage: e = M.basis('e') ; f = N.basis('f') 
sage: H = Hom(M,N) ; H 
Set of Morphisms from Rank-3 free module M over the Integer Ring to Rank-2 free module N over the Integer Ring in Category of finite dimensional modules over Integer Ring
sage: phi = H([[2,-1,3], [1,0,-4]], name='phi', latex_name=r'\phi') ; phi 
Generic morphism:
  From: Rank-3 free module M over the Integer Ring
  To:   Rank-2 free module N over the Integer Ring
sage: phi.matrix(e, f)                                                   
[ 2 -1  3]
[ 1  0 -4]
sage: f_prime = N.basis('f_prime', from_family=(f[0]-f[1], -2*f[0]+3*f[1]
# pretend that this basis is the right orthogonal matrix of the SVD...
sage: phi.matrix(e, f_prime)                                             
[ 8 -3  1]
[ 3 -1 -1]

CC: @egourgoulhon @mjungmath @honglizhaobob

Component: linear algebra

Branch/Commit: u/gh-honglizhaobob/finiterankfreemodulemorphism__add_method_svd__singular_value_decomposition_ @ 38bcecc

Issue created by migration from https://trac.sagemath.org/ticket/31992

[mkoeppe]

Description changed:

--- 
+++ 
@@ -2,3 +2,24 @@
 
 The version for `FiniteRankFreeModuleMorphism` would define a new basis on each of the domain and codomain (with an orthogonal change of basis) so that the matrix of the morphism in the new bases is the diagonal of the singular values.
 
+Example:
+
+```
+sage: M = FiniteRankFreeModule(ZZ, 3, name='M') 
+sage: N = FiniteRankFreeModule(ZZ, 2, name='N') 
+sage: e = M.basis('e') ; f = N.basis('f') 
+sage: H = Hom(M,N) ; H 
+Set of Morphisms from Rank-3 free module M over the Integer Ring to Rank-2 free module N over the Integer Ring in Category of finite dimensional modules over Integer Ring
+sage: phi = H([[2,-1,3], [1,0,-4]], name='phi', latex_name=r'\phi') ; phi 
+Generic morphism:
+  From: Rank-3 free module M over the Integer Ring
+  To:   Rank-2 free module N over the Integer Ring
+sage: phi.matrix(e, f)                                                   
+[ 2 -1  3]
+[ 1  0 -4]
+sage: f_prime = N.basis('f_prime', from_family=(f[0]-f[1], -2*f[0]+3*f[1]
+# pretend that this basis is the right orthogonal matrix of the SVD...
+sage: phi.matrix(e, f_prime)                                             
+[ 8 -3  1]
+[ 3 -1 -1]
+```

[honglizhaobob]

comment:2

Replying to @mkoeppe:

Sorry for asking a naive question... In the description, what does RDF refer to ("this operation is currently available for matrices over RDF")?

[mkoeppe]

comment:3

RDF is the "real double-precision field" - equivalent of double in C.
https://doc.sagemath.org/html/en/reference/rings_numerical/sage/rings/real_double.html

[honglizhaobob]

comment:4

Replying to @mkoeppe:

RDF is the "real double-precision field" - equivalent of double in C.
https://doc.sagemath.org/html/en/reference/rings_numerical/sage/rings/real_double.html

Thanks professor. Additionally, in implementing the SVD, should we resort to external libraries such as numpy, or is it a better practice to use sage methods whenever possible?

[honglizhaobob]

comment:5

If I'm understanding correctly after reading the code, we are still working with matrices (instead of high-dimensional tensors). Suppose A is the matrix obtained from M.hom(N, some_matrix, basis=(e, f)). Let x.parent() == M, y.parent() == N. Then we have:

y = A*x
y = U * S * V_dagger * x            # do SVD on A
(U_dagger * y) = S * (V_dagger * x) # invert U
y_prime = S*x_prime                 # renamed x and y

In code, we will reassign the default bases of M and N:

M._def_basis = V
N._def_basis = U                    # certainly, this reassignment needs more
                                    # work in terms of actual implementations

such that S would be the result of hom.matrix() in this new basis.

[mkoeppe]

comment:6

Replying to @honglizhaobob:

in implementing the SVD, should we resort to external libraries such as numpy, or is it a better practice to use sage methods whenever possible?

In Sage we support a wide range of rings for the arithmetic, and consequently often use multiple libraries for the implementation and also have an in-Sage generic implementation.

The existing implementation of SVD in src/sage/matrix/matrix_double_dense.pyx just calls scipy.linalg.svd.

It would be quite useful to add a generic implementation that only uses the general methods that Sage matrices offer. In that way, we could have access to arbitrary precision and symbolic computations.

It would also be good to check what other libraries that we already use or could easily use also provide implementations of SVD.

[honglizhaobob]

comment:7

Replying to @mkoeppe:

Replying to @honglizhaobob:

in implementing the SVD, should we resort to external libraries such as numpy, or is it a better practice to use sage methods whenever possible?

In Sage we support a wide range of rings for the arithmetic, and consequently often use multiple libraries for the implementation and also have an in-Sage generic implementation.

The existing implementation of SVD in src/sage/matrix/matrix_double_dense.pyx just calls scipy.linalg.svd.

It would be quite useful to add a generic implementation that only uses the general methods that Sage matrices offer. In that way, we could have access to arbitrary precision and symbolic computations.

It would also be good to check what other libraries that we already use or could easily use also provide implementations of SVD.

I think currently SVD() doesn't support other fields such as ZZ, or QQ. We may be restricted to RDF if we choose to use matrix_double_dense.pyx.

[mkoeppe]

comment:8

Yes, that's right -- and that's I think also exactly what scipy.linalg.svd provides. That's why I said that it will be good to check what other implementations are available in other software. A quick search for "arbitrary precision svd" did not lead me to good results though

[honglizhaobob]

comment:9

Replying to @mkoeppe:

Yes, that's right -- and that's I think also exactly what scipy.linalg.svd provides. That's why I said that it will be good to check what other implementations are available in other software. A quick search for "arbitrary precision svd" did not lead me to good results though

Could we consider symbolic SVD such as in MATLAB or sympy？ This way we can preserve the QQ structure, however I feel svd on ZZ will be a bit demanding..

[mkoeppe]

comment:10

Replying to @honglizhaobob:

Could we consider symbolic SVD such as in MATLAB or sympy？ This way we can preserve the QQ structure

Do you have an example or link that illustrates what you have in mind?

[honglizhaobob]

Branch: u/gh-honglizhaobob/finiterankfreemodulemorphism__add_method_svd__singular_value_decomposition_