Re: SVD, EIG and CHOL of a matrix

[Jaroslav Hajek]

On Wed, Jan 20, 2010 at 1:56 PM, Alberto Frigerio
<address@hidden> wrote:
>
> I looked for a code to find, for example, the SVD decomposition of a matrix .
> I found it on the web and I tried to use it, but I got different results
> from the Octave function SVD ... as you would see if you try ro run the
> attached file, I got NaN results, which are not so good I believe .
>
> http://old.nabble.com/file/p27241250/svdmio.m svdmio.m
> http://old.nabble.com/file/p27241250/pythag.m pythag.m
>
>

Oh, I see. I didn't look at the code thoroughly, but one thing that
poked me into my eyes was
abs(sqrt(s)*sign(f)) --> you probably wanted abs(sqrt(s))*sign(f).

Octave uses the LAPACK library for doing the SVD, so you can start here:
http://www.netlib.org/lapack/double/dgesvd.f

If you are an Octave newbie, I suggest you first learn how to
vectorize your code and make use of high-level functions.
Good programs in Octave are typically quite different from loopy code
in Fortran or C.
It's not a bad idea to try reimplementing some of the basic
factorization, although you'd need to try very hard to outmatch the
built-in functions in accuracy and it's quite impossible to outperform
them in speed (as long as you code in m-files). I would, however,
start with something simpler than the SVD.

For instance, you can write this:

                scale = 0;
                for k=i:m
                        scale=scale+abs(A(k,i));
                endfor

better as

                 scale = sum (abs (A(i:m,i)));

or even better

                 scale = norm (A(i:m, i), 1);

Both ways will be much faster.

Similarly, this:

                        for k=i:m
                                A(k,i)=A(k,i)/scale;
                                s=s+A(k,i)*A(k,i);
                        endfor

could be written as

                        A(i:m) /= scale;
                        s = sumsq (A(i:m));

etc. It's a bit of an art, and also requires some knowledge about
available functions. But the results are normally faster, sometimes
more accurate, typically more readable and easier to debug, so it's
worth the effort.

happy Octaving & best regards

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz