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. #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). 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Mon, 20 Aug 2001, Prof Brian Ripley wrote:> 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....Ok. Thank you.> > #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?Thanks for the tips about eigen(). My attempt to reduce the use of solve() and to use eigen to resolve rank deficiencies appears to be successful. Jerome Asselin> > 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._