Add reduced svd interface and tests (#65)

Author: Cole Scott

src/high-level/lapack-generics.lisp

@@ -10,7 +10,7 @@
 
 (defgeneric lapack-csd (matrix p q))
 
-(defgeneric lapack-svd (matrix))
+(defgeneric lapack-svd (matrix &key reduced))
 
 (defgeneric lapack-ql (matrix)
   (:documentation "Find the LAPACK intermediate representation of ql of a matrix"))

src/high-level/lapack-macros.lisp

@@ -133,43 +133,46 @@
 
 (defmacro def-lapack-svd (class type svd-function &optional real-type)
   `(progn
-     (defmethod svd ((m ,class))
-       (lapack-svd m))
+     (defmethod svd ((m ,class) &key reduced)
+       (lapack-svd m :reduced reduced))
 
-     (defmethod lapack-svd ((m ,class))
-       "Find the SVD of a matrix M. Return (VALUES U SIGMA Vt) where M = U*SIGMA*Vt"
-       (let ((jobu "A")
-             (jobvt "A")
-             (rows (nrows m))
-             (cols (ncols m))
-             (a (alexandria:copy-array (storage (if (eql :row-major (order m)) (transpose m) m))))
-             (lwork -1)
-             (info 0))
-         (let ((lda rows)
-               (s (make-array (min rows cols) :element-type ',(or real-type type)))
-               (ldu rows)
-               (ldvt cols)
-               (work1 (make-array (max 1 lwork) :element-type ',type))
-               (work nil)
-               ,@(when real-type
-                   `((rwork (make-array (* 5 (min rows cols)) :element-type ',real-type)))))
-           (let ((u (make-array (* ldu rows) :element-type ',type))
-                 (vt (make-array (* ldvt cols) :element-type ',type)))
-             ;; run it once as a workspace query
-             (,svd-function jobu jobvt rows cols a lda s u ldu vt ldvt
-                            work1 lwork ,@(when real-type `(rwork)) info)
-             (setf lwork (round (realpart (aref work1 0))))
-             (setf work (make-array (max 1 lwork) :element-type ',type))
-             ;; run it again with optimal workspace size
-             (,svd-function jobu jobvt rows cols a lda s u ldu vt ldvt
-                            work lwork ,@(when real-type `(rwork)) info)
-             (let ((smat (make-array (* rows cols) :element-type ',(or real-type type))))
-               (dotimes (i (min rows cols))
-                 (setf (aref smat (column-major-index (list i i) (shape m)))
-                       (aref s i)))
-               (values (from-array u (list rows rows) :order :column-major)
-                       (from-array smat (list rows cols) :order :column-major)
-                       (from-array vt (list cols cols) :order :column-major)))))))))
+     (defmethod lapack-svd ((m ,class) &key reduced)
+       "Find the SVD of a matrix M. Return (VALUES U SIGMA Vt) where M = U*SIGMA*Vt. If REDUCED is non-NIL, return the reduced SVD (where either U or V are just partial isometries and not necessarily unitary matrices)."
+       (let* ((jobu (if reduced "S" "A"))
+              (jobvt (if reduced "S" "A"))
+              (rows (nrows m))
+              (cols (ncols m))
+              (a (alexandria:copy-array (storage (if (eql :row-major (order m)) (transpose m) m))))
+              (lwork -1)
+              (info 0)
+              (k (min rows cols))
+              (u-cols (if reduced k rows))
+              (vt-rows (if reduced k cols))
+              (lda rows)
+              (s (make-array (min rows cols) :element-type ',(or real-type type)))
+              (ldu rows)
+              (ldvt vt-rows)
+              (work1 (make-array (max 1 lwork) :element-type ',type))
+              (work nil)
+              ,@(when real-type
+                  `((rwork (make-array (* 5 (min rows cols)) :element-type ',real-type)))))
+         (let ((u (make-array (* ldu rows) :element-type ',type))
+               (vt (make-array (* ldvt cols) :element-type ',type)))
+           ;; run it once as a workspace query
+           (,svd-function jobu jobvt rows cols a lda s u ldu vt ldvt
+                          work1 lwork ,@(when real-type `(rwork)) info)
+           (setf lwork (round (realpart (aref work1 0))))
+           (setf work (make-array (max 1 lwork) :element-type ',type))
+           ;; run it again with optimal workspace size
+           (,svd-function jobu jobvt rows cols a lda s u ldu vt ldvt
+                          work lwork ,@(when real-type `(rwork)) info)
+           (let ((smat (make-array (* u-cols vt-rows) :element-type ',(or real-type type))))
+             (dotimes (i k)
+               (setf (aref smat (column-major-index (list i i) (list u-cols vt-rows)))
+                     (aref s i)))
+             (values (from-array u (list rows u-cols) :order :column-major)
+                     (from-array smat (list u-cols vt-rows) :order :column-major)
+                     (from-array vt (list vt-rows cols) :order :column-major))))))))
 
 ;; TODO: This returns only the real parts when with non-complex numbers. Should do something different?
 (defmacro def-lapack-eig (class type eig-function &optional real-type)

src/high-level/matrix.lisp

@@ -523,10 +523,10 @@ If fast is t then just change order. Fast can cause problems when you want to mu
     (declare (ignore matrix))
     (error "INVERSE is not defined for the generic matrix type.")))
 
-(defgeneric svd (matrix)
+(defgeneric svd (matrix &key reduced)
   (:documentation "Find the SVD of a matrix M. Return (VALUES U SIGMA Vt) where M = U*SIGMA*Vt")
-  (:method ((matrix matrix))
-    (declare (ignore matrix))
+  (:method ((matrix matrix) &key reduced)
+    (declare (ignore matrix reduced))
     (error "SVD is not defined for the generic matrix type.")))
 
 (defgeneric qr (matrix)

tests/high-level-tests.lisp

@@ -67,3 +67,62 @@
     (is (cl:= (magicl:nrows xxx) (magicl:ncols xxx) (expt matrix-dim 3)))
     ))
 
+
+(deftest test-svd ()
+  "Test the full and reduced SVDs."
+  (labels ((mul-diag-times-gen (diag matrix)
+             "Returns a newly allocated matrix resulting from the product of DIAG (a diagonal real matrix) with MATRIX (a complex matrix)."
+             #+ignore
+             (declare (type matrix diag matrix)
+                      (values matrix))
+             (let* ((m (magicl:nrows diag))
+                    (k (magicl:ncols matrix))
+                    (result (magicl:empty (list m k))))
+               (dotimes (i (min m (magicl:ncols diag)) result)
+                 (let ((dii (magicl:tref diag i i)))
+                   (dotimes (j k)
+                     (setf (magicl:tref result i j)
+                           (* dii (magicl:tref matrix i j))))))))
+
+           (norm-inf (matrix)
+             "Return the infinity norm of vec(MATRIX)."
+             (let ((data (magicl::storage matrix)))
+               (reduce #'max data :key #'abs)))
+
+           (zero-p (matrix &optional (tolerance 1.0e-14))
+             "Return T if MATRIX is close to zero (within TOLERANCE)."
+             (< (norm-inf matrix) tolerance))
+
+           (check-full-svd (matrix)
+             "Validate full SVD of MATRIX."
+             (let ((m (magicl:nrows matrix))
+                   (n (magicl:ncols matrix)))
+               (multiple-value-bind (u sigma vh)
+                   (magicl:svd matrix)
+                 (is (= (magicl:nrows u) (magicl:ncols u) m))
+                 (is (and (= (magicl:nrows sigma) m) (= (magicl:ncols sigma) n)))
+                 (is (= (magicl:nrows vh) (magicl:ncols vh) n))
+                 (is (zero-p (magicl:- matrix (magicl:@ u (mul-diag-times-gen sigma vh))))))))
+
+           (check-reduced-svd (matrix)
+             "Validate reduced SVD of MATRIX."
+             (let* ((m (magicl:nrows matrix))
+                    (n (magicl:ncols matrix))
+                    (k (min m n)))
+
+               (multiple-value-bind (u sigma vh)
+                   (magicl:svd matrix :reduced t)
+                 (is (and (= (magicl:nrows u) m)
+                          (= (magicl:ncols u) k)))
+                 (is (= (magicl:nrows sigma) (magicl:ncols sigma) k))
+                 (is (and (= (magicl:nrows vh) k)
+                          (= (magicl:ncols vh) n)))
+                 (is (zero-p (magicl:- matrix (magicl:@ u (mul-diag-times-gen sigma vh)))))))))
+
+    (let ((tall-thin-matrix (magicl:rand '(8 2))))
+      (check-full-svd tall-thin-matrix)
+      (check-reduced-svd tall-thin-matrix))
+
+    (let ((short-fat-matrix (magicl:rand '(2 8))))
+      (check-full-svd short-fat-matrix)
+      (check-reduced-svd short-fat-matrix))))