SVD Only With Matrix of Full Rank?

[mark4flies]

I am testing some code for multiple correspondence analysis, which relies on the SVD for dimension reduction. The example is 5x4 (5 observations on 4 responses) that is essentially two dimensions (last two singular values are essentially 0):

M1 = [ [8 2 2 -6], [5 0 3 -4], [-2 -3 3 1], [2 3 -3 -1], [4 6 -6 -2] ]

The result of SVD( M1 ) is this error (Textbook mode):

"0 as singular value, not implemented Error: Bad Argument Value"

(1) So the SVD() function does not allow less than full rank matrices?

Also, putting M1 on the stack and evaluating SVD results in this error (RPN Mode):

"SVD(NULL) Error: Bad Argument Type"

(2) How is SVD supposed to work in RPN mode then?

Thanks for help. I searched documentation (HP Prime and CAS) and other threads mentioning SVD first but did not see these questions yet.

[Han]

(04-03-2015 11:20 AM)mark4flies Wrote:  I am testing some code for multiple correspondence analysis, which relies on the SVD for dimension reduction. The example is 5x4 (5 observations on 4 responses) that is essentially two dimensions (last two singular values are essentially 0):

M1 = [ [8 2 2 -6], [5 0 3 -4], [-2 -3 3 1], [2 3 -3 -1], [4 6 -6 -2] ]

The result of SVD( M1 ) is this error (Textbook mode):

"0 as singular value, not implemented Error: Bad Argument Value"

(1) So the SVD() function does not allow less than full rank matrices?

Also, putting M1 on the stack and evaluating SVD results in this error (RPN Mode):

"SVD(NULL) Error: Bad Argument Type"

(2) How is SVD supposed to work in RPN mode then?

Thanks for help. I searched documentation (HP Prime and CAS) and other threads mentioning SVD first but did not see these questions yet.

For (2), you would type: COMMAND(n) where n is the number of arguments to take from the stack. So in your case, SVD(1) in RPN mode.

As for (1), this is due to an incomplete implementation in which the vectors of the unitary matrices are computed using the inverse of the singular values.

Once I have some free time to finish up my own implementation of the SVD, I'll post it here. It's basically an implementation of the "standard" algorithm (Householder reflections + Givens rotations)

[parisse]

I have modified the implementation in Xcas. It will handle 0 as singular value.

Code:


m1:=[ [8 ,2 ,2 ,-6], [5, 0, 3, -4], [-2 ,-3 ,3 ,1], [2 ,3 ,-3, -1], [4, 6 ,-6 ,-2] ];
u,s,q:=SVD(m1);
m:=u*[op(diag(s)*trn(q)),[0$4]]

Then m is equal to m1 up to numerical errors, but the quick and dirty algorithm I' using is unfortunately far from optimal for numerical stability if 0 is a singular value, especially here where it has multiplicity >1.

[mark4flies]

Thanks for the replies!

In the meantime, I tested this computation on the HP-50G, SAS, JMP, and R. All give the correct answer for U, D, and V for this case. I was hoping to port some HP-50G programs to the new HP Prime, but the newer ROM is not up to it.

[parisse]

Writing a good implementation of SVD takes time and it's not really part of a CAS. Inside Xcas, it is not required because LAPACK is called. On the calc, you will hopefully have my homemade quick and dirty implementation in a future firmware. If one can convince me that there is a real need for a state of the art implementation on a calc, I might improve this implementation later, otherwise it will stay like this, unless Han or someone else writes a user program for that.

[mark4flies]

The SVD is one of the most important results of matrix theory for practical applications, such as statistics. Almost all of the multivariate statistical techniques, such as principle components analysis and multiple correspondence analysis, use SVD at their core. This single function would support many real-world applications with a hand-held device.

[parisse]

For

Code:

m1:=[ [8 ,2 ,2 ,-6], [5, 0, 3, -4], [-2 ,-3 ,3 ,1], [2 ,3 ,-3, -1], [4, 6 ,-6 ,-2] ];u,s,q:=svd(evalf(m1),-1);m:=u*[op(diag(s)*trn(q)),[0$4]];

I get currently

Code:

