+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: HP Prime (/forum-5.html)
+--- Thread: SVD problem.?bug (/thread-10024.html)
SVD problem.?bug - toshk - 01-26-2018 10:18 AM
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+%7B%7B1,2,3,4%7D,%7B4,3,2,1%7D,%7B-2,1,4,7%7D%7D
RE: SVD problem.?bug - parisse - 01-26-2018 04:53 PM
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...).
RE: SVD problem.?bug - toshk - 01-29-2018 08:43 PM
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;RE: SVD problem.?bug - Han - 01-30-2018 06:17 PM
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().
RE: SVD problem.?bug - mark4flies - 02-03-2018 05:47 PM
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?
RE: SVD problem.?bug - parisse - 02-03-2018 07:05 PM
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.
RE: SVD problem.?bug - toshk - 02-03-2018 09:20 PM
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;
RE: SVD problem.?bug - parisse - 02-04-2018 06:49 AM
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.
RE: SVD problem.?bug - mark4flies - 02-04-2018 01:36 PM
Thanks for the explanation and I look forward to the future correction. I will use SVD2 by Han.
RE: SVD problem.?bug - compsystems - 02-04-2018 03:38 PM
please BP, incorporate HAN-SVD2 into the CAS core
Thanks