SVD problem.?bug - Printable Version +- 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:
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 |