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).

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.