Singular value decomposition: Difference between revisions
Content added Content deleted
Tag: Manual revert |
|||
Line 438: | Line 438: | ||
Vt factor: |
Vt factor: |
||
([-0.7071067812 -0.7071067812] [-0.7071067812 0.7071067812])</pre> |
([-0.7071067812 -0.7071067812] [-0.7071067812 0.7071067812])</pre> |
||
=={{header|Scheme}}== |
|||
{{works with|CHICKEN Scheme|5.3.0}} |
|||
{{libheader|r7rs}} |
|||
{{libheader|srfi-42}} |
|||
{{libheader|srfi-132}} |
|||
{{libheader|srfi-143}} |
|||
{{libheader|srfi-144}} |
|||
{{libheader|srfi-179}} |
|||
{{libheader|format}} |
|||
<syntaxhighlight lang="scheme"> |
|||
;; This is code I wrote mostly in 2021. |
|||
;; |
|||
;; No doubt the code can be made faster for Scheme, by for instance |
|||
;; reducing the use of ‘set!’. That has sometimes seemed to help, at |
|||
;; least with some Scheme implementations. But the code should at |
|||
;; least work. |
|||
;; |
|||
;; If you are doing heavy duty vector processing, though, interface to |
|||
;; LAPACK instead. :) |
|||
;; |
|||
;; The Golub-Reinsch algorithm for the singular value decomposition, |
|||
;; based mostly on the EISPACK implementation. |
|||
;; |
|||
;; References: |
|||
;; |
|||
;; [1] G.H. Golub and C. Reinsch. 1970. Singular value decomposition |
|||
;; and least squares solutions. Numer. Math. 14 (1970), |
|||
;; 403–420. DOI: https://doi.org/10.1007/BF02163027 |
|||
;; |
|||
;; [2] EISPACK. https://www.netlib.org/eispack/ |
|||
;; |
|||
(define-library (svd) |
|||
(export svd) |
|||
(import (scheme base)) |
|||
;; This is part of R7RS-large, under the name (scheme sort). |
|||
(import (srfi 132)) ; Sorting. |
|||
;; These two are part of R7RS-large, under the names (scheme fixnum) |
|||
;; and (scheme flonum). |
|||
(import (srfi 143)) ; Fixnums. |
|||
(import (srfi 144)) ; Flonums. |
|||
;; |
|||
;; Arrays are not yet standardized in R7RS-large. I wrote this code |
|||
;; for CHICKEN Scheme, which supported SRFI-179 via an ‘egg’. So I |
|||
;; use SRFI-179. But SRFI-179 is not currently supported in Gauche |
|||
;; Scheme. Thus this code, as written, will not work in Gauche. |
|||
;; |
|||
;; Gauche DOES have ‘specialized arrays’, but they do not conform to |
|||
;; SRFI-179. At the time of this writing (10 May 2023), neither |
|||
;; CHICKEN nor Gauche supports SRFI-231 (a revision of SRFI-179), or |
|||
;; I would have switched to that. |
|||
;; |
|||
;; Arrays are on the docket for R7RS-large Orange edition: |
|||
;; https://github.com/johnwcowan/r7rs-work/blob/master/ColorDockets.md#user-content-orange-docket-portable-srfis |
|||
;; |
|||
(import (srfi 179)) ; Arrays. |
|||
;; SRFI-42 is not part of any Scheme standard, but is widely |
|||
;; supported. Both CHICKEN and Gauche support it. |
|||
(import (srfi 42)) ; Eager comprehensions. |
|||
(begin |
|||
(define (flpythag a b) |
|||
;; |
|||
;; Find sqrt(a**2+b**2) without overflow or destructive |
|||
;; underflow. |
|||
;; |
|||
;; Based on the EISPACK function PYTHAG(A,B). |
|||
;; |
|||
(let* ((abs-a (flabs a)) |
|||
(abs-b (flabs b)) |
|||
(p (flmax abs-a abs-b))) |
|||
(if (fl=? p 0.0) |
|||
p |
|||
(let ((sqrt-r (fl/ (flmin abs-a abs-b) p))) |
|||
(let loop ((p p) |
|||
(r (* sqrt-r sqrt-r))) |
|||
(let ((t (fl+ 4.0 r))) |
|||
(if (fl=? t 4.0) |
|||
p |
|||
(let* ((s (fl/ r t)) |
|||
(u (fl+ 1.0 (fl* 2.0 s))) |
|||
(p (fl* u p)) |
|||
(s/u (fl/ s u)) |
|||
(r (fl* (fl* s/u s/u) r))) |
|||
(loop p r))))))))) |
|||
(define svd-based-on-eispack-maximum-iterations |
|||
(make-parameter 30 (lambda (n) |
|||
(and (fixnum? n) (fxpositive? n))))) |
|||
(define (svd-based-on-eispack A with-U? with-V?) |
|||
;; |
|||
;; Computes the singular value decomposition A = UΣVᵀ of a real |
|||
;; matrix A, as implemented in EISPACK, where A and U are m×n, Σ |
|||
;; and V are n×n. The algorithm is based on that of Golub and |
|||
;; Reinsch. |
|||
;; |
|||
;; UᵀU = VᵀV = VVᵀ = I, the n×n identity matrix. |
|||
;; |
|||
;; Σ is a diagonal matrix of the singular values (the |
|||
;; non-negative square roots of the eigenvalues of AᵀA). In |
|||
;; addition, U comprises orthonormalized eigenvectors associated |
|||
;; with the n largest eigenvalues of AAᵀ, and V comprises the |
|||
;; orthonormalized eigenvectors of AᵀA. |
|||
;; |
|||
;; The input array A is assumed to be a two-dimensional SRFI-179 |
|||
;; array. The array is not modified, and the rows and columns |
|||
;; can be indexed starting at zero or one or however you like. |
|||
;; |
|||
;; So, for example, if there are to be m rows and n columns, |
|||
;; this will do: |
|||
;; |
|||
;; (make-specialized-array |
|||
;; (make-interval #(1 1) (vector (fx+ m 1) (fx+ n 1))) |
|||
;; f64-storage-class) |
|||
;; |
|||
;; The return values are |
|||
;; |
|||
;; (values ierr w U V) |
|||
;; |
|||
;; If ierr=0, then the procedure succeeded; otherwise ierr=k |
|||
;; where convergence failed while computing the kth singular |
|||
;; value. |
|||
;; |
|||
;; w is a single-dimensional f64-storage-class SRFI-179 array of |
|||
;; the singular values, indexed 1..n. It is the diagonal of |
|||
;; Σ. Note that the singular values are *not* sorted by size; |
|||
;; neither are they set to zero if below some threshold. |
|||
;; |
|||
;; If with-U? is #f, then U is #f. Otherwise it is a |
|||
;; two-dimensional f64-storage-class SRFI-179 array of the |
|||
;; matrix U, indexed 1..m,1..n. |
|||
;; |
|||
;; If with-V? is #f, then V is #f. Otherwise it is a |
|||
;; two-dimensional f64-storage-class SRFI-179 array of the |
|||
;; matrix V, indexed 1..n,1..n. |
|||
;; |
|||
(define (svd A m n) |
|||
(define (flsign a b) |
|||
;; Do what the SIGN function does in Fortran: if 0 <= b then |
|||
;; return abs(a), else -abs(a). |
|||
(if (fl<=? 0.0 b) |
|||
(flabs a) |
|||
(fl- (flabs a)))) |
|||
(define m+1 (fx+ m 1)) |
|||
(define n+1 (fx+ n 1)) |
|||
;; Some intervals. |
|||
(define 1:n (make-interval #(1) (vector n+1))) |
|||
(define 1:n×1:n (make-interval #(1 1) (vector n+1 n+1))) |
|||
;; The singular values. |
|||
(define w (make-specialized-array 1:n f64-storage-class)) |
|||
;; The left-hand result matrix, or workspace if the left-hand |
|||
;; matrix is not requested. |
|||
(define U (array-copy A f64-storage-class #f #t)) |
|||
;; The right-hand result matrix, if requested. |
|||
(define V (if with-V? |
|||
(make-specialized-array 1:n×1:n f64-storage-class) |
|||
#f)) |
|||
;; Some workspace. |
|||
(define rv1 (make-specialized-array 1:n f64-storage-class)) |
|||
;; Flonum variables. |
|||
(define c 0.0) |
|||
(define f 0.0) |
|||
(define g 0.0) |
|||
(define h 0.0) |
|||
(define s 0.0) |
|||
(define x 0.0) |
|||
(define y 0.0) |
|||
(define z 0.0) |
|||
(define scale 0.0) |
|||
;; The maximum number of iterations. |
|||
(define maxit (svd-based-on-eispack-maximum-iterations)) |
|||
;; Fixnum variables. |
|||
(define l 0) |
|||
(define (results ierr w U V) |
|||
(values ierr w (and with-U? U) (and with-V? V))) |
|||
(define (test-for-convergence iterations tst1 k l) |
|||
(define k1 (fx- k 1)) |
|||
;; |
|||
;; Test for convergence. |
|||
;; |
|||
(set! z (array-ref w k)) |
|||
(if (fx=? l k) |
|||
(begin |
|||
;; |
|||
;; Convergence succeeded. |
|||
;; |
|||
(unless (fl<=? 0.0 z) |
|||
;; w(k) is made non-negative. |
|||
(array-set! w (fl- z) k) |
|||
(when with-V? |
|||
(do-ec (:range j 1 n+1) |
|||
(array-set! V (fl- (array-ref V j k)) j k)))) |
|||
(diagonalization tst1 (fx- k 1))) |
|||
(begin |
|||
;; Shift from bottom 2-by-2 minor. |
|||
(cond |
|||
((fx=? iterations maxit) |
|||
;; |
|||
;; Convergence failed. |
|||
;; |
|||
(results k w U V)) |
|||
(else |
|||
(set! iterations (fx+ iterations 1)) |
|||
(set! x (array-ref w l)) |
|||
(set! y (array-ref w k1)) |
|||
(set! g (array-ref rv1 k1)) |
|||
(set! h (array-ref rv1 k)) |
|||
(set! f (fl* 0.5 |
|||
(fl- (fl+ (fl* (fl/ (fl+ g z) h) |
|||
(fl/ (fl- g z) y)) |
|||
(fl/ y h)) |
|||
(fl/ h y)))) |
|||
(set! g (flpythag f 1.0)) |
|||
(set! f (fl+ (fl- x (fl* (fl/ z x) z)) |
|||
(fl* (fl/ h x) |
|||
(fl- (fl/ y (fl+ f (flsign g f))) |
|||
h)))) |
|||
;; |
|||
;; The next QR transformation. |
|||
;; |
|||
(set! c 1.0) |
|||
(set! s 1.0) |
|||
(do-ec (:range i1 l k) |
|||
(let ((i (fx+ i1 1))) |
|||
(set! g (array-ref rv1 i)) |
|||
(set! y (array-ref w i)) |
|||
(set! h (fl* s g)) |
|||
(set! g (fl* c g)) |
|||
(set! z (flpythag f h)) |
|||
(array-set! rv1 z i1) |
|||
(set! c (fl/ f z)) |
|||
(set! s (fl/ h z)) |
|||
(set! f (fl+ (fl* g s) (fl* x c))) |
|||
(set! g (fl- (fl* g c) (fl* x s))) |
|||
(set! h (fl* y s)) |
|||
(set! y (fl* y c)) |
|||
(when with-V? |
|||
(do-ec (:range j 1 n+1) |
|||
(begin |
|||
(set! x (array-ref V j i1)) |
|||
(set! z (array-ref V j i)) |
|||
(array-set! V (fl+ (fl* z s) (fl* x c)) |
|||
j i1) |
|||
(array-set! V (fl- (fl* z c) (fl* x s)) |
|||
j i)))) |
|||
(set! z (flpythag f h)) |
|||
(array-set! w z i1) |
|||
(unless (fl=? z 0.0) |
|||
;; Rotation can be arbitrary if z is zero. |
|||
(set! c (fl/ f z)) |
|||
(set! s (fl/ h z))) |
|||
(set! f (fl+ (fl* s y) (fl* c g))) |
|||
(set! x (fl- (fl* c y) (fl* s g))) |
|||
(when with-U? |
|||
(do-ec (:range j 1 m+1) |
|||
(begin |
|||
(set! y (array-ref U j i1)) |
|||
(set! z (array-ref U j i)) |
|||
(array-set! U (fl+ (fl* z s) (fl* y c)) |
|||
j i1) |
|||
(array-set! U (fl- (fl* z c) (fl* y s)) |
|||
j i)))) )) |
|||
(array-set! rv1 0.0 l) |
|||
(array-set! rv1 f k) |
|||
(array-set! w x k) |
|||
(test-for-splitting iterations tst1 k)))))) |
|||
(define (cancellation iterations tst1 k l) |
|||
(define l1 (fx- l 1)) |
|||
;; |
|||
;; Cancellation of rv1(l) if l > 1. |
|||
;; |
|||
(set! c 0.0) |
|||
(set! s 1.0) |
|||
(do-ec (:range i l (fx+ k 1)) |
|||
(let ((rv1@i (array-ref rv1 i))) |
|||
(set! f (fl* s rv1@i)) |
|||
(array-set! rv1 (fl* c rv1@i) i) |
|||
(if (fl=? (fl+ tst1 (flabs f)) tst1) |
|||
(test-for-convergence iterations k l) |
|||
(begin |
|||
(set! g (array-ref w i)) |
|||
(set! h (flpythag f g)) |
|||
(array-set! w h i) |
|||
(set! c (fl/ g h)) |
|||
(set! s (fl/ (fl- f) h)) |
|||
(when with-U? |
|||
(do-ec (:range j 1 m+1) |
|||
(begin |
|||
(set! y (array-ref U j l1)) |
|||
(set! z (array-ref U j i)) |
|||
(array-set! U (fl+ (fl* z s) (fl* y c)) |
|||
j l1) |
|||
(array-set! U (fl- (fl* z c) (fl* y s)) |
|||
j i)))))))) |
|||
(test-for-convergence iterations tst1 k l)) |
|||
(define (test-for-splitting iterations tst1 k) |
|||
;; |
|||
;; Test for splitting. |
|||
;; |
|||
(let loop ((l k)) |
|||
(cond ((fx=? l 0) |
|||
;; rv1(1) is always zero; therefore l never reaches |
|||
;; zero. |
|||
(error 'svd-based-on-eispack |
|||
(string-append |
|||
"THERE IS A BUG: " |
|||
"this branch should never have been run"))) |
|||
((fl=? (fl+ tst1 (flabs (array-ref rv1 l))) |
|||
tst1) |
|||
(test-for-convergence iterations tst1 k l)) |
|||
((fl=? (fl+ tst1 (flabs (array-ref w (fx- l 1)))) |
|||
tst1) |
|||
(cancellation iterations tst1 k l)) |
|||
(else |
|||
(loop (fx- l 1)))))) |
|||
(define (diagonalization tst1 k) |
|||
;; |
|||
;; Diagonalization of the bidiagonal form. |
|||
;; |
|||
(if (fx=? k 0) |
|||
(results 0 w U V) |
|||
(test-for-splitting 0 tst1 k))) |
|||
;; |
|||
;; Householder reduction to bidiagonal form. |
|||
;; |
|||
(set! g 0.0) |
|||
(set! scale 0.0) |
|||
(set! x 0.0) |
|||
(do-ec (:range i 1 n+1) |
|||
(begin |
|||
(set! l (fx+ i 1)) |
|||
(array-set! rv1 (fl* scale g) i) |
|||
(set! g 0.0) |
|||
(set! s 0.0) |
|||
(set! scale 0.0) |
|||
(unless (fx<? m i) |
|||
(do-ec (:range k i m+1) |
|||
(set! scale (fl+ scale (flabs (array-ref U k i))))) |
|||
(unless (fl=? scale 0.0) |
|||
(do-ec (:range k i m+1) |
|||
(let ((U@ki (fl/ (array-ref U k i) scale))) |
|||
(array-set! U U@ki k i) |
|||
(set! s (fl+ s (fl* U@ki U@ki))))) |
|||
(set! f (array-ref U i i)) |
|||
(set! g (fl- (flsign (flsqrt s) f))) |
|||
(set! h (fl- (fl* f g) s)) |
|||
(array-set! U (fl- f g) i i) |
|||
(unless (fx=? i n) |
|||
(do-ec (:range j l n+1) |
|||
(begin |
|||
(set! s 0.0) |
|||
(do-ec (:range k i m+1) |
|||
(set! s (fl+ s (fl* (array-ref U k i) |
|||
(array-ref U k j))))) |
|||
(set! f (fl/ s h)) |
|||
(do-ec (:range k i m+1) |
|||
(array-set! U (fl+ (array-ref U k j) |
|||
(fl* f (array-ref U k i))) |
|||
k j)))) |
|||
) ;; end (unless (fx=? i n) ...) |
|||
(do-ec (:range k i m+1) |
|||
(array-set! U (fl* scale (array-ref U k i)) k i)) |
|||
) ;; end (unless (fl=? scale 0.0) ...) |
|||
) ;; end (unless (fx<? m i) ...) |
|||
(array-set! w (fl* scale g) i) |
|||
(set! g 0.0) |
|||
(set! s 0.0) |
|||
(set! scale 0.0) |
|||
(unless (or (fx<? m i) (fx=? i n)) |
|||
(do-ec (:range k l n+1) |
|||
(set! scale (fl+ scale (flabs (array-ref U i k))))) |
|||
(unless (fl=? scale 0.0) |
|||
(do-ec (:range k l n+1) |
|||
(let ((U@ik (fl/ (array-ref U i k) scale))) |
|||
(array-set! U U@ik i k) |
|||
(set! s (fl+ s (fl* U@ik U@ik))))) |
|||
(set! f (array-ref U i l)) |
|||
(set! g (fl- (flsign (flsqrt s) f))) |
|||
(set! h (fl- (fl* f g) s)) |
|||
(array-set! U (fl- f g) i l) |
|||
(do-ec (:range k l n+1) |
|||
(array-set! rv1 (fl/ (array-ref U i k) h) k)) |
|||
(unless (fx=? i m) |
|||
(do-ec (:range j l m+1) |
|||
(begin |
|||
(set! s 0.0) |
|||
(do-ec (:range k l n+1) |
|||
(set! s (fl+ s (fl* (array-ref U j k) |
|||
(array-ref U i k))))) |
|||
(do-ec (:range k l n+1) |
|||
(array-set! U (fl+ (array-ref U j k) |
|||
(fl* s (array-ref rv1 k))) |
|||
j k)))) |
|||
) ;; end (unless (fx=? i m) ... ) |
|||
(do-ec (:range k l n+1) |
|||
(array-set! U (fl* scale (array-ref U i k)) i k)) |
|||
) ;; end (unless (fl=? scale 0.0) ... ) |
|||
) ;; end (unless (or (fx<? m i) (fx=? i n)) ... ) |
|||
(set! x (flmax x (fl+ (flabs (array-ref w i)) |
|||
(flabs (array-ref rv1 i)))))) |
|||
) ;; end (do-ec (:range i 1 n+1) ... ) |
|||
;; |
|||
;; Accumulation of right-hand transformations. |
|||
;; |
|||
(when with-V? |
|||
(do-ec (:range i n 0 -1) |
|||
(begin |
|||
(unless (fx=? i n) |
|||
(unless (fl=? g 0.0) |
|||
(do-ec (:range j l n+1) |
|||
;; The double division avoids a possible underflow. |
|||
(array-set! V (fl/ (fl/ (array-ref U i j) |
|||
(array-ref U i l)) |
|||
g) |
|||
j i)) |
|||
(do-ec (:range j l n+1) |
|||
(begin |
|||
(set! s 0.0) |
|||
(do-ec (:range k l n+1) |
|||
(set! s (fl+ s (fl* (array-ref U i k) |
|||
(array-ref V k j))))) |
|||
(do-ec (:range k l n+1) |
|||
(array-set! V (fl+ (array-ref V k j) |
|||
(fl* s (array-ref V k i))) |
|||
k j))))) |
|||
(do-ec (:range j l n+1) |
|||
(begin |
|||
(array-set! V 0.0 i j) |
|||
(array-set! V 0.0 j i)))) |
|||
(array-set! V 1.0 i i) |
|||
(set! g (array-ref rv1 i)) |
|||
(set! l i)))) |
|||
;; |
|||
;; Accumulation of left-hand transformations. |
|||
;; |
|||
(when with-U? |
|||
(do-ec (:range i (fxmin m n) 0 -1) |
|||
(begin |
|||
(set! l (fx+ i 1)) |
|||
(set! g (array-ref w i)) |
|||
(unless (fx=? i n) |
|||
(do-ec (:range j l n+1) |
|||
(array-set! U 0.0 i j))) |
|||
(if (fl=? g 0.0) |
|||
(do-ec (:range j i m+1) |
|||
(array-set! U 0.0 j i)) |
|||
(begin |
|||
(unless (fx=? i (fxmin m n)) |
|||
(do-ec (:range j l n+1) |
|||
(begin |
|||
(set! s 0.0) |
|||
(do-ec (:range k l m+1) |
|||
(set! s (fl+ s (fl* (array-ref U k i) |
|||
(array-ref U k j))))) |
|||
;; The double division avoids a possible |
|||
;; underflow. |
|||
(set! f (fl/ (fl/ s (array-ref U i i)) g)) |
|||
(do-ec (:range k i m+1) |
|||
(array-set! |
|||
U (fl+ (array-ref U k j) |
|||
(fl* f (array-ref U k i))) |
|||
k j))) )) |
|||
(do-ec (:range j i m+1) |
|||
(array-set! U (fl/ (array-ref U j i) g) j i)))) |
|||
(array-set! U (fl+ (array-ref U i i) 1.0) i i) ))) |
|||
(diagonalization x n) |
|||
) ;; end of procedure svd. |
|||
(let* ((domain_A (array-domain A)) |
|||
(i^_min (interval-lower-bound domain_A 0)) |
|||
(i^_max+1 (interval-upper-bound domain_A 0)) |
|||
(j^_min (interval-lower-bound domain_A 1)) |
|||
(j^_max+1 (interval-upper-bound domain_A 1)) |
|||
;; m and n. |
|||
(m (fx- i^_max+1 i^_min)) |
|||
(n (fx- j^_max+1 j^_min)) |
|||
;; Reindex A, if necessary. |
|||
(A (if (and (fx=? i^_min 1) (fx=? j^_min 1)) |
|||
A |
|||
(let ((getter^ (array-getter A)) |
|||
(i^_min-1 (fx- i^_min 1)) |
|||
(j^_min-1 (fx- j^_min 1))) |
|||
(make-array |
|||
(make-interval #(1 1) (vector (fx+ m 1) |
|||
(fx+ n 1))) |
|||
(lambda (i j) |
|||
(getter^ (fx+ i i^_min-1) |
|||
(fx+ j j^_min-1)))))))) |
|||
(svd A m n))) |
|||
(define (sort-svd w U V) |
|||
;; |
|||
;; Sort w and the optional U and V in order of descending size |
|||
;; of the singular value. The reordered arrays are new (and |
|||
;; read-only), but share their storage space with the originals. |
|||
;; |
|||
;; FIXME: I tried to use ‘vector-sort!’ instead of |
|||
;; ‘vector-stable-sort!’, but, at the time of this |
|||
;; writing (9 Mar 2021), it seemed to be broken for the |
|||
;; compiler I was using. |
|||
;; |
|||
(let* ((n (- (interval-upper-bound (array-domain w) 0) 1)) |
|||
(n+1 (fx+ n 1)) |
|||
(indices (make-vector n+1))) |
|||
;; The first entry of ‘indices’ is ignore. Do not bother to set |
|||
;; it. Set the others to 1, 2, 3, ... , n. |
|||
(do-ec (:range j 1 n+1) |
|||
(vector-set! indices j j)) |
|||
(vector-stable-sort! (lambda (i j) |
|||
(let ((i^ (vector-ref indices i)) |
|||
(j^ (vector-ref indices j))) |
|||
(< (array-ref w j^) (array-ref w i^)))) |
|||
indices 1) |
|||
(values (make-array (array-domain w) |
|||
(lambda (j) |
|||
(let ((j^ (vector-ref indices j))) |
|||
(array-ref w j^)))) |
|||
(and U (make-array (array-domain U) |
|||
(lambda (i j) |
|||
(let ((j^ (vector-ref indices j))) |
|||
(array-ref U i j^))))) |
|||
(and V (make-array (array-domain V) |
|||
(lambda (i j) |
|||
(let ((j^ (vector-ref indices j))) |
|||
(array-ref V i j^)))))))) |
|||
(define (svd A with-U? with-V?) |
|||
;; |
|||
;; Compute the singular value decomposition (SVD) of the real |
|||
;; matrix stored in SRFI-179 array A. A is not altered. |
|||
;; |
|||
;; The indexing base of A may be 0, 1, some other number, or |
|||
;; different for the two dimensions. However, all the output |
|||
;; arrays are 1-indexed, read-only SRFI-179 arrays. |
|||
;; |
|||
;; Our definition of the singular value decomposition of a real |
|||
;; matrix A is A = UΣVᵀ where A and U are m×n, Σ and V are n×n. |
|||
;; |
|||
;; UᵀU = VᵀV = VVᵀ = I, the n×n identity matrix. |
|||
;; |
|||
;; Σ is a diagonal matrix of the singular values (the |
|||
;; non-negative square roots of the eigenvalues of AᵀA). In |
|||
;; addition, U comprises orthonormalized eigenvectors associated |
|||
;; with the n largest eigenvalues of AAᵀ, and V comprises the |
|||
;; orthonormalized eigenvectors of AᵀA. |
|||
;; |
|||
;; Return values: |
|||
;; |
|||
;; ierr : If the SVD was successfully computed, then ierr=0; |
|||
;; otherwise ierr=k, where convergence failed on the |
|||
;; kth singular value. |
|||
;; |
|||
;; w : The diagonal of Σ. In other words, the singular values. |
|||
;; The values will be sorted in descending order. |
|||
;; |
|||
;; U : If with-U? is false, then U=#f. Otherwise, it is the |
|||
;; matrix U. |
|||
;; |
|||
;; V : If with-V? is false, then V=#f. Otherwise, it is the |
|||
;; matrix V. |
|||
;; |
|||
(let-values (((ierr w U V) |
|||
(svd-based-on-eispack A with-U? with-V?))) |
|||
(if (zero? ierr) |
|||
(let-values (((w U V) (sort-svd w U V))) |
|||
(values ierr w U V)) |
|||
(values ierr w U V)))) |
|||
)) ;; end library. |
|||
(import (scheme base) |
|||
(scheme write) |
|||
(srfi 179) |
|||
(srfi 42) |
|||
(svd)) |
|||
;; Common Lisp-style formatting. |
|||
(import (format)) |
|||
(define m (read)) |
|||
(define n (read)) |
|||
(define A |
|||
(make-specialized-array |
|||
(make-interval #(1 1) (vector (+ m 1) (+ n 1))) |
|||
f64-storage-class)) |
|||
;; Read the elements of the input array. They are assumed to be in |
|||
;; row-major order. |
|||
(do-ec (nested (:range i 1 (+ m 1)) |
|||
(:range j 1 (+ n 1))) |
|||
(array-set! A (inexact (read)) i j)) |
|||
(define-values (ierr w U V) (svd A #t #t)) |
|||
(format #t "U:~%") |
|||
(do-ec (:range i 1 (+ m 1)) |
|||
(begin |
|||
(do-ec (:range j 1 (+ n 1)) |
|||
(format #t "~30,14F" (array-ref U i j))) |
|||
(format #t "~%"))) |
|||
(format #t "~%") |
|||
(format #t "Σ:~%") |
|||
(do-ec (:range i 1 (+ n 1)) |
|||
(begin |
|||
(do-ec (:range j 1 (+ n 1)) |
|||
(format #t "~30,14F" (if (= j i) (array-ref w i) 0.0))) |
|||
(format #t "~%"))) |
|||
(format #t "~%") |
|||
(format #t "V:~%") |
|||
(do-ec (:range i 1 (+ n 1)) |
|||
(begin |
|||
(do-ec (:range j 1 (+ n 1)) |
|||
(format #t "~30,14F" (array-ref V i j))) |
|||
(format #t "~%"))) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>$ csc -O5 -X r7rs -R r7rs svd_task.scm && echo '2 2 3 0 4 5' | ./svd_task |
|||
U: |
|||
-0.31622776601684 -0.94868329805051 |
|||
-0.94868329805051 0.31622776601684 |
|||
Σ: |
|||
6.70820393249937 0.00000000000000 |
|||
0.00000000000000 2.23606797749979 |
|||
V: |
|||
-0.70710678118655 -0.70710678118655 |
|||
-0.70710678118655 0.70710678118655 |
|||
</pre> |
|||
=={{header|Wren}}== |
=={{header|Wren}}== |