Singular value decomposition

$A$ 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 $X$ have three big problems: They are usually not orthogonal, there are not always enough eigenvectors, and $Ax$ = Failed to parse (syntax error): {\displaystyle λx} requires $A$ to be a square matrix. The singular vectors of $A$ solve all those problems in a perfect way.

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

The Singular Value Decomposition (SVD)

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

Task Description

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

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

Finally, output Failed to parse (syntax error): {\displaystyle U,Σ,V} with respect to $A$.

Example

Sample Input

2 2
3 0
4 5

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

C

Library: GNU Scientific Library

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⎦

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

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: GNU Scientific Library

Library: Wren-fmt

An embedded solution so we can use GSL to perform SVD on any m x n matrix htough 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|