[R] Singularity problem

[William Dunlap]

You can get fairly bad results from solve() and
other linear algebra routines without warning.

E.g., the following function makes a 2 by 2 matrix which
would have determinate 1.0 if we had infinite precision
arithmetic and actually will produce one of determinant
1 on a computer if you use the "right" formula, even for
fairly large n:
  f <- function (n) {
    matrix(c(n, n - 1, n + 1, n), 2, 2)
  }
Here is the school formula for the determinant:
  det2by2 <- function (x) {
    stopifnot(ncol(x)==2, nrow(x)==2)
    x[1, 1] * x[2, 2] - x[1, 2] * x[2, 1]
  }
Try it for n=10^7:
  > x <- f(10^7)
  > det2by2(x) == 1
  [1] TRUE
However, the general linear algebra routines give fairly
poor results without any warnings:
  > solve(x, x) # expect identity, wrong in 3rd decimal place
              [,1] [,2]
  [1,] 0.994194930    0
  [2,] 0.005805069    1
  > det(x) # c. 0.5% error again
  [1] 1.005837
Usually rcond() or eigen() or svd() will give you
a hint that the matrix is ill conditioned
  > rcond(x)
  [1] 2.514593e-15
but knowing you are dealing with odd matrices and
using appropriate procedures can help also.  (E.g.,
you could create a class of 2x2 integral matrices
of determinant 1 and use the school formulas for
dealing with them.)
  
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com  

> -----Original Message-----
> From: r-help-bounces at r-project.org 
> [mailto:r-help-bounces at r-project.org] On Behalf Of David Winsemius
> Sent: Wednesday, March 16, 2011 1:12 PM
> To: Berend Hasselman
> Cc: r-help at r-project.org
> Subject: Re: [R] Singularity problem
> 
> 
> 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 at 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-tp3382093p33
82465.html
> > Sent from the R help mailing list archive at Nabble.com.
> >
> > ______________________________________________
> > R-help at 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 at 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.
>