Vectorized scale_matrix_exp_multiply



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.