Singular value decomposition

From Rosetta Code
Task
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 . The eigenvectors in have three big problems: They are usually not orthogonal, there are not always enough eigenvectors, and = 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

Task Description

Firstly, input two numbers "m" and "n".

Then, input a square/rectangular matrix .

Finally, output 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[edit]

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++[edit]

Translation of: 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  |

Go[edit]

Library: gonum

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⎦

Java[edit]

Works with: Java version 10+
Library: 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"

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[edit]

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

Phix[edit]

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[edit]

Library: numpy
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[edit]

# 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])

Wren[edit]

Translation of: C
Library: Wren-fmt

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.

/* svd_emdedded.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 svd_embedded.c -o svd_embedded -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 = "svd_embedded.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|