|
[[-0.70710678 -0.70710678]
[-0.70710678 0.70710678]]
</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>
|