Singular value decomposition
You are encouraged to solve this task according to the task description, using any language you may know.
is any m by n matrix, square or rectangular. Its rank is r. We will diagonalize this A, but not by Failed to parse (syntax error): {\displaystyle X^{−1}AX} . The eigenvectors in have three big problems: They are usually not orthogonal, there are not always enough eigenvectors, and = Failed to parse (syntax error): {\displaystyle λx} requires to be a square matrix. The singular vectors of solve all those problems in a perfect way.
The Singular Value Decomposition (SVD)
According to the web page above, for any rectangular matrix , we can decomposite it as Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle A=UΣV^T}
Task Description
Firstly, input two numbers "m" and "n".
Then, input a square/rectangular matrix .
Finally, output Failed to parse (syntax error): {\displaystyle U,Σ,V} with respect to .
Example
Sample Input
2 2 3 0 4 5
From the input above we can know that is a 2 by 2 matrix.
Sample Output
0.31622776601683794 -0.9486832980505138 0.9486832980505138 0.31622776601683794 6.708203932499369 0 0 2.23606797749979 0.7071067811865475 -0.7071067811865475 0.7071067811865475 0.7071067811865475
The output may vary depending your choice of the data types.
Remark
1. It’s encouraged to implement the algorithm by yourself while using libraries is still acceptible.
2. The algorithm should be applicable for general case().
C
The gsl_linalg_SV_decomp function can decompose any m x n matrix though the example here is for a 2 x 2 matrix.
Requires a C99 or later compiler.
#include <stdio.h>
#include <gsl/gsl_linalg.h>
/* Custom function for printing a gsl_matrix in matrix form. */
void gsl_matrix_print(const gsl_matrix *M) {
int rows = M->size1;
int cols = M->size2;
for (int i = 0; i < rows; i++) {
printf("|");
for (int j = 0; j < cols; j++) {
printf("% 12.10f ", gsl_matrix_get(M, i, j));
}
printf("\b|\n");
}
printf("\n");
}
int main(){
double a[] = {3, 0, 4, 5};
gsl_matrix_view A = gsl_matrix_view_array(a, 2, 2);
gsl_matrix *V = gsl_matrix_alloc(2, 2);
gsl_vector *S = gsl_vector_alloc(2);
gsl_vector *work = gsl_vector_alloc(2);
/* V is returned here in untransposed form. */
gsl_linalg_SV_decomp(&A.matrix, V, S, work);
gsl_matrix_transpose(V);
double s[] = {S->data[0], 0, 0, S->data[1]};
gsl_matrix_view SM = gsl_matrix_view_array(s, 2, 2);
printf("U:\n");
gsl_matrix_print(&A.matrix);
printf("S:\n");
gsl_matrix_print(&SM.matrix);
printf("VT:\n");
gsl_matrix_print(V);
gsl_matrix_free(V);
gsl_vector_free(S);
gsl_vector_free(work);
return 0;
}
- Output:
U: |-0.3162277660 -0.9486832981| |-0.9486832981 0.3162277660| S: | 6.7082039325 0.0000000000| | 0.0000000000 2.2360679775| VT: |-0.7071067812 -0.7071067812| |-0.7071067812 0.7071067812|
C++
Requires C++11 or higher.
#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;
}
- Output:
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 |
Alternative using Eigen
#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;
}
- Output:
0.316227766016838 0.948683298050514 0.948683298050514 -0.316227766016838 6.708203932499368 2.236067977499789 0.707106781186547 0.707106781186548 0.707106781186548 -0.707106781186547
Go
The SVD.Factorize method can decompose any m x n matrix though the example here is for a 2 x 2 matrix.
<package main
import (
"fmt"
"gonum.org/v1/gonum/mat"
"log"
)
func matPrint(m mat.Matrix) {
fa := mat.Formatted(m, mat.Prefix(""), mat.Squeeze())
fmt.Printf("%13.10f\n", fa)
}
func main() {
var svd mat.SVD
a := mat.NewDense(2, 2, []float64{3, 0, 4, 5})
ok := svd.Factorize(a, mat.SVDFull)
if !ok {
log.Fatal("Something went wrong!")
}
u := mat.NewDense(2, 2, nil)
svd.UTo(u)
fmt.Println("U:")
matPrint(u)
values := svd.Values(nil)
sigma := mat.NewDense(2, 2, []float64{values[0], 0, 0, values[1]})
fmt.Println("\nΣ:")
matPrint(sigma)
vt := mat.NewDense(2, 2, nil)
svd.VTo(vt)
fmt.Println("\nVT:")
matPrint(vt)
}
- Output:
U: ⎡-0.3162277660 -0.9486832981⎤ ⎣-0.9486832981 0.3162277660⎦ Σ: ⎡6.7082039325 0.0000000000⎤ ⎣0.0000000000 2.2360679775⎦ VT: ⎡-0.7071067812 -0.7071067812⎤ ⎣-0.7071067812 0.7071067812⎦
J
Using https://github.com/jsoftware/math_misc and providing a test matrix inline.
Assumed pre-requisite:
install'all'
Task example:
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│
└──────────────────┴──────────────┴──────────────────┘
Asymmetric example:
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│ │ │
└──────────────────────────────┴────────────────────────┴─────────────────────────────┘
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:
({.,(* =@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│ │ │
└──────────────────────────────┴────────────────────────┴─────────────────────────────┘
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.
Java
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"
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);
}
}
- Output:
0.3162277660 0.9486832981 0.9486832981 -0.3162277660 6.7082039325 0.0000000000 0.0000000000 2.2360679775 0.7071067812 0.7071067812 0.7071067812 -0.7071067812
Julia
Julia has an svd() function as part of its built-in LinearAlgebra package.
julia> using LinearAlgebra
julia> function testsvd()
rows, cols = [parse(Int, s) for s in split(readline())]
arr = zeros(rows, cols)
for row in 1:rows
arr[row, :] .= [tryparse(Float64, s) for s in split(readline())]
end
display(svd(arr))
end
testsvd (generic function with 1 method)
julia> testsvd()
2 2
3 0
4 5
SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}
U factor:
2×2 Matrix{Float64}:
-0.316228 -0.948683
-0.948683 0.316228
singular values:
2-element Vector{Float64}:
6.70820393249937
2.2360679774997894
Vt factor:
2×2 Matrix{Float64}:
-0.707107 -0.707107
-0.707107 0.707107
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)
- Output:
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|
Phix
with javascript_semantics requires("1.0.2") -- builtins/svd.e added, haha include builtins/svd.e sequence a = {{3,0}, {4,5}} sequence {u,w,v} = svd(a) ?u ?w ?v
- Output:
{{0.9486832981,0.316227766},{-0.316227766,0.9486832981}}
{2.236067977,6.708203932}
{{0.7071067812,0.7071067812},{-0.7071067812,0.7071067812}}
Python
from numpy import *
A = matrix([[3, 0], [4, 5]])
U, Sigma, VT = linalg.svd(A)
print(U)
print(Sigma)
print(VT)
- Output:
[[-0.31622777 -0.9486833 ] [-0.9486833 0.31622777]] [6.70820393 2.23606798] [[-0.70710678 -0.70710678] [-0.70710678 0.70710678]]
Raku
# 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' }
- Output:
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])
RPL
Singular values decomposition has a built-in instruction in RPL versions running on HP-48G and newer models.
[[ 3 0 ][ 4 5 ]] SVD
- Output:
3: [[-0.316227766017 -0.948683298051][-0.948683298051 0.316227766017]] 2: [[-0.707106781187 -0.707106781187][-0.707106781187 -0.707106781187]] 1: [6.70820393245 2.2360679775]
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 "~%")))
- Output:
$ 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
Wren
An embedded solution so we can use GSL to perform SVD on any m x n matrix though the example here is for a 2 x 2 matrix.
/* Singular_value_decomposition.wren */
import "./fmt" for Fmt
var matrixPrint = Fn.new { |r, c, m|
for (i in 0...r) {
System.write("|")
for (j in 0...c) {
Fmt.write("$13.10f ", m[i*c + j])
}
System.print("\b|")
}
System.print()
}
class GSL {
// returns 'v' in transposed form
foreign static svd(r, c, a, v, s)
}
var r = 2
var c = 2
var l = r * c
var a = [3, 0, 4, 5]
var v = List.filled(l, 0)
var s = List.filled(l, 0)
GSL.svd(r, c, a, v, s)
System.print("U:")
matrixPrint.call(r, c, a)
System.print("Σ:")
matrixPrint.call(r, c, s)
System.print("VT:")
matrixPrint.call(r, c, v)
We now embed this Wren script in the following C program, compile and run it.
/* gcc Singular_value_decomposition.c -o Singular_value_decomposition -lgsl -lgslcblas -lwren -lm */
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_linalg.h>
#include "wren.h"
void C_svd(WrenVM* vm) {
int r = (int)wrenGetSlotDouble(vm, 1);
int c = (int)wrenGetSlotDouble(vm, 2);
int l = r * c;
double a[l];
for (int i = 0; i < l; ++i) {
wrenGetListElement(vm, 3, i, 1);
a[i] = wrenGetSlotDouble(vm, 1);
}
gsl_matrix_view A = gsl_matrix_view_array(a, r, c);
gsl_matrix *V = gsl_matrix_alloc(r, c);
gsl_vector *S = gsl_vector_alloc(r);
gsl_vector *work = gsl_vector_alloc(r);
/* V is returned here in untransposed form. */
gsl_linalg_SV_decomp(&A.matrix, V, S, work);
gsl_matrix_transpose(V);
double s[] = {S->data[0], 0, 0, S->data[1]};
gsl_matrix_view SM = gsl_matrix_view_array(s, 2, 2);
for (int i = 0; i < r; ++i) {
for (int j = 0; j < c; ++j) {
int ix = i*c + j;
wrenSetSlotDouble(vm, 1, gsl_matrix_get(&A.matrix, i, j));
wrenSetListElement(vm, 3, ix, 1);
wrenSetSlotDouble(vm, 1, gsl_matrix_get(V, i, j));
wrenSetListElement(vm, 4, ix, 1);
wrenSetSlotDouble(vm, 1, gsl_matrix_get(&SM.matrix, i, j));
wrenSetListElement(vm, 5, ix, 1);
}
}
gsl_matrix_free(V);
gsl_vector_free(S);
gsl_vector_free(work);
}
WrenForeignMethodFn bindForeignMethod(
WrenVM* vm,
const char* module,
const char* className,
bool isStatic,
const char* signature) {
if (strcmp(module, "main") == 0) {
if (strcmp(className, "GSL") == 0) {
if (isStatic && strcmp(signature, "svd(_,_,_,_,_)") == 0) return C_svd;
}
}
return NULL;
}
static void writeFn(WrenVM* vm, const char* text) {
printf("%s", text);
}
void errorFn(WrenVM* vm, WrenErrorType errorType, const char* module, const int line, const char* msg) {
switch (errorType) {
case WREN_ERROR_COMPILE:
printf("[%s line %d] [Error] %s\n", module, line, msg);
break;
case WREN_ERROR_STACK_TRACE:
printf("[%s line %d] in %s\n", module, line, msg);
break;
case WREN_ERROR_RUNTIME:
printf("[Runtime Error] %s\n", msg);
break;
}
}
char *readFile(const char *fileName) {
FILE *f = fopen(fileName, "r");
fseek(f, 0, SEEK_END);
long fsize = ftell(f);
rewind(f);
char *script = malloc(fsize + 1);
fread(script, 1, fsize, f);
fclose(f);
script[fsize] = 0;
return script;
}
static void loadModuleComplete(WrenVM* vm, const char* module, WrenLoadModuleResult result) {
if( result.source) free((void*)result.source);
}
WrenLoadModuleResult loadModule(WrenVM* vm, const char* name) {
WrenLoadModuleResult result = {0};
if (strcmp(name, "random") != 0 && strcmp(name, "meta") != 0) {
result.onComplete = loadModuleComplete;
char fullName[strlen(name) + 6];
strcpy(fullName, name);
strcat(fullName, ".wren");
result.source = readFile(fullName);
}
return result;
}
int main(int argc, char **argv) {
WrenConfiguration config;
wrenInitConfiguration(&config);
config.writeFn = &writeFn;
config.errorFn = &errorFn;
config.bindForeignMethodFn = &bindForeignMethod;
config.loadModuleFn = &loadModule;
WrenVM* vm = wrenNewVM(&config);
const char* module = "main";
const char* fileName = "Singular_value_decomposition.wren";
char *script = readFile(fileName);
WrenInterpretResult result = wrenInterpret(vm, module, script);
switch (result) {
case WREN_RESULT_COMPILE_ERROR:
printf("Compile Error!\n");
break;
case WREN_RESULT_RUNTIME_ERROR:
printf("Runtime Error!\n");
break;
case WREN_RESULT_SUCCESS:
break;
}
wrenFreeVM(vm);
free(script);
return 0;
}
- Output:
U: |-0.3162277660 -0.9486832981| |-0.9486832981 0.3162277660| Σ: | 6.7082039325 0.0000000000| | 0.0000000000 2.2360679775| VT: |-0.7071067812 -0.7071067812| |-0.7071067812 0.7071067812|