Hi Guys

My task is to run a regression (i.e. solve for a vector x from Ax = b) by using the SVD.

However, I already have calculated a design matrix A from my other routine (i.e. choose a basis function I want and obtain A accordingly) and also have a vector b ready.

Therefore, I skipped the first part of a class Fitsvd (actually I didn’t use Fitsvd at all) and started by constructing an SVD object (i.e. SVD mySVD(const & A)) and then called a solve member function (i.e. mySVD.solve(...)).

However, my result are incorrect (i.e. chisq = 100000000 or something!!!!). Hence my question is whether it is possible to skip the first part of Fitsvd by supplying a design matrix A directly for a decomposition and do the solve..etc. (To be honest, the line b[i] = y[i]*tmp in a Fitsvd worries me a lot because it seems there are some modification on the vector b while I supplied it directly).

Cheers

PChira

...the line b[i] = y[i]*tmp in a Fitsvd worries me a lot because it seems there are some modification on the vector b while I supplied it directly).

Actually, that does not happen.

The vector parameters to the constructor are xx, yy, and ssig.

In fitsvd.h:

Fitsvd(VecDoub_I &xx, VecDoub_I &yy, VecDoub_I &ssig,
.
.
.

Now, refer to nr3.h:

typedef const NRvector<Doub> VecDoub_I;



Since they are declared VecDoub_I, these parameters are const, and it is absolutely, positively, guaranteed that nothing in this function (the Fitsvd constructor) can modify values in the corresponding arguments back in the calling function (your main program or wherever you created the object).


In particular, the line you quoted from the Fitsvd constructor is calculating values in a locally-declared vector, b, which is one of the arguments that it feeds to the solve() member function of the locally declared SVD object.

I hate to repeat myself, but the Fitsvd constructor does not, and can not, modify any of the values of its vector arguments back in the calling function.


I didn’t use Fitsvd at all...

It's certainly possible to do "certain things" and feed your matrix to an SVD constructor so that calling SVD::solve() can work for you. Section 15.4.2 describes those "things" as it leads up to Fitsvd.

But...

Why would you want to do the work of Fitsvd, when it has already been done for you? I mean, I completely understand that it is a learning experience (and, maybe, a Good Idea) to try to reproduce the functionality of Fitsvd with your own code, but is that what you really want (need) to do for your project?

If you have specific questions about why your code fails, then, if you can post the code, maybe someone can spot (and give back) some helpful clues.


Regards,

Dave

Hi Dave

Thank you so much for your response. The reason that I didn't use Fitsvd because I already got everything in place; namely, a class that handles basis function, a class that does the regression, etc, so I just want the SVD algorithm to give me back the required vector x.

In particular, the line you quoted from the Fitsvd constructor is calculating values in a locally-declared vector, b, which is one of the arguments that it feeds to the solve() member function of the locally declared SVD object.


Just to confirm that I understand this correctly, is this vector b that will be entering a 'solve' function the usual dependent variable vector in the regression equation?


Anyway, I will implement the Fitsvd and see what my results look like and I may post my code if the result is still incorrect.

Cheers and Thank you again

PChira

...I already got everything in place...
Yes, you did say that. It just didn't register. I'm really sorry; sometimes I simply miss the most important points. Thanks for your patience.

Just to confirm that I understand this correctly...The input vector (or matrix) [b] to SVD::solve(b,x, thresh) is the right-hand side of the equation [A]*[x] = [b]; the output vector (or matrix) [x] is the "best" solution that SVD can get with the given matrix and right-hand side.

Reference: Section 2.6 in the book. Equations 2.6.5, 2.6.7


Regards,

Dave

Perhaps another angle on this may help. I was interested in trying my own network adjustment to a typical surveying problem. Since SVD can solve over-determined systems of equations, this is what you typically get when you run a level survey between a set of points and you collect more than the required number of shots to solve for all the elevations.

A simple example is a piece of land with 4 corner pins: a,b,c,d. You set up your instrument and site each corner, including the 2 diagonal directions. You end up with 6 shots and 4 unknowns. If you set up your A matrix where the columns represent the 4 unknowns a,b,c,d, and the 6 rows represent the 6 shots, then the elevations differences between the point shot from and the point shot to is the vector b. Each equation is coded with a -1 at the point shot from and a 1 at the point shot to.
A =
-1 1 0 0
0-1 1 0
0 0 -1 1
-1 0 1 0
0 -1 1 0

b= [dEab,dEbc,dEcd,dEda,dEac,dEbd]T

where dEab = elevation at b minus elevation at a.

Solve for x.

// read in m,n prior to this
VecDoub b(m);
VecDoub x(n);
MatDoub a(m,n,0.0);

// read in A and b
// Create an SVD object
SVD svd(a);
svd.solve(b,x,-1);

// print out x

I tried this with a project involving 21 unknown points and 50 equations and it worked very well. I got a solution vector where unknown 1 was offset; it seems my straight-forward setup required a shift to put [b] on track, but otherwise, all singular rows and columns were dealt with appropriately. After working with other solvers this is nothing short of a dream.

Yes, this is a nice example of a set of equations that is both overdetermined (more equations than unknowns) and underdetermined (arbitrary overall constant can be added to all the elevations)!

You might want to add a row like
1 0 0 0
to your matrix, and an entry
0
to the right hand side vector. Then, all elevations will be relative to the first point, which will be at elevation 0.