Vectorized scale_matrix_exp_multiply

Developers
1

I would like to compute the solution of an ODE at several (irregularly spaced) times ts. This computation can currently be done in a loop:

for (t in ts)
  y[t] = matrix_exp(t * A) * y0;

My question is: Would it be possible to reduce the computational load by vectorizing this computation? I assume @yizhang is probably the most knowledgeable on this topic.

This could naturally be included in the scale_matrix_exp_multiply function that recently was exposed to Stan, i.e. it would get the signature:

scale_matrix_exp_multiply(vector t, matrix A, matrix B)

If it is expected that this can reduce computationally load, I can submit a feature request for this.

2

@charlesm93 wrote the matrix_exp function and @yizhang wrote the extensions.

I don’t know if there’s anything to be gained by a repeated A argument for varying t in matrix_exp(t * A), but one thing you can do for more efficiency is:

for (t in ts)
  y[t] = matrix_exp(t * A);
y *= y0;

We have matrix_exp(A + A) = matrix_exp(A) * matrix_exp(A) because A commutes with itself, but I don’t know if that means

matrix_exp(t * A) = matrix_exp(A)^t.

other than for integer powers of 2.

If so, then there’s an obvious savings.

Wow, you’re already using our foreach syntax? It was just released in CmdStan a couple days ago!

Otherwise, for multiple time points, the ODE solvers might be competitive. And they have the advantage of allowing you to specify tolerances.

With multiple time points, the ODE solver may also be efficient enough.

3

For each t, matrix_exp_multiply would go through different iteration paths for exp(tA)*b so it’s not obvious how we can vectorize the function for effects more than aesthetic.

4

In order to use matrix_exp_multiply to repeatedly integrate ODE through a set of steps efficiently, one needs to incorporate previous steps’ solution to move forward, something like

B_next = scale_matrix_exp_multiply(h, A, B_old);

This is based on the nature of ODE solution instead of simple vectorization, though the final UI would be what you proposed. I plan to work on this after getting PDE interface in.

5

Thanks for your answers! I will use your advice to incorporate the previous steps’ solutions. It seems that there is no obvious acceleration to be had through vectorization because of the iteration paths.

@Bob_Carpenter: Yes, the following holds for n integer:

matrix_exp(A*n*t0) = matrix_exp(A*t0)^n

This can be used to accelerate the matrix exponential example in the documentation. The only possible caveat I see is accumulating errors, which perhaps are less when using the current implementation.

Anyway, for the problem I am considering at the moment, I am dealing with irregularly spaced times, so I can’t exploit this rule.

When using the computation for diagonalizable matrices, there is an obvious benefit to vectorization because the eigenvalue decomposition has to be computed only once. However, I did read in a previous thread that computation of derivatives is a challenge here.

6

Cool!

Thanks. It’s a little above my matrix algebra skills to work this out myself and I couldn’t find the answer anywhere after a quick look other than that it was obvious for powers of two since A * A = A * A.