Disagreement between R's eigen() and nimEigen() ?

[danielturek]

There appears to be disagreement between the numerical results of R's eigen(), and our nimEigen().

Is this a problem?

Reproducible example:

M <- 1:4
M1 <- diag(M^-0.5)
M2 <- diag(M^ 0.5)

Cmatrix <- t(array(c(0,1,0,0,  .5,0,.5,0,  0,.5,0,.5,   0,0,1,0), c(4,4)))

x <- M1 %*% Cmatrix %*% M2

eigen(x)

library(nimble)
nimEigen(x)

[NLMichaud]

Good question! I think this is occurring because here x is a non-symmetric matrix, and nimble's version of the eigen() function is not equipped to handle non-symmetric matrices. I recall this design choice being made because non-symmetric matrices can have complex eigenvalues, which nimble can't work with. This is noted in the help file for nimEigen(), but the function itself does not emit a warning or error message when a non-symmetric matrix is supplied.

Adding in a check for symmetry and an error message could be a good idea, since nimEigen() will return incorrect / useless results if a non-symmetric matrix is supplied.

[perrydv]

Sounds like that is in the class of "silent errors" where we are trying to improve our practices. If there is much computational burden to checking symmetry, it could be a user option to do so that defaults to TRUE, so that unless the user is willing to guarantee symmetry, it will be checked. What do you think?

[NLMichaud]

I think including a user option, defaulting to TRUE, is a good idea. I'm not sure of how expensive symmetry checking is, but I imagine for large matrices it could be noticeable, so giving the option to disable checking would be nice.

[danielturek]

So stepping back, is it the case then that we've decided that NIMBLE's use of eigen will be restricted to symmetric matrices? This concerns me, since there are certainly valid uses for an eigen decomposition of a non-symmetric matrix. In this case, it's necessary for finding the bounds of the gamma parameter of the car_proper distribution.

I've obviously missed some of the design discussions & decisions here, but it seems to make more sense to do something like:

• If a user guarantees symmetry, then we do our current processing.
• If they don't guarantee symmetry, then we check for it.
  • If it's symmetric, do our current processing
  • If it's non-symmetric, then do the correct processing for non-symmetric matrices. If this results in a complex eigenvalue, then we return NaN (for example) and issue a warning.

Again, I know I'm late to this discussion, but it does seem there are use cases for non-symmetric matrices.

[NLMichaud]

Good point @danielturek, I hadn't considered uses for non-symmetric eigendecompositons but there definitely are many, and it seems like implementing that functionality would be useful and straightforward. One question about the proposed logic flow:

• Should we check for symmetry if it's not guaranteed by the user, as you propose, or should we just do the more general (and slower) eigendecomposition that will work for both symmetric and non-symmetric matrices without symmetry checking?
  • On the one hand, checking for symmetry will allow faster, more accurate eigendecompositions if the matrix is determined to be symmetric.
  • On the other hand, checking for symmetry will impose a computational burden each time the eigendecomposition is taken for a non-symmetric matrix.

It's not clear to me which is the better idea here. I may have a slight preference for defaulting to the general eigendecomposition, without symmetry checking, as it seems like usually a user would know whether their matrix was symmetric or not before running an algorithm and thus could set symmetric = TRUE. Any thoughts or preferences?

[fritzo]

Fact Check: It is very cheap to check symmetry (O(n^2) to check symmetry vs O(n^omega) for eigendecomposition). Even when computing only a single eigenvector, there is only small overhead to check symmetry.

[danielturek]

Good points all around. So it could be wise to:

• Have a flag for symmetry (default = FALSE). If a user specifies eigen(x, symmetric = TRUE), then we do the faster eigen decomposition.
• When symmetric = FALSE (the default), we check for symmetry ourselves:
  • If it's symmetric, we do the faster eigen decomposition.
  • If it's not symmetric, we do the slower, correct eigen decomposition.

Come to think of it, I think this logic is basically the same that I proposed earlier.

Any alternative suggestions to this idea?

[NLMichaud]

@danielturek you're right, I think that's definitely the way to go. I'm working on a PR that will fix this issue, should be ready shortly!

[danielturek]

Fantastic @NLMichaud . I'll look forward to the update.