# Vectorized scale\_matrix\_exp\_multiply

**URL:** <https://discourse.mc-stan.org/t/vectorized-scale-matrix-exp-multiply/4802>\
**Category:** Developers\
**Tags:** features\
**Created:** [July 9, 2018, 4:52pm UTC](https://discourse.mc-stan.org/t/vectorized-scale-matrix-exp-multiply/4802 "2018-07-09T16:52:28Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Stijn\_de\_Waele](https://avatars.discourse-cdn.com/v4/letter/s/3e96dc/32.png) [@Stijn\_de\_Waele](https://discourse.mc-stan.org/u/Stijn_de_Waele)\
**Post date:** [July 9, 2018, 4:52pm UTC](https://discourse.mc-stan.org/t/vectorized-scale-matrix-exp-multiply/4802/1 "2018-07-09T16:52:28Z")

</div>

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:

```stan
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](https://github.com/stan-dev/stan/pull/2548), i.e. it would get the signature:

```stan
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.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 16, 2018, 1:19pm UTC](https://discourse.mc-stan.org/t/vectorized-scale-matrix-exp-multiply/4802/2 "2018-07-16T13:19:25Z")

</div>

@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:

```no-highlight
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

```no-highlight
matrix_exp(t * A) = matrix_exp(A)^t.

```

other than for integer powers of 2.

If so, then there’s an obvious savings.

> [@Stijn\_de\_Waele](#):
>
> for (t in ts)

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.

---

<div class="post-metadata">

**Author:** ![yizhang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/yizhang/32/15673_2.png) [@yizhang](https://discourse.mc-stan.org/u/yizhang)\
**Post date:** [July 16, 2018, 1:25pm UTC](https://discourse.mc-stan.org/t/vectorized-scale-matrix-exp-multiply/4802/3 "2018-07-16T13:25:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![yizhang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/yizhang/32/15673_2.png) [@yizhang](https://discourse.mc-stan.org/u/yizhang)\
**Post date:** [July 16, 2018, 1:37pm UTC](https://discourse.mc-stan.org/t/vectorized-scale-matrix-exp-multiply/4802/4 "2018-07-16T13:37:55Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Stijn\_de\_Waele](https://avatars.discourse-cdn.com/v4/letter/s/3e96dc/32.png) [@Stijn\_de\_Waele](https://discourse.mc-stan.org/u/Stijn_de_Waele)\
**Post date:** [July 17, 2018, 3:41am UTC](https://discourse.mc-stan.org/t/vectorized-scale-matrix-exp-multiply/4802/5 "2018-07-17T03:41:37Z")

</div>

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:

```stan
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](https://github.com/stan-dev/stan/issues/680) that computation of derivatives is a challenge here.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 17, 2018, 3:07pm UTC](https://discourse.mc-stan.org/t/vectorized-scale-matrix-exp-multiply/4802/6 "2018-07-17T15:07:00Z")

</div>

> [@yizhang](#):
>
> 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,

Cool!

> [@Stijn\_de\_Waele](#):
>
> Yes, the following holds for `n` integer:

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