Changeset 8023:3a19fbc01ae5

Timestamp:

12/11/07 04:17:47 (5 years ago)

Author:

Mike Hansen <mhansen@…>

Branch:

default

Message:

Added SVD to matrix_complex_double_dense.pyx and numpy SVD to matrix_real_double_dense.pyx

Location:

sage/matrix

Files:

• matrix_complex_double_dense.pyx (modified) (2 diffs)
• matrix_real_double_dense.pyx (modified) (7 diffs)

sage/matrix/matrix_complex_double_dense.pyx

@@ -175 +175 @@
                 for i from 0<=i<self._ncols:
                     gsl_matrix_complex_set(self._matrix,i,i,z._complex)
+                    
     cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value):
         cdef ComplexDoubleElement z
@@ -826 +827 @@
         return ans
 
-
+    def SVD(self):
+        r"""
+        Return the singular value decomposition of this matrix.
+
+        INPUT:
+            A -- a matrix
+        OUTPUT:
+            U, S, V -- matrices such that A = U * S * V^t, where
+                       U and V are orthogonal and S is diagonal.
+
+
+        EXAMPLES:
+            sage: m = matrix(CDF,4,range(16))
+            sage: U,S,V = m.SVD()
+            sage: U*S*V.transpose()    # slightly random output (due to computer architecture)
+            [3.45569519412e-16               1.0               2.0               3.0]
+            [4.0               5.0               6.0               7.0]
+            [8.0               9.0              10.0              11.0]
+            [12.0              13.0              14.0              15.0]
+
+        A non-square example:
+            sage: m = matrix(CDF, 2, range(6)); m
+            [  0 1.0 2.0]
+            [3.0 4.0 5.0]
+            sage: U, S, V = m.SVD()
+            sage: U
+            [-0.274721127897 -0.961523947641]
+            [-0.961523947641  0.274721127897]
+            sage: S
+            [7.34846922835             0             0]
+            [            0           1.0             0]
+            sage: V
+            [-0.392540507864  0.824163383692  0.408248290464]
+            [-0.560772154092  0.137360563949 -0.816496580928]
+            [ -0.72900380032 -0.549442255795  0.408248290464]
+            sage: U*S*V.transpose()           # random low bits
+            [7.62194127257e-17               1.0               2.0]
+            [              3.0               4.0               5.0]
+            sage: m = matrix(CDF,3,2,range(6)); m
+            [  0 1.0]
+            [2.0 3.0]
+            [4.0 5.0]
+            sage: U,S,V = m.SVD()
+            sage: U*S*V.transpose()   # random low order bits
+            [-8.13151629364e-19                1.0]
+            [               2.0                3.0]
+            [               4.0                5.0]            
+
+        TESTS:
+            sage: m = matrix(CDF, 3, 0, []); m
+            []
+            sage: m.SVD()
+            ([], [], [])
+            sage: m = matrix(CDF, 0, 3, []); m
+            []        
+            sage: m.SVD()
+            ([], [], [])
+        """
+        if self._nrows == 0 or self._ncols == 0:
+            U_t = self.new_matrix(self._nrows, self._ncols)
+            S_t = self.new_matrix(self._nrows, self._ncols)
+            V_t = self.new_matrix(self._ncols, self._nrows)
+            return U_t, S_t, V_t
+
+        import numpy.linalg
+        cdef int i, s_dim
+        P = self.parent()
+        CDF = P.base_ring()
+        
+        U,_S,V = numpy.linalg.svd(self.numpy())
+
+        #Create the inner diagonal matrix
+        s_dim = len(_S)
+        S = matrix(CDF, self._nrows, self._ncols, 0)
+        for i from 0 <= i < s_dim:
+            S[(i,i)] = _S[i]
+
+        return (matrix(U),S,matrix(V).transpose())
+                                 
+        

sage/matrix/matrix_real_double_dense.pyx

@@ -629 +629 @@
         return trans
 
