(defun size (&rest x)
"
;; Description: function size (&rest x)
;; Args: object(s) such as vectors, lists or arrays. Calculates the array
;; dimensions. Vectorized.
;; Examples:
(size #2a((1 2 3)(4 5 6)))
(size #(1 2 3))
(size (iseq 4))
"
(if (= (length x) 1)
(let ((obj (first x)))
(if (or (listp obj) (vectorp obj))
(length obj)
(array-dimensions obj)
)
)
(mapcar #'size x)
)
)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Description for SVD:
;;
(defun svd (a)
"
;; Description: function svd
performs the Singular Value Decomposition of a matrix. The native svd code
(sv-decomp) requires that the number of rows be >= the number of columns. This
being unnecessary, svd handles either case.
A = U L V'
where U and V are orthogonal matrices, and L is a 'diagonal' matrix containing
the singular values on the diagonal (largest to smallest). For an mxn matrix,
the U returned is m by (min m,n), L is (min m,n) by (min m,n), and V is
m by (min m,n).
It also assures that the singular values are returned as real (if a
complex matrix is given, they are returned as complex with zero imaginary part
by sv-decomp).
Returns a list of the matrix U, a list of the singular values, the matrix
V, and a flag indicating either that the algorithm converged (T), or NIL
otherwise.
;; See Also: sv-decomp, pseudo-inverse
;; Examples:
(def a (matrix '(3 2) (list 1 -1 -1 2 1 2)))
(def svd-a (svd a))
;; as promised works on matrices where the number of rows is less than the
;; number of columns:
(def a-transpose (transpose a))
(def svd-a-transpose (svd a-transpose))
;; Here's how to put the matrix a back together again:
(def u (first svd-a))
(def lambda (second svd-a))
(def v (third svd-a))
(mprint (matmult u (diagonal lambda) (transpose v))) ;; = a
;; works for complex matrices, too:
(def a #2a((2 #c(0 1))(#c(0 1) 2)))
(def svd-a (svd a))
(def u (first svd-a))
(def lambda (second svd-a))
(def v (third svd-a))
;; Notice that to put a complex matrix back together you need to conjugate the
;; (possibly complex) matrix v :
(mprint (matmult u (diagonal lambda) (transpose (conjugate v))))
"
(if (matrixp a)
(let* ( ; true, so get the matrix dimensions:
(size (size a))
(swap (< (first size) (second size)))
;; if swap is true, then transpose:
(svd (if swap (sv-decomp (transpose a))
(sv-decomp a)))
(u (if swap (third svd) (first svd)))
;; make sure that the lambda are real numbers, and return them as
;; a list rather than as a vector:
(lambda (mapcar #'realpart (coerce (second svd) 'list)))
(v (if swap (first svd) (third svd)))
(converged (fourth svd))
)
(list u lambda v converged)
)
)
)
(defun mprint(A &rest args)
(apply
#'print-matrix
(if (matrixp A) A (bind-columns A))
args)
)