Hi!
I just read in the depth of the matrix exponential pull that Ben has a Schur decomposition on a branch… this would be extremely useful as we often need the repeated evaluation of the matrix exponential (solutions at different times). Having the individual pieces to do the matrix exponential from the Schur decomposition would be very beneficial for speeding up the repeated evaluation.
So… what’s missing? Can this be exposed to the language? Is help needed?
Sebastian
What do you mean by repeated evaluation of the matrix exponential? I wrapped the real Schur decomposition from Eigen, but I’m not sure how that would avoid repeated evaluation.
Hi!
The Schur decomposition gives me Q and T for A:
A = Q T Q’
when I want exp(A * t) as needed for linear ODE systems for which I want the solution at different times t, then given the Schur decomposition leads me to:
exp(A t) = Q exp(T t) Q’
So I am left with the problem of solving repeatedly the Matrix exponential of a upper tridiagonal matrix. That is a far simpler problem and the Q and T matrices never changes when evaluating different time-points t… even better would be a Jordan decomposition.
I am not 100% sure if this solves problems of arbitrary size, but one should be able to work out formulas for 4x4 or even 5x5 cases which take such a tridiagonal matrix and give the exponential of it.
Sebastian
Doing a matrix exponential for T * t for an upper Hessenberg matrix T is not much faster than doing a matrix exponential for A = Q T Q’. T and A have the same dimensions and the Pade approximation algorithm for the matrix exponential that Charles is implementing / wrapping is not specialized for structured matrices. The upper Hessenberg structure might make things more numerically stable because scaling and squaring does not alter the zeros. For a 4x4 or smaller matrix, we could probably do some specializations.
Hi!
Even if we can only work out the smaller cases 4x4 or 5x5 in analytical form, the speedup should be substantial due to the repeated evaluation.
These specializations can at the moment only be done in the Stan language itself as only then one can at least try to take advantage of the tridiagonal structure - which one should be able to attack in a recursive manner (or I am happy to throw Mathematica onto this problem).
Having 4x4 or even 5x5 coded in a super-fast way would be very helpful, I think.
Sebastian
What do you mean that it can only be done in the Stan language itself?
We don’t have a matrix type representing tridiagonal matrices, but if you operate within the Stan language, then I think one can take at least some shortcuts. For example, one could throw the upper triangle into a appropiatley sized vector which contains only the non-zeros or/and construct loops which only access those non-zeros, etc. One could possibly do this in C++ as well, but I am not sure if we want those hacks inside stan-math, I don’t think so.
We would need that tridiagonal matrix type in Stan to get the full benefit (I am not saying we should implement this - implementing new types sounds like a scary amount of crazy work).
Sebastian
The T matrix in the real Schur decomposition is only going to be tridiagonal if the A matrix is antisymmetric. Is that the case in your models?