[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

From |
[email protected] (Jean Marie Linhart, StataCorp LP) |

To |
[email protected] |

Subject |
Re: st: bug in matrix svd |

Date |
Wed, 21 Dec 2005 13:21:14 -0600 |

AbdelRahmen "U.A.QU.AP" <uaquap at gnet at tn> writes: > when I perform a singular value decomposition of a matrix G (say) I > expect that stata create 3 matrix U W V such that U W V =G where W > contains the singular value og G in DECREASING order (this is > important to perform later a test of rank of such marix), but Stata > put the smallest singular value in the top diagonal of W while the > rest of singular value are ordred as expected. > > How to fix such problem. > below the code: > mat svd U W V=G > mat list W > W[1,9] > y1 y1 y1 y1 y1 > y1 > r1 3.538e-18 .61151468 .61151468 .31692701 .31692701 > .16776301 > > y1 y1 y1 > r1 .16776301 .02892478 .02892478 This is not a bug -- -matrix svd- does not necessarily return the singular values in decreasing order. There's an example in [P] matrix svd showing this. Abdel Rahmen still wants SVDs in decreasing order, and this can be accomplished without too much sweat. When one sorts SVDs (the w matrix) one must also sort the columns in matrices U and V to reflect the sort order of w. Since I do not have Abdel Rahmen's matrix, I created an example of my own. This code should be pretty easy to modify for another matrix. Here is the business end of the code. I assume you already have defined matrix M, locals R and C containing its number of rows and number of columns respectively, and singular value decomposition matrices U,w and V so that M=U*w*V'. So w is the matrix of singular values, and it is not in decreasing order. This code segment will make newV, newU, and newW for which newU*newW*newV' = M is a singular value decomposition for M, and newW will be in decreasing order. ------Here's the code segment: mat newV = V mat newU = U mat newW = w /* I can sort the svds, but in order to keep the same decomposition, I have to also swap columns on the matricies U and W */ forv i = 1/`C' { local i1 = `i'+1 forv j = `i1'/`C' { if newW[1,`j'] > newW[1,`i'] { /* sort the svd values */ scalar temp = newW[1,`i'] mat newW[1,`i'] = newW[1,`j'] mat newW[1,`j'] = temp forv k= 1/`R' { /* exchange cols on newU to reflect the sorting */ scalar temp = newU[`k',`i'] mat newU[`k',`i'] = newU[`k',`j'] mat newU[`k',`j'] = temp } forv k = 1/`C' { /* exchange cols on newV to reflect the sorting */ scalar temp = newV[`k',`i'] mat newV[`k',`i'] = newV[`k',`j'] mat newV[`k',`j'] = temp } } } } ------End of the code segment. Here's the results I got when using it: . matlist M | c1 c2 c3 c4 -------------+-------------------------------------------- r1 | .0309208 .5281447 -.265614 -.0775027 r1 | .0538678 .3738923 -.3125677 .0430579 r2 | .011597 .1006803 .0999158 -.1563584 r3 | .1994718 -.0759514 -.02415 -.0380916 r4 | -.2715102 -.2585969 .1569542 .1472474 . local R = rowsof(M) . local C = colsof(M) . mat svd U w V = M . /* not in sorted order */ . mat list w w[1,4] c1 c2 c3 c4 r1 .84156286 .33406167 6.954e-17 .24098829 . mat newV = V . mat newU = U . mat newW = w . /* I can sort the svds, but in order to keep the same decomposition, > I have to also swap columns on the matricies U and W */ . forv i = 1/`C' { 2. local i1 = `i'+1 3. forv j = `i1'/`C' { 4. if newW[1,`j'] > newW[1,`i'] { /* sort the svd values */ 5. scalar temp = newW[1,`i'] 6. mat newW[1,`i'] = newW[1,`j'] 7. mat newW[1,`j'] = temp 8. forv k= 1/`R' { /* exchange cols on newU to > reflect the sorting */ 9. scalar temp = newU[`k',`i'] 10. mat newU[`k',`i'] = newU[`k',`j'] 11. mat newU[`k',`j'] = temp 12. } 13. forv k = 1/`C' { /* exchange cols on newV to > reflect the sorting */ 14. scalar temp = newV[`k',`i'] 15. mat newV[`k',`i'] = newV[`k',`j'] 16. mat newV[`k',`j'] = temp 17. } 18. } 19. } 20. } . /* Now I have the svds in sorted order */ . matlist newW | c1 c2 c3 c4 -------------+-------------------------------------------- r1 | .8415629 .3340617 .2409883 6.95e-17 . /* Now I check that I did the right thing -- that my new matrix is > the same as the old one ... */ . mat newM = newU*diag(newW)*newV' . /* assert that I got back what I started with, to a reasonable numeric > precision */ . assert mreldif(newM, M) < 1e-13 . mat Id = I(`C') . mat Ic = newV*newV' . assert mreldif(Ic, Id) < 1e-13 /* newV*newV' ~ Identity matrix */ . . /* in case you'd rather see this ... */ . matlist newM | c1 c2 c3 c4 -------------+-------------------------------------------- r1 | .0309208 .5281447 -.265614 -.0775027 r1 | .0538678 .3738923 -.3125677 .0430579 r2 | .011597 .1006803 .0999158 -.1563584 r3 | .1994718 -.0759514 -.02415 -.0380916 r4 | -.2715102 -.2585969 .1569542 .1472474 . /* It should be the same as the old one */ . matlist M | c1 c2 c3 c4 -------------+-------------------------------------------- r1 | .0309208 .5281447 -.265614 -.0775027 r1 | .0538678 .3738923 -.3125677 .0430579 r2 | .011597 .1006803 .0999158 -.1563584 r3 | .1994718 -.0759514 -.02415 -.0380916 r4 | -.2715102 -.2585969 .1569542 .1472474 . /* and it is! */ Hope that helps, --Jean Marie [email protected] * * For searches and help try: * http://www.stata.com/support/faqs/res/findit.html * http://www.stata.com/support/statalist/faq * http://www.ats.ucla.edu/stat/stata/

- Prev by Date:
**st: Testing model improvement for intreg, cluster?** - Next by Date:
**st: LaTeX production files for the Stata Journal have been updated** - Previous by thread:
**st: bug in matrix svd** - Next by thread:
**Re: st: bug in matrix svd** - Index(es):

© Copyright 1996–2024 StataCorp LLC | Terms of use | Privacy | Contact us | What's new | Site index |