Rosetta Code

Singular value decomposition: Difference between revisions

Page Discussion

Singular value decomposition (view source)

Revision as of 12:09, 6 February 2024

34,182 bytes added ,  3 months ago
m
→‎{{header|Wren}}: Minor tidy
(→‎{{header|Wren}}: Replaced existing solution with one which uses GSL.)
m (→‎{{header|Wren}}: Minor tidy)
 
(15 intermediate revisions by 8 users not shown)
Line 113:
|-0.7071067812 -0.7071067812|
|-0.7071067812 0.7071067812|
</pre>
 
=={{header|C++}}==
{{trans|C}}
{{libheader|GNU Scientific Library}}
Requires C++11 or higher.
<syntaxhighlight lang="C++">
#include <iostream>
#include <iomanip>
#include <gsl/gsl_linalg.h>
 
void gsl_matrix_print(const gsl_matrix *M) {
auto rows = M->size1;
auto cols = M->size2;
std::cout.precision(10);
std::cout.setf(std::ios::fixed);
for (auto i = 0; i < rows; i++) {
std::cout << "|";
for (auto j = 0; j < cols; j++) {
std::cout << gsl_matrix_get(M, i, j) << " ";
}
std::cout << "|" << std::endl;
}
std::cout << std::endl;
}
 
int main(int argc, char** argv) {
double a[] = {3, 0, 4, 5};
auto A = gsl_matrix_view_array(a, 2, 2);
auto *V = gsl_matrix_alloc(2, 2);
auto *S = gsl_vector_alloc(2);
auto *work = gsl_vector_alloc(2);
 
gsl_linalg_SV_decomp(&A.matrix, V, S, work);
double s[] = {S->data[0], 0.0, 0.0, S->data[1]};
auto SM = gsl_matrix_view_array(s, 2, 2);
std::cout << "U:" << std::endl;
gsl_matrix_print(&A.matrix);
 
std::cout << "Sigma:" << std::endl;
gsl_matrix_print(&SM.matrix);
 
std::cout << "V:" << std::endl;
gsl_matrix_print(V);
 
gsl_matrix_free(V);
gsl_vector_free(S);
gsl_vector_free(work);
return 0;
}
</syntaxhighlight>
 
{{out}}
<pre>
U:
|-0.3162277660 -0.9486832981 |
|-0.9486832981 0.3162277660 |
 
Sigma:
|6.7082039325 0.0000000000 |
|0.0000000000 2.2360679775 |
 
V:
|-0.7071067812 -0.7071067812 |
|-0.7071067812 0.7071067812 |
 
</pre>
 
===Alternative using Eigen===
{{libheader|Eigen}}
<syntaxhighlight lang="C++">
#include <Eigen/Dense>
#include <cstdlib>
#include <iomanip>
#include <iostream>
 
int main() {
using namespace Eigen;
int rows, columns;
if (!(std::cin >> rows >> columns))
return EXIT_FAILURE;
MatrixXd matrix(rows, columns);
for (int row = 0; row < rows; ++row) {
for (int column = 0; column < columns; ++column) {
if (!(std::cin >> matrix(row, column)))
return EXIT_FAILURE;
}
}
auto svd = matrix.bdcSvd().compute(matrix, ComputeFullU | ComputeFullV);
std::cout << std::setprecision(15) << std::fixed << svd.matrixU() << "\n\n"
<< svd.singularValues() << "\n\n"
<< svd.matrixV() << '\n';
return EXIT_SUCCESS;
}
</syntaxhighlight>
 
{{out}}
<pre>
0.316227766016838 0.948683298050514
0.948683298050514 -0.316227766016838
 
6.708203932499368
2.236067977499789
 
0.707106781186547 0.707106781186548
0.707106781186548 -0.707106781186547
</pre>
 
Line 165 ⟶ 271:
⎡-0.7071067812 -0.7071067812⎤
⎣-0.7071067812 0.7071067812⎦
</pre>
 
=={{header|J}}==
 
Using https://github.com/jsoftware/math_misc and providing a test matrix inline.
 
Assumed pre-requisite:
 
<syntaxhighlight lang=J> install'all'</syntaxhighlight>
 
Task example:
 
<syntaxhighlight lang=J> require'math/misc/svd'
svd 3 0,:4 5
┌──────────────────┬──────────────┬──────────────────┐
│0.316228 _0.948683│6.7082 2.23607│0.707107 _0.707107│
│0.948683 0.316228│ │0.707107 0.707107│
└──────────────────┴──────────────┴──────────────────┘</syntaxhighlight>
 
Asymmetric example:
 
<syntaxhighlight lang=J> svd 2 3 5,7 11 13,17 19 23,:29 31 37
┌──────────────────────────────┬────────────────────────┬─────────────────────────────┐
│0.0872159 _0.426839 _0.875752│68.6864 2.90306 0.868048│0.499426 0.860756 _0.0983486│
│ 0.265687 _0.830077 0.466712│ │0.554592 _0.230425 0.799582│
│ 0.499894 _0.0635067 _0.12233│ │0.665583 _0.453876 _0.592449│
│ 0.819701 0.353195 0.0165088│ │ │
└──────────────────────────────┴────────────────────────┴─────────────────────────────┘</syntaxhighlight>
 
Note that Σ is just the diagonal here. If we wish it in matrix form we can multiply that diagonal by the corresponding identity matrix. For example:
 
