On Mar 16, 2011, at 1:32 PM, Berend Hasselman wrote:


Peter Langfelder wrote:

On Wed, Mar 16, 2011 at 8:28 AM, Feng Li <m...@feng.li> wrote:
Dear R,

If I have remembered correctly, a square matrix is singular if and only
if
its determinant is zero. I am a bit confused by the following code error.
Can someone give me a hint?

a <- matrix(c(1e20,1e2,1e3,1e3),2)
det(a)
[1] 1e+23
solve(a)
Error in solve.default(a) :
system is computationally singular: reciprocal condition number = 1e-17


You are right, a matrix is mathematically singular iff its determinant
is zero. However, this condition is useless in practice since in
practice one cares about the matrix being "computationally" singular,
i.e. so close to singular that it cannot be inverted using the
standard precision of real numbers. And that's what your matrix is
(and the error message you got says so).

You can write your matrix as

a = 1e20 * matrix (c(1, 1e-18, 1e-17, 1e-17), 2, 2)

Compared to the first element, all of the other elements are nearly
zero, so the matrix is numerically nearly singular even though the
determinant is 1e23. A better measure of how numerically unstable the
inversion of a matrix is is the condition number which IIRC is
something like the largest eigenvalue divided by the smallest
eigenvalue.


svd(a) indicates the problem.

largest singular value / smallest singular value=1e17 (condition number)
--> reciprocal condition number=1e-17
and the standard solve can't handle that.

Actually it can if you relax the default tolerance settings:

> solve(a, tol=1e-21)
      [,1]   [,2]
[1,]  1e-20 -1e-20
[2,] -1e-21  1e-03
> a%*%solve(a, tol=1e-21)
    [,1] [,2]
[1,]    1    0
[2,]    0    1

(I posted a similar reply to the OP soon after his complaint hit the list, but the system seems to have eaten it.)

--

David Winsemius, MD
West Hartford, CT

(pivoted) QR decomposition does help. And so does SVD.

Berend

--
View this message in context: 
http://r.789695.n4.nabble.com/Singularity-problem-tp3382093p3382465.html
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

David Winsemius, MD
West Hartford, CT

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to