Rosetta Code

Singular value decomposition: Difference between revisions

Page Discussion

Singular value decomposition (view source)

Revision as of 22:08, 8 December 2022

3,180 bytes added ,  1 year ago
→‎{{header|Wren}}: Replaced existing solution with one which uses GSL.
m (→‎{{header|C}}: V -> VT in output.)
(→‎{{header|Wren}}: Replaced existing solution with one which uses GSL.)
Line 240:
 
=={{header|Wren}}==
{{trans|C}}
{{libheader|Wren-matrix}}
{{libheader|GNU Scientific Library}}
{{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
 
var svd2x2matrixPrint = Fn.new { |ar, widthc, precm|
varfor at(i =in a0.transpose..r) {
var p = a * atSystem.write("|")
for (j in 0...c) {
// the eigenvalues, λ, are the solutions of the characteristic polynomial:
// λ^2 - (p[0,0] + p[1,1]) * λ +Fmt.write("$13.10f p[0",0]*p[1,1] - pm[0,1]i*p[1,0]c =+ 0j])
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 = System.print(-b + d) / 2
}
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
varforeign sigmastatic = Matrix.newsvd([[s1r, c, 0]a, [0v, s2]]s)
}
// now find the eigenvectors
 
p = at * a
var r = 2
var i = Matrix.identity(p.numRows)
var m1c = p - i*l12
var m2l = pr -* i*l2c
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}}
<pre>
U:
| -0.316227766016843162277660 -0.948683298050519486832981|
| -0.948683298050519486832981 0.316227766016843162277660|
 
Σ:
| 6.708203932499377082039325 0.000000000000000000000000|
| 0.000000000000000000000000 2.236067977499792360679775|
 
VT:
| -0.707106781186557071067812 -0.707106781186557071067812|
| -0.707106781186557071067812 0.707106781186557071067812|
</pre>
PureFox
9,476

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/334168"