Singular value decomposition: Difference between revisions

Content added Content deleted
m (→‎{{header|C}}: V -> VT in output.)
(→‎{{header|Wren}}: Replaced existing solution with one which uses GSL.)
Line 240: Line 240:


=={{header|Wren}}==
=={{header|Wren}}==
{{trans|C}}
{{libheader|Wren-matrix}}
{{libheader|GNU Scientific Library}}
{{libheader|Wren-fmt}}
{{libheader|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.
We don't have a built-in method for this so we need to work from first principles.
<syntaxhighlight lang="ecmascript">/* svd_emdedded.wren */
The following only decomposes 2 x 2 matrices.

<syntaxhighlight lang="ecmascript">import "./matrix" for Matrix
import "./fmt" for Fmt
import "./fmt" for Fmt


var svd2x2 = Fn.new { |a, width, prec|
var matrixPrint = Fn.new { |r, c, m|
var at = a.transpose
for (i in 0...r) {
var p = a * at
System.write("|")
for (j in 0...c) {
// the eigenvalues, λ, are the solutions of the characteristic polynomial:
// λ^2 - (p[0,0] + p[1,1]) * λ + p[0,0]*p[1,1] - p[0,1]*p[1,0] = 0
Fmt.write("$13.10f ", m[i*c + j])
var b = -(p[0, 0] + p[1, 1])
}
System.print("\b|")
var c = p[0, 0]*p[1, 1] - p[0, 1]*p[1, 0]
}
var d = (b * b - 4 * c).sqrt
var l1 = (-b + d) / 2
System.print()
}
var l2 = (-b - d) / 2

// hence the singular values and 'sigma' are:
class GSL {
var s1 = l1.sqrt
// returns 'v' in transposed form
var s2 = l2.sqrt
var sigma = Matrix.new([[s1, 0], [0, s2]])
foreign static svd(r, c, a, v, s)
}
// now find the eigenvectors

p = at * a
var r = 2
var i = Matrix.identity(p.numRows)
var m1 = p - i*l1
var c = 2
var m2 = p - i*l2
var l = r * c
var a = [3, 0, 4, 5]
m1.toReducedRowEchelonForm
var v = List.filled(l, 0)
m2.toReducedRowEchelonForm
var s = List.filled(l, 0)
// solve for null space and convert to unit vector

var ev1 = Matrix.new([[-m1[0,0] / m1[0, 1]], [1]])
GSL.svd(r, c, a, v, s)
ev1 = ev1 / ev1.norm
System.print("U:")
var ev2 = Matrix.new([[-m2[0,0] / m2[0, 1]], [1]])
matrixPrint.call(r, c, a)
ev2 = ev2 / ev2.norm
System.print("Σ:")
// we can now put these vectors together to get 'v' transpose
matrixPrint.call(r, c, s)
var vt = Matrix.new([[ev1[0, 0], ev2[0, 0]], [ev1[1, 0], ev2[1, 0]]])
System.print("VT:")
// and since 'v' is orthogonal its transpose and inverse are the same
matrixPrint.call(r, c, v)</syntaxhighlight>
var u = a * vt * (sigma.inverse)

System.print("U:")
We now embed this Wren script in the following C program, compile and run it.
Fmt.mprint(u, width, prec)
<syntaxhighlight lang="c">/* gcc svd_embedded.c -o svd_embedded -lgsl -lgslcblas -lwren -lm */
System.print("\nΣ:")

Fmt.mprint(sigma, width, prec)
#include <stdio.h>
System.print("\nVT:")
#include <string.h>
Fmt.mprint(vt, width, prec)
#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) {
var a = Matrix.new([[3, 0], [4, 5]])
WrenConfiguration config;
svd2x2.call(a, 17, 14)</syntaxhighlight>
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;
}</syntaxhighlight>


{{out}}
{{out}}
<pre>
<pre>
U:
U:
| 0.31622776601684 -0.94868329805051|
|-0.3162277660 -0.9486832981|
| 0.94868329805051 0.31622776601684|
|-0.9486832981 0.3162277660|


Σ:
Σ:
| 6.70820393249937 0.00000000000000|
| 6.7082039325 0.0000000000|
| 0.00000000000000 2.23606797749979|
| 0.0000000000 2.2360679775|


VT:
VT:
| 0.70710678118655 -0.70710678118655|
|-0.7071067812 -0.7071067812|
| 0.70710678118655 0.70710678118655|
|-0.7071067812 0.7071067812|
</pre>
</pre>
Retrieved from "https://rosettacode.org/wiki/Singular_value_decomposition"