s=[9.58068546525e-08,1.74628850627e-07,9.82257667369,14.1250482296]
m= [[8.0,2.0,2.0,-6.0],[5.0,-1.15463194561e-14,3.0,-4.0],[-1.99999994329,-2.99999993961,3.00000005965,1.00000011563],[1.99999995381,3.00000005855,-2.9999999624,-1.00000002953],[4.00000005145,6.00000000092,-5.99999998897,-1.99999992742]]

In other words you have at least 7 correct digits for a matrix of rank 2 (0 is singular value of multiplicity 2), of course it's more accurate for non degenerate matrices. I believe it is sufficient on a calculator for such kind of computations. I may be wrong, but if you want to convince me you need to provide better argumentation than generalities, because implementing a better conditionned algorithm is not trivial and is far away from my field of expertise, i.e. CAS.

[parisse]

I just discovered I can cheat on degenerated matrices, I can replace the small computed singular values by 0, then precision becomes much better. However it would have a drawback, really small svl would be replaced by 0. Is it a problem?

[mark4flies]

Thank you (Parisse) for looking into this matter. I am not an expert in numerical analysis. I have only started to learn HPP. I am just a user. I found this very recent article that might provide some help: it suggests a practical criterion for deciding if the singular value should be set to 0.

Again, thank you.

[mark4flies]

Link to article: http://blogs.sas.com/content/iml/2015/04...of-matrix/

[Han]

The svd2 program gives:

Code:


m1:=[ [8 ,2 ,2 ,-6], [5, 0, 3, -4], [-2 ,-3 ,3 ,1], [2 ,3 ,-3, -1], [4, 6 ,-6 ,-2] ];
svd2(m1);

{
  [ [−0.663425500474,  −0.457402727004, −0.592156525464,    3.5904651478e−13,  −8.40216785036e−13],
    [−0.364141989004,  −0.493987842429,  0.789542033952,    5.44011574989e−13,  4.52975068342e−13],
    [ 0.266854272705,  −0.301871594352, −6.57951694972e−2,  0.581376370892,     0.703800290354],
    [−0.266854272705,   0.301871594352,  6.57951694972e−2, −0.573305362697,     0.710390240949],
    [−0.53370854541,    0.603743188705,  0.131590338991,    0.577340866795,    −3.29497529792e−3] ],

  [ 14.1250482296, 9.82257667368, 6.3812064393e−12, 1.79164513489e−23],

  [ [−0.731350792814,   −0.433996954204,   0.168785260197,  0.498281243408],
    [−0.255197997351,    0.460050696925,  −0.797189774107,  0.296168537264],
    [−0.628778070935,    0.390531480179,   0.186669569957, −0.645970411205], 
    [−6.81038729371e−2, −0.668943318219,  −0.548775429161, −0.496711413046] ]
}

for U, S, and V respectively, where U*diag(S)*V = original matrix. The difference between U*diag(S)*V and m1 is

Code:

[ [−1.45803369378e−11, −2.13162820728e−12, −1.79767312147e−12,  4.64694949187e−11],
  [−5.17275111633e−12,  9.8119794999e−12,  −3.41060513165e−12,  2.09752215596e−11],
  [ 8.17124146124e−13,  4.41957581643e−12, −2.30215846386e−12, −3.47100126419e−12],
  [−8.17124146124e−13, −4.41957581643e−12,  2.30215846386e−12,  3.47100126419e−12],
  [−4.13535872212e−12, −4.32009983342e−12, −3.24007487507e−12,  9.85522774499e−12] ]

[mark4flies]

This result looks quite good compared to the result from other programs such as SAS IML or SAS JMP.

Is 'svd2' above the replacement for the built-in SVD function or a stand-alone program? I don't see 'svd2' in the catalog on my HP Prime with the latest ROM.

Thanks for all your help!

[Han]

(11-14-2015 03:37 PM)mark4flies Wrote:  This result looks quite good compared to the result from other programs such as SAS IML or SAS JMP.

Is 'svd2' above the replacement for the built-in SVD function or a stand-alone program? I don't see 'svd2' in the catalog on my HP Prime with the latest ROM.

Thanks for all your help!

It's a standalone program: http://www.hpmuseum.org/forum/thread-4976.html