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 |
var matrixPrint = Fn.new { |r, c, m| |
||
for (i in 0...r) { |
|||
System.write("|") |
|||
for (j in 0...c) { |
|||
// the eigenvalues, λ, are the solutions of the characteristic polynomial: |
|||
Fmt.write("$13.10f ", m[i*c + j]) |
|||
} |
|||
System.print("\b|") |
|||
var c = p[0, 0]*p[1, 1] - p[0, 1]*p[1, 0] |
|||
} |
|||
var d = (b * b - 4 * c).sqrt |
|||
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 |
|||
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 c = 2 |
|||
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.3162277660 -0.9486832981| |
||
| |
|-0.9486832981 0.3162277660| |
||
Σ: |
Σ: |
||
| 6. |
| 6.7082039325 0.0000000000| |
||
| 0. |
| 0.0000000000 2.2360679775| |
||
VT: |
VT: |
||
| |
|-0.7071067812 -0.7071067812| |
||
| |
|-0.7071067812 0.7071067812| |
||
</pre> |
</pre> |