merci beaucoup Jean Marie.
---- Original message ----
>Date: Wed, 21 Dec 2005 13:21:14 -0600
>From: [email protected] (Jean Marie Linhart, StataCorp LP)
>Subject: Re: st: bug in matrix svd
>To: [email protected]
>
>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/
*
* 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/
