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/
