[Rd] svd error (PR#631)

[Douglas Bates]

polzehl@wias-berlin.de writes:

> SVD-Error on
> 
> R 1.1.0
> Windows 98
> 
> 
> I get the following error applying  svd  on a positive definite matrix :
> 
>  > sk2
>                [,1]         [,2]          [,3]          [,4]          [,5]
> [1,]  1.0460139783  0.084356992 -2.810553e-04 -1.303237e-03  2.011809e-03
> [2,]  0.0843569924  1.154650879 -5.152560e-04 -2.389213e-03  3.688230e-03
> [3,] -0.0002810553 -0.000515256  1.000002e+00  7.960227e-06 -1.228821e-05
> [4,] -0.0013032374 -0.002389213  7.960227e-06  1.000037e+00 -5.697973e-05
> [5,]  0.0020118087  0.003688230 -1.228821e-05 -5.697973e-05  1.000088e+00

If you know that the matrix is positive definite then you are probably
better off taking the svd of the Cholesky decomposition of the matrix.

> svd(chol(sk2), nu = 0)
$d
[1] 1.095806 1.000000 1.000000 1.000000 1.000000

$v
             [,1]          [,2]         [,3]        [,4]        [,5]
[1,] -0.478709768  0.0011572932  0.013159984  0.87312442 -0.09119310
[2,] -0.877614104 -0.0003224026  0.006761341 -0.47348540  0.07456035
[3,]  0.002923979 -0.4202666121  0.844168101 -0.04500196 -0.32973090
[4,]  0.013558325  0.7925962305  0.519924862  0.03162204  0.31667897
[5,] -0.020929998  0.4417756749 -0.129766919 -0.10217730 -0.88154213

The singular values of sk2 will be the squares of the singular values
of its Cholesky decomposition.  The left and right singular vectors of
sk2 are identical and can be expressed as the V shown here.  Note that
since one singular value has multiplicity 4, the last four columns of
v are not unique.

>  > svd(sk2)
> Error in svd(sk2) : error  4  in dsvdc

This is related to the convergence in the Linpack subroutine dsvdc.
It is a reasonably good singular value routine but it will
occasionally be unable to converge.  Svd routines often have
difficulty with matrices with repeated singular values, like sk2.

The Lapack routine used in SVD from the Matrix library is more
reliable

> SVD(sk2)
$d
[1] 1.200791 1.000000 1.000000 1.000000 1.000000

$u
             [,1]          [,2]          [,3]          [,4]         [,5]
[1,] -0.478709768  5.814978e-19 -1.734663e-19 -2.917201e-20  0.877973210
[2,] -0.877614104 -2.109493e-02 -1.889409e-02 -3.985595e-03 -0.478513967
[3,]  0.002923979 -4.102994e-02  7.970942e-02 -9.959678e-01  0.001594282
[4,]  0.013558325  1.644232e-01 -9.825683e-01 -8.535898e-02  0.007392597
[5,] -0.020929998  9.853104e-01  1.668804e-01 -2.731482e-02 -0.011411959

$vt
           [,1]         [,2]         [,3]         [,4]        [,5]
[1,] -0.4787098 -0.877614104  0.002923979  0.013558325 -0.02093000
[2,]  0.0000000 -0.021094933 -0.041029944  0.164423228  0.98531038
[3,]  0.0000000 -0.018894086  0.079709424 -0.982568251  0.16688036
[4,]  0.0000000 -0.003985595 -0.995967801 -0.085358979 -0.02731482
[5,]  0.8779732 -0.478513967  0.001594282  0.007392597 -0.01141196

A new release of the Matrix library will be uploaded to CRAN shortly.  I
would advise waiting for it to be compiled for Windows before using
the SVD function extensively.  In looking at PR#630 by Torsten Hothorn
I found a subtle bug in SVD.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._