# Singular value decomposition

Singular value decomposition
You are encouraged to solve this task according to the task description, using any language you may know.

${\displaystyle A}$ is any m by n matrix, square or rectangular. Its rank is r. We will diagonalize this A, but not by $\displaystyle X^{−1}AX$ . The eigenvectors in ${\displaystyle X}$ have three big problems: They are usually not orthogonal, there are not always enough eigenvectors, and ${\displaystyle Ax}$ = $\displaystyle λx$ requires ${\displaystyle A}$ to be a square matrix. The singular vectors of ${\displaystyle A}$ solve all those problems in a perfect way.

According to the web page above, for any rectangular matrix ${\displaystyle A}$, we can decomposite it as $\displaystyle A=UΣV^T$

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

Then, input a square/rectangular matrix ${\displaystyle A^{m\times n}}$.

Finally, output $\displaystyle U,Σ,V$ with respect to ${\displaystyle A}$.

Example

Sample Input

2 2
3 0
4 5


From the input above we can know that ${\displaystyle A}$ 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(${\displaystyle m\times n}$).

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


## Go

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

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

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

with javascript_semantics
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

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


## Wren

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;
}
}

FILE *f = fopen(fileName, "r");
fseek(f, 0, SEEK_END);
long fsize = ftell(f);
rewind(f);
char *script = malloc(fsize + 1);
fclose(f);
script[fsize] = 0;
return script;
}

if( result.source) free((void*)result.source);
}

if (strcmp(name, "random") != 0 && strcmp(name, "meta") != 0) {
char fullName[strlen(name) + 6];
strcpy(fullName, name);
strcat(fullName, ".wren");
}
return result;
}

int main(int argc, char **argv) {
WrenConfiguration config;
wrenInitConfiguration(&config);
config.writeFn = &writeFn;
config.errorFn = &errorFn;
config.bindForeignMethodFn = &bindForeignMethod;
WrenVM* vm = wrenNewVM(&config);
const char* module = "main";
const char* fileName = "svd_embedded.wren";
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|