[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: SVD, EIG and CHOL of a matrix
From: |
Jaroslav Hajek |
Subject: |
Re: SVD, EIG and CHOL of a matrix |
Date: |
Wed, 20 Jan 2010 15:44:45 +0100 |
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
- SVD, EIG and CHOL of a matrix, Alberto Frigerio, 2010/01/19
- Re: SVD, EIG and CHOL of a matrix, Carlo de Falco, 2010/01/19
- Re: SVD, EIG and CHOL of a matrix, Jaroslav Hajek, 2010/01/20
- Re: SVD, EIG and CHOL of a matrix, Alberto Frigerio, 2010/01/20
- Re: SVD, EIG and CHOL of a matrix,
Jaroslav Hajek <=
- Re: SVD, EIG and CHOL of a matrix, Alberto Frigerio, 2010/01/22
- Re: SVD, EIG and CHOL of a matrix, Carlo de Falco, 2010/01/22
- Re: SVD, EIG and CHOL of a matrix, Jaroslav Hajek, 2010/01/22
- Re: SVD, EIG and CHOL of a matrix, Alberto Frigerio, 2010/01/25
- Re: SVD, EIG and CHOL of a matrix, Jaroslav Hajek, 2010/01/25
- Re: SVD, EIG and CHOL of a matrix, Alberto Frigerio, 2010/01/27