<syntaxhighlight lang=J> ({.,(* =@i.@#)&.>@(1&{),{:) svd 2 3 5,7 11 13,17 19 23,:29 31 37
┌──────────────────────────────┬────────────────────────┬─────────────────────────────┐
│0.0872159 _0.426839 _0.875752│68.6864 0 0│0.499426 0.860756 _0.0983486│
│ 0.265687 _0.830077 0.466712│ 0 2.90306 0│0.554592 _0.230425 0.799582│
│ 0.499894 _0.0635067 _0.12233│ 0 0 0.868048│0.665583 _0.453876 _0.592449│
│ 0.819701 0.353195 0.0165088│ │ │
└──────────────────────────────┴────────────────────────┴─────────────────────────────┘</syntaxhighlight>
 
<span style="background-color: black">File i/o is also implementable, though the current task specification (where the input matrix format (explicit shape) is different from the output matrix format (implicit shape)) seems unnecessary for svd.</span>
 
=={{header|Java}}==
{{works with|Java|10+}}
{{libheader|Jama}}
The library "Jama" can decompose m x n matrix though the example here is a 2 x 2 matrix.
 
If you want to use lower Java version, you need to replace the first "var" with "Matrix" and the second with "Singular
ValueDecomposition". You need also import "Jama.SingularValueDecomposition" and change the name of "public class"
<syntaxhighlight lang="java">
import Jama.Matrix;
public class SingularValueDecomposition {
public static void main(String[] args) {
double[][] matrixArray = {{3, 0}, {4, 5}};
var matrix = new Matrix(matrixArray);
var svd = matrix.svd();
svd.getU().print(0, 10); // The number of digits after the decimal is 10.
svd.getS().print(0, 10);
svd.getV().print(0, 10);
}
}
</syntaxhighlight>
{{out}}
<pre>
 
0.3162277660 0.9486832981
0.9486832981 -0.3162277660
 
 
6.7082039325 0.0000000000
0.0000000000 2.2360679775
 
 
0.7071067812 0.7071067812
0.7071067812 -0.7071067812
 
</pre>
 
Line 200 ⟶ 381:
-0.707107 0.707107
</syntaxhighlight>
 
=={{header|Nim}}==
{{libheader|arraymancer}}
<syntaxhighlight lang="Nim">import arraymancer, arraymancer/linear_algebra
 
var m = [[3, 0], [4, 5]].toTensor().asType(float)
let (u, s, vt) = m.svd()
 
# With "$", floats are displayed with 6 digits.
# So we use "pretty" to display 8 digits.
 
echo "U:"
echo u.pretty(8), '\n'
 
echo "Σ:"
echo s.pretty(8), '\n'
 
echo "V:"
echo vt.transpose().pretty(8)
</syntaxhighlight>
 
{{out}}
<pre>U:
Tensor[system.float] of shape "[2, 2]" on backend "Cpu"
|-0.31622777 -0.94868330|
|-0.94868330 0.31622777|
 
Σ:
Tensor[system.float] of shape "[2]" on backend "Cpu"
6.7082039 2.2360680
 
V:
Tensor[system.float] of shape "[2, 2]" on backend "Cpu"
|-0.70710678 -0.70710678|
|-0.70710678 0.70710678|
</pre>
 
=={{header|Phix}}==
Line 237 ⟶ 454:
[[-0.70710678 -0.70710678]
[-0.70710678 0.70710678]]
</pre>
 
=={{header|Raku}}==
<syntaxhighlight lang="raku" line># 20230108 Raku programming solution
 
use Math::Libgsl::Matrix;
use Math::Libgsl::LinearAlgebra;
 
my @M = <3 0>, <4 5>;
my Math::Libgsl::Matrix \M .= new: @M.elems, @M.first.elems;
(^M.size1)>>.&{ M.set-row: $_, @M[$_;*] }
 
my (\V,\S) = SV-decomp M;
 
say "U factor: " and say (^M.size1)>>.&{ M.get-row($_)>>.fmt: '%.10g' }
say "singular values: " and say (^S.size )>>.&{ S.get($_)>>.fmt: '%.10g' }
say "Vt factor: " and say (^V.size1)>>.&{ V.get-row($_)>>.fmt: '%.10g' }</syntaxhighlight>
{{out}}
<pre>U factor:
([-0.316227766 -0.9486832981] [-0.9486832981 0.316227766])
singular values:
((6.708203932) (2.236067977))
Vt factor:
([-0.7071067812 -0.7071067812] [-0.7071067812 0.7071067812])</pre>
 
=={{header|RPL}}==
Singular values decomposition has a built-in instruction in RPL versions running on HP-48G and newer models.
[[ 3 0 ][ 4 5 ]] SVD
{{out}}
<pre>
3: [[-0.316227766017 -0.948683298051][-0.948683298051 0.316227766017]]
2: [[-0.707106781187 -0.707106781187][-0.707106781187 -0.707106781187]]
1: [6.70820393245 2.2360679775]
</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>
 
Line 243 ⟶ 1,166:
{{libheader|GNU Scientific Library}}
{{libheader|Wren-fmt}}
An embedded solution so we can use GSL to perform SVD on any m x n matrix htoughthough the example here is for a 2 x 2 matrix.
<syntaxhighlight lang="ecmascriptwren">/* svd_emdeddedSingular_value_decomposition.wren */
 
import "./fmt" for Fmt
Line 280 ⟶ 1,203:
 
We now embed this Wren script in the following C program, compile and run it.
<syntaxhighlight lang="c">/* gcc svd_embeddedSingular_value_decomposition.c -o svd_embeddedSingular_value_decomposition -lgsl -lgslcblas -lwren -lm */
 
#include <stdio.h>
Line 393 ⟶ 1,316:
WrenVM* vm = wrenNewVM(&config);
const char* module = "main";
const char* fileName = "svd_embeddedSingular_value_decomposition.wren";
char *script = readFile(fileName);
WrenInterpretResult result = wrenInterpret(vm, module, script);
PureFox
9,476

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/334168...358315"