search

sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.48k stars 488 forks source link

Matrix: Add generic SVD and singular_values methods using sympy #32171

Open mkoeppe opened 3 years ago

mkoeppe commented 3 years ago

(from #31992)

31942 added _sympy_ methods that provide conversion of Sage matrices and vectors to Sympy, as well as the conversion from Sympy to Sage using monkey-patched _sage_ methods.

Based on this, we can now implement various matrix methods using Sympy.

In this ticket, we add SVD (singular value decomposition) https://docs.sympy.org/latest/modules/matrices/matrices.html

This provides exact and symbolic SVD. Sage currently only has one implementation, for RDF, using scipy.linalg.

CC: @honglizhaobob

Component: linear algebra

Branch/Commit: u/mkoeppe/matrix__add_generic_svd_and_singular_values_methods_using_sympy @ 130033d

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

mkoeppe commented 3 years ago

Description changed:

--- 
+++ 
@@ -1,6 +1,7 @@
 (from #31992)

-#31942 added `_sympy_` methods for matrices and vectors.
+#31942 added `_sympy_` methods that provide conversion of Sage matrices and vectors to Sympy, as well as the conversion from Sympy to Sage using monkey-patched `_sage_` methods.
+

 Based on this, we can now implement various matrix methods using Sympy.
mkoeppe commented 3 years ago
comment:2

This would go into src/sage/matrix/matrix2.pyx

mkoeppe commented 3 years ago

Branch: u/mkoeppe/matrix__add_generic_svd_and_singular_values_methods_using_sympy

mwageringel commented 3 years ago

Commit: 130033d

mwageringel commented 3 years ago
comment:5

I think the documentation should make very clear that in 99% of the cases, it is better to use scipy and work over RDF/CDF if one expects a numerical result, rather than using sympy with RR. The sympy implementation gets very slow for larger matrices.

Also, even though sympy seems to support singular value decompositions up to precision 83 (or at least does not raise an error), the result is not really much more accurate compared to RDF. In fact, the result has 53-bit precision as well.

sage: a = matrix.random(RealField(83), 10)
sage: b = a._sympy_().as_mutable()
sage: U,S,V = b.singular_value_decomposition()
sage: u,s,v = a.change_ring(RDF).SVD()
sage: (U * S * V.H - b).norm('fro'), (u * s * v.H - a).norm('frob')
(1.73609121130844744221236e-15, 5.9963160390748384e-15)
sage: S[0,0]._prec
53

New commits:

130033dMatrix.SVD: Create stub; Matrix_double_dense.SVD: Add algorithm keyword
mkoeppe commented 3 years ago
comment:6

That's a great point, thanks for looking into this. The main use case for using sympy here should be symbolic SVD.