SVD problem.?bug

01262018, 10:18 AM
Post: #1




SVD problem.?bug
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 

01262018, 04:53 PM
Post: #2




RE: SVD problem.?bug
The builtin SVD command does not work for underranked matrices, here you have rank 2 (instead of 3 for 3x4 matrices). There is a clear warning about that (0 as singular value...).


01292018, 08:43 PM
(This post was last modified: 01292018 08:51 PM by toshk.)
Post: #3




RE: SVD problem.?bug
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:


01302018, 06:17 PM
Post: #4




RE: SVD problem.?bug
A general SVD program has been implemented here:
http://www.hpmuseum.org/forum/thread4976.html and also provides a command for the pseudo inverse via pinv(). Graph 3D  QPI  SolveSys 

02032018, 05:47 PM
Post: #5




RE: SVD problem.?bug
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? 

02032018, 07:05 PM
Post: #6




RE: SVD problem.?bug
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 fullranked 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. 

02032018, 09:20 PM
(This post was last modified: 02032018 09:25 PM by toshk.)
Post: #7




RE: SVD problem.?bug
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; 

02042018, 06:49 AM
(This post was last modified: 02042018 07:11 AM by parisse.)
Post: #8




RE: SVD problem.?bug
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. 

02042018, 01:36 PM
Post: #9




RE: SVD problem.?bug
Thanks for the explanation and I look forward to the future correction. I will use SVD2 by Han.


02042018, 03:38 PM
Post: #10




RE: SVD problem.?bug
please BP, incorporate HANSVD2 into the CAS core
Thanks http://www.hpcalc.org 20 years online 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)