-    def SVD(self):
+    def SVD(self, algorithm='gsl'):
         r"""
         Return the singular value decomposition of this matrix.
@@ -635 +635 @@
         INPUT:
             A -- a matrix
+            algorithm -- 'numpy' or 'gsl'
         OUTPUT:
             U, S, V -- matrices such that A = U * S * V^t, where
@@ -642 +643 @@
         EXAMPLES:
             sage: m = matrix(RDF,4,range(16))
-            sage: U,S,V = m.SVD()
+            sage: U,S,V = m.SVD(algorithm='gsl')
             sage: U*S*V.transpose()    # slightly random output (due to computer architecture)
             [3.45569519412e-16               1.0               2.0               3.0]
@@ -653 +654 @@
             [0.0 1.0 2.0]
             [3.0 4.0 5.0]
+            sage: U, S, V = m.SVD(algorithm='numpy')
+            sage: U*S*V.transpose()           # random low bits
+            [7.62194127257e-17               1.0               2.0]
+            [              3.0               4.0               5.0]            
             sage: U, S, V = m.SVD()
             sage: U
@@ -671 +676 @@
             [2.0 3.0]
             [4.0 5.0]
+            sage: U,S,V = m.SVD(algorithm='numpy')
+            sage: U*S*V.transpose()   # random low order bits
+            [-8.13151629364e-19                1.0]
+            [               2.0                3.0]
+            [               4.0                5.0]  
             sage: U,S,V = m.SVD()
             sage: U*S*V.transpose()   # random low order bits
             [-8.13151629364e-19                1.0]
             [               2.0                3.0]
-            [               4.0                5.0]            
+            [               4.0                5.0]
+
 
         TESTS:
@@ -687 +698 @@
             ([], [], [])
         """
+        if algorithm == 'numpy':
+            return self._SVD_numpy()
+        elif algorithm == 'gsl':
+            return self._SVD_gsl()
+        else:
+            raise ValueError, "unknown algorithm"
+        
+    def _SVD_gsl(self):
+        """
+        Return the singular value decomposition of this matrix
+        using GSL.  Note that the matrices the dimensions of the
+        matrices that this returns are (m,p), (p,p), and (n, p)
+        where p = min(m,n).
+
+        EXAMPLES:
+            sage: def shape(x): return (x.nrows(), x.ncols())
+            sage: m = matrix(RDF, 2, 3, range(6))
+            sage: map(shape, m._SVD_gsl())
+            [(2, 2), (2, 2), (3, 2)]
+
+        """
         if self._ncols > self._nrows:
             m = self.transpose()
             V_t,S_t,U_t = m.SVD()
             return U_t,S_t,V_t
+        
         if self._nrows == 0 or self._ncols == 0:
             U_t = self.new_matrix(self._nrows, self._ncols)
@@ -719 +752 @@
         return A,_S,V
 
+    def _SVD_numpy(self):
+        """
+        Return the singular value decomposition of this matrix
+        using GSL.  Note that the matrices the dimensions of the
+        matrices that this returns are (m,m), (m,n), and (n, n).
+
+        EXAMPLES:
+            sage: def shape(x): return (x.nrows(), x.ncols())
+            sage: m = matrix(RDF, 2, 3, range(6))
+            sage: map(shape, m._SVD_numpy())
+            [(2, 2), (2, 3), (3, 3)]
+
+        """
+        if self._nrows == 0 or self._ncols == 0:
+            U_t = self.new_matrix(self._nrows, self._ncols)
+            S_t = self.new_matrix(self._nrows, self._ncols)
+            V_t = self.new_matrix(self._ncols, self._nrows)
+            return U_t, S_t, V_t
+
+        import numpy.linalg
+        cdef int i, s_dim
+        P = self.parent()
+        CDF = P.base_ring()
+        
+        U,_S,V = numpy.linalg.svd(self.numpy())
+
+        #Create the inner diagonal matrix
+        s_dim = len(_S)
+        S = matrix(CDF, self._nrows, self._ncols, 0)
+        for i from 0 <= i < s_dim:
+            S[(i,i)] = _S[i]
+
+        return (matrix(U),S,matrix(V).transpose())
+                                 
+        
+
+
+
     def QR(self):
         """