[Rd] Catching errors from solve() with near-singular matrices

From: David Sterratt <david.c.sterratt_at_ed.ac.uk>
Date: Tue, 11 Dec 2012 15:13:35 +0000


Dear all,

The background is that I'm trying to fix this bug in the geometry package:
https://r-forge.r-project.org/tracker/index.php?func=detail&aid=1993&group_id=1149&atid=4552

Boiled down, the problem is that there exists at least one matrix X for which det(X) != 0 and for which solve(X) fails giving the error "system is computationally singular: reciprocal condition number = ..." (see appended code & attached file). I don't want my function that calls solve(X) to return an error.

I can think of two strategies for dealing with this problem:

Strategy 1: Some code like this:

   if (det(X) < epsilon) {

      warning("Near singular matrix") 
      return(NULL) 

   }

   return(solve(X))

The problem is then to find what epsilon should be.

Strategy 2: Catch the error thrown by solve(X) like this:

        f <- function(X) {
          invX <- tryCatch(solve(X), error=function(e) {
            warning(e)
            error.flag <<- TRUE})
          if (error.flag) return(NULL)
          return(invX)
        }
        
This works OK if called without a surrounding try() 
        
        ret  <- f(matrix(0, 2, 2))      ## Gives warning
        

However, if I encase the call to f() in a try statement, I get an error:

        ret1 <- try(f(matrix(0, 2, 2))) ## Gives error "Lapack routine dgesv: system is exactly singular"

This is undesirable.

Advice on how to solve the problem with either strategy would be much appreciated - as indeed would be a completely different solution.

All the best,

David.

Code to throw an error in solve():

> X = as.matrix(read.csv("X.csv"))
> det(X)

[1] 2.32721e-21
> solve(X)

Error in solve.default(X) :
  system is computationally singular: reciprocal condition number = 1.79977e-16

-- 
David C Sterratt, Research Fellow http://homepages.inf.ed.ac.uk/sterratt
Institute for Adaptive and Neural Computation     tel: +44 131 651 1739
School of Informatics, University of Edinburgh    fax: +44 131 650 6899
Informatics Forum, 10 Crichton Street, Edinburgh EH8 9AB, Scotland  
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
NEW BOOK: Principles of Computational Modelling in Neuroscience
Sterratt, Graham, Gillies & Willshaw (CUP, 2011).
http://www.compneuroprinciples.org  

The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.


______________________________________________ R-devel_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel

Received on Tue 11 Dec 2012 - 15:18:50 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

All messages

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 12 Dec 2012 - 09:22:48 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive