Handling complex eigenvalues/-vectors in model calculations

[Chloe Nater]

Hello dev-team and fellow nimble-users,

I've recently been working with a Bayesian population model involving the "classical" approach of calculating growth rates and stable age distributions via eigen analysis of a projection matrix. 

It seemed like a perfect scenario to work with nimbleFunctions, and I successfully wrote a function for building the projection matrix ('make.sealPM' in attached code).

The function works as it should, but nimble spams the following warning messages whenever it is used: 

"Run-time warning: matrix used in call to nimEigen() has a complex eigenvalue."

"Run-time warning: matrix matrix used in call to nimEigen() has a complex valued 

eigenvector."

Is there a way to turn these off?

 

When written directly into the MCMC model, the eigen cacluations worked fine, presumably because nimble converts the complex numbers into real ones at some point. 
Ideally though, I'd write the eigen calculations as nimbleFunctions as well (to give a bit more flexibility and account for the possibility that the dominant eigenvalue may not always have the same index). 
Unfortunately, my attempt to do that did not go as well, again because of the complex numbers. So far, it looks to me like there is no nimble equivalent of R's "Re()" function (which extracts the real part of a complex number), but maybe I have misse something?

I am attaching an example to re-create/illustrate the issue. Any tips/suggestions welcome. 

Thank you!

[Chris Paciorek]

Hi Chloe,

On quick glance, I see that your A matrix has non-negligible imaginary parts for the eigenvalues so this is not just a floating point precision issue where some rounding error causes negligible but non-zero imaginary parts. Are you sure it's ok from a mathematical perspective to just throw away the imaginary parts? I haven't thought about matrix population models in a while, so that may be a naive question.

But if you did really want to do it, you could use nimbleRcall embedded within the nimbleFunction (or possibly directly within your model code) to call out to an R function (that you would write) that does the eigendecomposition (using R's eigen()) and the use of Re() on top of that.  Happy to give more details if that is a direction you want to go.

-chris

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/1c30409c-4ea8-4fc9-9ff3-dee43b393ecfn%40googlegroups.com.

[Daniel Eacker]

Hi Chloe,

I was curious if you've seen any of the functions available in the 'popbio' package in R (sorry I just saw this used in your code)? I've used some basic matrix algebra to produce asymptotic growth rates using the MCMC sampling in an integrated population model, but I think it would be much easier to just process the samples outside of nimble if you're seeking to derive these types of quantities.

I realize it's not as satisfying as including everything in the model, but I don't see anything wrong with a little post-processing to derive these quantities and I've seen it published.

Thanks,

Dan

[Chloe Nater]

Hi Chris, 

thanks a lot for getting back to me so quickly. 

As I have been told, positive definite orthogonal matrices (such as population projection matrices) should always have at least one non-complex dominant eigenvalue. I do believe that that is always the one we are looking to extract for the sake of the calculations here. So it's a complex number, but the imaginary part is actually 0 and nothing should be lost by dropping it. 

On that note, I am afraid I do really want to do it. nimbleRCall is a new one for me, and I'd be very happy to get some more details for using it :-) 

Cheers, 

Chloé

[Chloe Nater]

Hi Dan!

Thanks a lot for the suggestion. It is indeed the 'popbio' package functions that I am hoping to use within my nimble model. 

I totally agree that if this was just about calculating some asymptotic growth rates "post-hoc", then it would make much more sense to just derive those from the posterior samples. 

However, in this particular case, I am working with a very data-deficient integrated model, and am trying to use a stable age distribution calculated via eigen analysis of a matrix with estimated elements as the starting point of a stochastic population model. The choice of Bayesian framework here is mostly motivated by the fact that I have quite a few parameters that will need "expert opinion" as informative priors...

(That being said, it has crossed my mind to just build a large simulation outside nimble instead...)

[Perry de Valpine]

Hi Chloe,

That makes sense, if long-term growth rate is all you need.

nimbleRcall is a way to pass control back to an R function, even from compiled nimble code.  It does impose the overhead of using the R evaluator, much slower than compiled C++, but possibly not a problem for a quick step.  If speed becomes an issue, you could also do it via a nimbleExternalCall, which can call separately compiled code, if you get set up to use an eigen decomposition outside of nimble.

Have you made a start with the example in help(nimbleRcall)?  Since you're a seasoned nimble user that might be enough to get you under sail.  Please give another holler if that doesn't do it.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/1b3acab2-a729-4e78-a9be-c1df7e41f408n%40googlegroups.com.

[Chloe Nater]

Hi Perry, 
thanks for the details. I will look into nimbleRcall and see how far I get. Will report back on how it went later!

Chloé