Domain Error from logpdf & occasional incorrect T matrix errors

[smith-garrett]

Error:

ERROR: DomainError with -4.5958295283709015e-27:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).

This happens occasionally. Most recently, I've noticed it using ReverseDiff.jl while fitting fpdistribs in Turing. I suspect that this might be due to proposals being sampled that cause the entries in the T matrix to be very close to zero. When that gets fed to pdf(), it might lead to numerical overflow.

I think something similar happens when I get feedback from the fpdistribution() fn. that the T matrix is incorrect, i.e., the columns don't sum to zero.

[smith-garrett]

The domain error hasn't happened yet with ForwardDiff.

It has happened with ReverseDiff even when I manually exclude very small scaling values ($\tau_i$). Might be an issue with ReverseDiff.

[smith-garrett]

Now the domain error has even happened with MH().

[smith-garrett]

This seems to be largely solved by using ExponentialAction.jl to compute pdfs, etc. It's really slow and doesn't seem to ever get to sampling with ReverseDiff.

For this class of model, it might work to go back to a more naive approach to matrix exponentials. The Schur decomposition works for any square matrix and might be less sensitive to weird eigenvalue spectra than an eigendecomposition. Could just sub in: $\exp(tA) = Q \exp(t * T) Q'$, where $Q$ and the upper-triangular $T$ are given by the Schur decomposition. Moler (2003, 19 dubious ways...) suggests the Schur decomposition could be promising, if somewhat less general, than the scaling and squaring method. Especially since the spectra of the matrices (the transient matrices) I intend to use generally have a single isolated maximal eigenvalue (so far that I've seen), eigenvalue clustering doesn't seem to be an issue.

[smith-garrett]

Another option would be to use an ODE solver to solve for the pdf(t) as $\sum \mathbf{A} \frac{dp_T}{dt}$. Would probably need to nest the solution to the ODE in another function that does the mulitplication by $A$.

[smith-garrett]

Yet another option: Switch to ExponentialUtilities b/c faster (speed and fewer allocations), and do sampling with MH, RWMH, emcee, ESS (elliptical slice sampling, seems to be more efficient than MH).

[smith-garrett]

Closing this until it comes up again....