HP Forums

CAS
SVD([[1,2,3,4],[4,3,2,1],[-2,1,4,7]])-->[undef]

Home
SVD([[1,2,3,4],[4,3,2,1],[-2,1,4,7]])-->

{[[−0.816496580928,−0.218217890236,0.534522483825],[0.408248290464,−0.872871560944,0.267261241912],[0.408248290464,0.436435780472,0.801783725737]],[0,5.47722557505,10],{{undef,−0.836660026534,3.5527136788ᴇ−16},{undef,−0.478091443734,0.267261241912},{undef,−0.119522860933,0.534522483825},{undef,0.239045721867,0.801783725737}}}

Han's svd2([[1,2,3,4],[4,3,2,1],[-2,1,4,7]])-->Correct
{[[−0.534522483825,−0.218217890236,0.816496580928],[−0.267261241912,−0.872871560944,−0.408248290464],[−0.801783725737,0.436435780472,−0.408248290464]],[10,5.47722557505,9.35472980649ᴇ−15],[[−6.2172489379ᴇ−15,−0.267261241912,−0.534522483825,−0.801783725737],[−0.836660026534,−0.478091443734,−0.119522860933,0.239045721867],[−0.540390815439,0.653947124252,0.313278197813,−0.426834506626],[−8.93183440801ᴇ−2,0.521874658019,−0.775794283799,0.343237969859]]}

as in
http://www.wolframalpha.com/input/?i=svd...,4,7%7D%7D

The built-in SVD command does not work for under-ranked matrices, here you have rank 2 (instead of 3 for 3x4 matrices). There is a clear warning about that (0 as singular value...).

Thanks...parisse;
this discussion led to this code; pseudoinverse

ppinv(matrix)
example; ppinv([[1,2,3,4],[4,3,2,1],[-2,1,4,7]])

Code:

EXPORT ppinv(MM1)

//PSEUDOINVERSE ppinv(matrix)

BEGIN

local LL,L1,MM2,MM11,aa;

local AA:=0;

local BB:=0;

aa:=DIM(MM1); 

IF aa[1]>aa[2] THEN MM1:=transpose(MM1); AA:=1  END;

IF RANK(MM1)<MIN(aa) THEN LL:=QR(MM1);  BB:=1; END;

IF aa[1]==aa[2] AND (RANK(MM1)<MIN(aa)) THEN MM1:=(LL[1]*LL[2]*transpose(LL[2])*transpose(LL[1])); END;

IF aa[1]==aa[2] AND (RANK(MM1)==MIN(aa)) THEN

LL:=SVD(MM1);

L1:=size(LL[2]);

FOR A FROM 1 TO L1[1] DO

IF LL[2,A]<>0 THEN LL[2,A]:=1/LL[2,A]; END;

END;

END;

CASE

IF aa[1]==aa[2] THEN return (LL[3]*diag(LL[2])*transpose(LL[1])); END;  

IF aa[1]<>aa[2] THEN 

CASE

IF BB==0 THEN MM2:=MM1*transpose(MM1); MM11:=MM1;  END;

IF BB==1 THEN MM2:=(LL[1]*LL[2]*transpose(LL[2])*transpose(LL[1])); MM11:=MM1; END;

END;

MM2:=(transpose(MM11)*ppinv(MM2)); 

IF AA==1 THEN return transpose(MM2); ELSE return (MM2); END;

END;

END;

END;

A general SVD program has been implemented here:

http://www.hpmuseum.org/forum/thread-4976.html

and also provides a command for the pseudo inverse via pinv().

Wondering a bit more about the built in SVD function. I used it with a M1 = 9x4 matrix. I expected to get L1 = { U = 9x4, D = [4] vector, and V = 4x4 }. By definition U*D*V` should reproduce M1 but I get U = 9x9, so that multiplication is impossible. What are the extra 5 columns in U?
Thanks?

They come from a qr factorization, and you are right, I will remove them.
Let me explain a little bit: numeric algorithms is not my main speciality, it's symbolics. On many platforms where giac is ported, I have access to lapack/atlas/etc. where SVD is available. On the Prime, these libs are not available. But I didn't want to take a lot of time and effort to implement a state of the art SVD algorithm (it's not something that students need on a calculator at least here in France, it's much work: look at Han's code, and imagine the additional work to implement all that code efficiently in C++). Instead I wrote something quick (and a little dirty).
The idea is the following: let m be a matrix with #rows>=#columns (otherwise tranpose), if m=u*diag(d)*trn(q) then trn(m)*m=q*diag(d)^2*trn(q), therefore you get q and d by diagonalization of trn(m)*m. Then you get u, except if 0 is eigenvalue of d (m not full rank), in that case I complete u by doing a qr factorization.
It's bad for matrices that are not full-ranked or more generally having a large condition number. But for that kind of matrices, you can use Han's svd program or switch to a port of giac with lapack/atlas support.

The same crude but effective way mention my parisse is used to in the above code ppinv();
it is effective but was wonder how effective and guarded numerically is the qr implemented on prime.
For ill mannered matrix a well guarded decomposition can be used to solve SVD even with the crude method implemented on prime.
SVD(A_ill) into qr(A_ill) into SVD(R)=U*S*V'
hence
A=(Q*U)*S*V'
i used Han's SVD (Tridiagonal Method) but too much to port into any application; hence i resulted to the above code which is faster;

Same condition number since it will diagonalize R^t R=A^t A, and cond(R^t R)=cond(A)^2 (for the Euclidean norm). This means that if say cond(R)=100, you could expect to loose 2 digits with a state of the art SVD and 4 with the quick and dirty method.
Fortunately the situation is not always that bad. I tried with hilbert(6), despite a high condition number 1.5e7, the SVD is accurate.

Thanks for the explanation and I look forward to the future correction. I will use SVD2 by Han.

please BP, incorporate HAN-SVD2 into the CAS core
Thanks

toshk

parisse

toshk

Han

mark4flies

parisse

toshk

parisse

mark4flies

compsystems