Svd


hedgehog
04-25-2009, 11:25 AM
Hi all - first post - be gentle ;)

I suspect my usage of svdcmp for mxn is wrong?

To test, I started with a simple mxm
2*x + 3*y = 18
5*x + 2*y = 23

x = 3, and y = 4.


float **a = matrix(1, 2, 1, 2);
float **v = matrix(1, 2, 1, 2);
int m = 2;
int n = 2;
float *w = vector(1, 2);

a[1][1] = 2;
a[1][2] = 3;
a[2][1] = 5;
a[2][2] = 2;

svdcmp(a, m, n, w, v);

as per NR: float wmin = wmax * 1.0e-6; etc. and zeroed out if necessary. N/A in this example.

float *b = vector(1, 2);
b[1] = 18.0f;
b[2] = 23.0f;

float *x = vector(1, 2);
x[1] = 0.0f;
x[2] = 0.0f;

svbksb(a, w, v, 2, 2, b, x);

The output for the x vector = (2.9999998, 3.9999986). So far so good...

In my application I need to solve mxn where m < n.
From my application: (there can be 10 - 20 eqs. where rows are always less than columns and RHS = 0).

2 * x – 4 * y + 0 * z = 0
3 * x + 0 * y - 7 * z = 0

I believe this is known as an underdetermined homogenous system and there are infinite solutions - E.g.

x = 14
y = 7
z = 6

x = 28
y = 14
z = 12

Using code:


float **a = matrix(1, 3, 1, 3);
float **v = matrix(1, 3, 1, 3);
int m = 3;
int n = 3;
float *w = vector(1, 3);

a[1][1] = 2;
a[1][2] = -4;
a[1][3] = 0;
a[2][1] = 3;
a[2][2] = 0;
a[2][3] = -7;
a[3][1] = 0;
a[3][2] = 0;
a[3][3] = 0;

svdcmp(a, m, n, w, v);

float wmax = 0.0f;

for (int ii = 1; ii <= 3; ii++)
if (w[ii] > wmax)
wmax = w[ii];

float wmin = wmax * 1.0e-6;

for (int ii = 1; ii <= 3; ii++)
if (w[ii] < wmin)
w[ii] = 0.0f;

float *b = vector(1, 3);
b[1] = 0.0f;
b[2] = 0.0f;
b[3] = 0.0f;

float *x = vector(1, 3);

svbksb(a, w, v, 3, 3, b, x);

/*
OUTPUT:

a[1][1] = -0.15234411
a[1][2] = 1.4359910e-025
a[1][3] = 0.98832738
a[2][1] = -0.98832750
a[2][2] = 9.3159454e-025
a[2][3] = -0.15234420
a[3][1] = -9.4259688e-025
a[3][2] = -0.99999994
a[3][3] = 0.00000000

w[1] = 7.6762538
w[2] = 2.0007453e-022 (prior to zeroing out)
w[3] = 4.3675098

v[1][1] = -0.42594624
v[1][2] = -0.83517003
v[1][3] = 0.34793806
v[2][1] = 0.079384640
v[2][2] = -0.41758484
v[2][3] = -0.90516347
v[3][1] = 0.90125895
v[3][2] = -0.35793003
v[3][3] = 0.24416845

x[1] = 0.00000000
x[2] = 0.00000000
x[3] = 0.00000000
*/


I'm really hoping some kind soul can take the time to tell me how to extract the solutions to this equation from NR, where I want to minimise x - i.e. x = 14 in the example given. Please help! C'mon - you know you want to! :p

Hedgehog.

hedgehog
04-27-2009, 05:24 AM
Okay, a bit more info:

2 * x – 4 * y + 0 * z = 0
3 * x + 0 * y - 7 * z = 0
-5 * x + 4 * y + 7 * z = 0
(linearly dependent?)

2 -4 0 x = 0
3 0 -7 y = 0
-5 4 7 z = 0

svdcmp gives me w[j] as before. So the 2nd column of v as before gives me a solution:

x = -0.83517003
y = -0.41758505
z = -0.35793000

Given that the solutions are a(x, y, z) where a = -16.76305362 would give x = 14, y = 7 and z = 6.

What I really need is to minimise x so that it is a positive whole number - in this case 14. y = 7, z = 6.

I'm really anxious and could really use your help. Even if you can post up just to say "you need to look in this area..." that will be good enough.

Thank you very much,

Hedgehog.

davekw7x
04-27-2009, 01:25 PM
I suspect my usage of svdcmp ... is wrong
I can't see where any techniques (svd or anything else) that involve floating point calculations can lead to the kind of integer solution that I think you require.

...if you can...say "you need to look in this area..."

Maybe start here:
http://en.wikipedia.org/wiki/Linear_programming
http://en.wikipedia.org/wiki/Linear_programming#Integer_unknowns

Regards,

Dave

hedgehog
04-28-2009, 05:05 AM
davekw7x,

That's extremely kind of you to point me in the right direction. Might be just what I'm looking for.

Many thanks again Dave,

Hedgehog.