[Rd] Catching errors from solve() with near-singular matrices
Ravi Varadhan
ravi.varadhan at jhu.edu
Tue Dec 11 16:46:09 CET 2012
I am not sure that this query is appropriate for r-devel, it seems to be more appropriate for r-help.
In any case, you might want to try MASS::ginv instead of solve(), if you expect ill-conditioning. Here is one possible solution:
f <- function(X) {
invX <- tryCatch(ginv(X, tol=.Machine$double.eps), error=function(e) {
warning(e)
error.flag <- TRUE}) # you should avoid the global assignment `<<-'
if (error.flag) return(NULL)
return(invX)
}
n <- 14
A <- matrix(NA, n, n)
for (i in 1:n) for (j in 1:n) A[i,j] <- 1/(i+j-1)
ret <- f(A)
ret1 <- try(f(A))
Hope this helps,
Ravi
-----Original Message-----
From: r-devel-bounces at r-project.org [mailto:r-devel-bounces at r-project.org] On Behalf Of David Sterratt
Sent: Tuesday, December 11, 2012 10:14 AM
To: r-devel at r-project.org
Subject: [Rd] Catching errors from solve() with near-singular matrices
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.
More information about the R-devel
mailing list