[Rd] error message in chol() (PR#1061)

[ripley@stats.ox.ac.uk]

On Sun, 19 Aug 2001 jerome@stat.ubc.ca wrote:

> Full_Name: Jerome Asselin
> Version: 1.3.0
> OS: Windows 98
> Submission from: (NULL) (24.77.112.193)
>
>
>
> I am having accuracy problems involving the computation of inverse of
> nonnegative definite matrices with solve(). I also have to compute the
> Choleski decomposition of matrices. My numerical problems involving solve()
> made me find a bug in the chol() function. Here is an example.
>
> #Please, load the (8x8) matrix "mat" given below.
> #The matrix "mat" is indeed invertible as
> solve(mat)
> #exists. However, the command
> chol(mat)
> #claims that "mat" is singular. The correct error message should be
> #that "mat" is not symmetric positive definite.

Yes.  Message changed.  Actually, it only uses the upper triangle,
so symmetry is assumed, as the help page says....

> #In fact, the matrix "mat" should be symmetric, but if
> #we adjust for rounding errors by
> chol((mat+t(mat))/2)
> #we have the same error message in R, whereas S-Plus works (with
> #a warning message saying that the rank is not full).

Well, the symmetrized matrix is not even close to positive definite, so
what does `works' mean?

> eigen((mat+t(mat))/2)
$values:
[1]  2.3768888  1.8976031  1.4630238  0.9286787  0.9159388  0.5205478  0.5018393
[8] -0.1343581
 ...

L <- chol((mat+t(mat))/2)
t(L) %*% L - (mat+t(mat))/2

has one element of about 0.27.  That's just undocumented in S-PLUS.

Generally, if you want good numerical control you should be using eigen
and not chol, and if you know your matrix is symmetric, you can do better
than solve().  Have you considered the Matrix package?


> mat <- structure(c(0.587765115222886, 0.167899163888134, -1.17517413388306e-05,
>   2.84014489816759e-05, -0.00314742260306333, 0.00957598738126672,
>  -0.000181283657365524, 0.000268404362197193, 0.167899163888134,
>  0.830166933867836, 2.84014489816759e-05, -6.37164529062413e-05,
> -0.00221693916036902, 0.00290212917834065, 9.0762242654496e-05,
>  -0.000148298146106774, -1.17517413388306e-05, 2.84014489816759e-05,
>  0.592681746779869, 0.155708697468761, -0.000181283657365506,
>  0.000268404362197185, -0.00203952216264033, 0.00850865103383899,
>  2.84014489816759e-05, -6.37164529062413e-05, 0.155708697468761,
>  0.856464279975445, 9.0762242654496e-05, -0.000148298146106784,
>  -0.00296804648853004, 0.00336588501549565, -0.00314742260306333,
>   -0.00221693916036901, -0.000181283657365512, 9.07622426544948e-05,
>   1.72579173368880, 0.238005942444999, 0.120912001671674, -0.141897803793986,
>   0.00957598738126672, 0.00290212917834065, 0.000268404362197191,
>  -0.000148298146106782, 0.238005942444999, 1.56986591648628,
> -0.141897803793988,  0.213872027136185, -0.000181283657365518,
> 9.07622426544981e-05,  -0.00203952216264033, -0.00296804648853004,
>  0.120912001671674,  -0.141897803793987, 1.53425760353976, 1.16174845591861,
>  0.000268404362197194,  -0.000148298146106774, 0.00850865103383899,
>  0.00336588501549565,  -0.141897803793985, 0.213872027136185,
> 1.16174845591861, 0.773168806527852 ), .Dim = c(8, 8))
>
>
>
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> 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
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>

-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._