Hi!
I have no clue how hard this is numerically, but would it make sense to enable a function in Stan which outputs the logarithm of the matrix exponential?
The most obvious problem is to me that this may not be applicable to all real square matrices.
The intention is to gain numerical stability such the usual tricks to perform calculations on the log scale throughout a Stan program can be applied (instead of working on the natural scale and then taking logs).
Sebastian
No idea about the algorithm or making it more stable.
Is it unstable now?
There’s a matrix log that inverts matrix exp:
Looks complicated. And only applies to invertible
matrices.
Looks complicated was also my impression; so I will do some research here.
I have been experiencing lots of numerical warnings whenever I fit PK models on the natural scale. This is why I have switched long ago all calculations of PK models into log-space. This was tedious, but it was well worth the effort I think. Now I am getting almost zero numerical warnings from my Stan runs. The same should hold true for the matrix exponential.
Just from my intuition it should be doable… I mean the log of any matrix will be hard, yes; but I want the log of the matrix exponential such that hopes are that this is better behaved. The motivation is the same as to why we want log_sum_exp stuff.
Just to clarify, are we talking about matrices which we know we can express in the form A = exp(B), but with B unknown?
To be specific what I have in mind is actually a log of a matrix; so this:
log( exp(A) * x0 )
with A being some matrix. Sorry that was probably misleading in the beginning. This is probably harder than it looks like.
Sebastian