inverting singular matrixes with non-zero determinants

[josherrickson]

solve() will sometimes produce an undesirable result when the matrix is singular but the determinant is extremely close but not equal to 0.

> a <- cov(model.matrix(~ as.factor(gear) + 0, data = mtcars))
> a[1] <- a[1] + 1e-10
> det(a)
[1] 2.926638e-12
> det(a) == 0
[1] FALSE
> all.equal(det(a), 0)
[1] TRUE
> solve(a)
                 as.factor(gear)3 as.factor(gear)4 as.factor(gear)5
as.factor(gear)3            1e+10            1e+10            1e+10
as.factor(gear)4            1e+10            1e+10            1e+10
as.factor(gear)5            1e+10            1e+10            1e+10

Our generalized inverse code produces a more sensible inverse:

> s <- svd(a)
> nz <- (s$d > sqrt(.Machine$double.eps) * s$d[1])
> s$v[, nz] %*% (t(s$u[, nz])/s$d[nz])
           [,1]       [,2]      [,3]
[1,]  1.8944444 -0.3444444 -1.550000
[2,] -0.3444444  2.0666667 -1.722222
[3,] -1.5500000 -1.7222222  3.272222

This is related to #230 in the sense that the issue in 230 arose when an externally calculated distance matrix had extreme values due to this bug, however this is a separate concern as this happens internally.

Fix is to drop down to generalized inverse if either solve() errors OR the determinant is "close" to 0. I have a rough solution coded up but don't really have time right now to finish testing and polishing.

How to determine close to 0? .Machine$double.eps is 2e-16 on my macs, but that feels too small. MatchIt will be using 1e-8 for their threshold in their Mahalanobis calculation. I'm tempted to go even higher like 1e4 (adjusting the nz calculation as well) - is there a big downside to carrying out the generalized inverse in situations like this?

[benthestatistician]

1. Do we have a sense of the additional time cost of a det() call over and above the solve()? If it gets us up to the time cost of an svd(), then we might just simplify things and always do a generalized inverse.
2. I wonder if we should be changing these calcs without also bumping the major version number. Maybe not. Maybe the changes we made in resolving factor misunderstood to be a numeric in match_on.formula  #220 were already a bit more than should go into a minor version bump release, given their potential to change peoples' matching results. Next major version (0.10 or 1.0) changes #134.
3. I think an adjustment of this nature should be paired with the switch to inverting and generalized inverting correlation rather than covariance matrices, as suggested in this comment to 220.
4. If there does continue to be a solve()/svd() branching step, the tolerance I'd be inclined toward is sqrt(.Machine$double.eps). In other, only weakly related work I've tried .Machine$double.eps^0.25 and left thinking that was going to far. Not a firm opinion.

[josherrickson]

1. Speed shouldn't be a concern; we're inverting a covariance matrix with dimension pxp so at worst its O(f(p)) whereas any of our O(g(n)) operations will be much slower. I ran a quick test; det() is about twice as slow as solve(), but still 3.5x faster than our inverse.
2. I have no strong thoughts about version number.
3. Makes sense to me.

I do wonder how much this matters - I haven't gotten a chance to test it, but putting aside the issue of extreme values in the inverse causing numerical overflow issues, are the relative values similar enough that the relative distances would be same regardless of nonsensical solve() vs a more "proper" generalized inverse?

Edit: While there may be cases where this could be true, it's not always true - in the example in the opening comment solve(a) was producing a matrix with 9 identical values.

> sinv <- solve(a)
> all.equal(as.numeric(sinv), rep(sinv[1], 9), check.attributes = FALSE)
[1] TRUE

[benthestatistician]

Update: I'm on board with either the det()-based branching logic Josh proposed here or categorically using svd() for a generalized inverse. Would like to see inversion/gen inversion of correlation rather than covariance matrices and would like to use sqrt(.Machine$double.eps) for tolerance, as discussed above.