# More efficient anisotropic calculation for Gaussian Process Kernel?

**URL:** <https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469>\
**Category:** Modeling\
**Tags:** gaussian-process, performance\
**Created:** [March 24, 2021, 4:17pm UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469 "2021-03-24T16:17:03Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![mathDR](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mathdr/32/1997_2.png) [@mathDR](https://discourse.mc-stan.org/u/mathDR)\
**Post date:** [March 24, 2021, 4:17pm UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469/1 "2021-03-24T16:17:03Z")

</div>

Hello! In Chapter 4 of Rasmussen and Williams (pg. 89) it mentions that  
“Anisotropic versions of these isotropic covariance functions can be created by setting  
`r^2(x, x') = (x − x')^T M (x − x')` for some positive semidefinite `M`.”

Assuming `M` is `DxD`, and `x` are `D-`dimensional vectors, and there are `N` of them, one could compute

```
for (i in 1:N-1) 
{
  Omega[i,i] = 0.0;
  for (j in (i+1):N) 
  {
    Omega[i,j] = quad_form_sym(M, (x[i]-x[j])');
    Omega[j,i] = Omega[i,j];
  }
}
Omega[N,N] = 0.0;

```

assuming `M` is symmetric positive definite.

My question is, how can this be vectorized? Is there a way of making a full matrix out of the `x` data and doing something more efficient than the double loop? Any suggestions are most welcome!

---

<div class="post-metadata">

**Author:** ![mathDR](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mathdr/32/1997_2.png) [@mathDR](https://discourse.mc-stan.org/u/mathDR)\
**Post date:** [March 25, 2021, 12:45am UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469/2 "2021-03-25T00:45:11Z")

</div>

Okay update: If I know that `M` has Cholesky factor `L` so `L L' = M`. Then would the following be faster:

```
matrix[D,D] L_transpose = L'; // I assume this costs flops so we should do it once and not in the loop?
for (i in 1:N-1) 
{
  Omega[i,i] = 0.0;
  for (j in (i+1):N) 
  {
    Omega[i,j] = dot_self(L_transpose * (x[i]-x[j])');
    Omega[j,i] = Omega[i,j];
  }
}
Omega[N,N] = 0.0;

```

That is, does Stan take advantage of the lower triangular structure in the matrix-vector multiply?

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [April 2, 2021, 9:17pm UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469/3 "2021-04-02T21:17:31Z")

</div>

Hi,  
sorry for not gettign to you earlier. I think this math stuff would be better handled by @stevebronder or @syclik, it goes a bit over my head. But I can answer a subquestion:

> [@mathDR](#):
>
> That is, does Stan take advantage of the lower triangular structure in the matrix-vector multiply?

AFAIK Stan doesn’t do any optimization like this “magically”. There are some specialized functions to handle triangular matrices like `mdivide_left_tri_low` but I don’t think there is anything for multiplication.

However, implementing this yourself should be easy. Stan is compiled to C++ and pays very little penalty for loops (unlike R). Obviously this may not help much, especially if this is not a bottleneck in your model.

Hope that helps at least a little and hope others will be able to provide a bit better advice.

---

<div class="post-metadata">

**Author:** ![mathDR](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mathdr/32/1997_2.png) [@mathDR](https://discourse.mc-stan.org/u/mathDR)\
**Post date:** [April 8, 2021, 8:13pm UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469/4 "2021-04-08T20:13:50Z")

</div>

So @martinmodrak you are suggesting something like:

```
for (i in 1:N-1) {
      Omega[i,i] = 0.0;
      for (j in (i+1):N) {
        real val = 0.;
        Omega[i,j] = 0.;
        for (k in 1:D) {
          for (l in k:D) {
            val += L[k,l]*(x[i,l]-x[j,l]);
          }
          Omega[i,j] += square(val);
        }
        Omega[j,i] = Omega[i,j];
      }
    }
    Omega[N,N] = 0.0;

```

(if I have my math right). Where `x` is `D` dimensional.

What impact (if any) does this have on the AD stack? Is the AD “smart enough” to know to aggregate the same values of `L[k,l]` for differentiation?

Thanks!

---

<div class="post-metadata">

**Author:** ![mathDR](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mathdr/32/1997_2.png) [@mathDR](https://discourse.mc-stan.org/u/mathDR)\
**Post date:** [April 10, 2021, 12:31am UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469/5 "2021-04-10T00:31:17Z")

</div>

Update: doing some timings (locally) this recent loop method is very slow.

Of the things I tried, in order of speed (from slowest to fastest–or least slow in this case):

`dot_self`, then `quad_form_sym` then `quad_form` version is fastest, but still very, very slow.

In python, I would precompute all of the `x[i]-x[j]` into a matrix, call it `X`.

Then `X` is `NxNxD` and the above computation (in python/numpy) would be:  
`Omega = (np.dot(X,M) * X).sum(axis=2).`

So question: is there a way of doing matrix multiplication with a tensor in stan?

I would like to compute `X` in the transformed parameters block as dimension `NxNxD` and do:

```
some_type Omega_pre_sum[N,N,D] = (X * M) .* X
for (i in 1:N) {
  for (j in 1:N) {
    Omega[i,j] = sum(Omega_pre_sum[i,j,:])
  }
}

```

is this possible?

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [April 13, 2021, 5:53pm UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469/6 "2021-04-13T17:53:27Z")

</div>

I admit, this is getting out of my expertise. Asking @stevebronder or @rok_cesnovar if they have time to discuss this.

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [April 20, 2021, 2:20pm UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469/7 "2021-04-20T14:20:14Z")

</div>

Bumping this up to hopefully get you more attention. Maybe @bgoodri has an idea?

I also think that partial indexing via `sum(Omega_pre_sum[i,j,])` or `sum(Omega_pre_sum[i,j])` should work out of the box.

---

<div class="post-metadata">

**Author:** ![mathDR](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mathdr/32/1997_2.png) [@mathDR](https://discourse.mc-stan.org/u/mathDR)\
**Post date:** [April 20, 2021, 2:37pm UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469/8 "2021-04-20T14:37:38Z")

</div>

Thanks for bumping this! Just an update:

Best so far (in terms of time): (assume `M = L*L'` in `DxD`)

```
for (i in 1:N-1) {
      Omega[i,i] = 0.0;
      for (j in (i+1):N) {
        Omega[i,j] = sum((X[i,j]*M) .* X[i,j]);
        Omega[j,i] = Omega[i,j];
      }
    }
    Omega[N,N] = 0.0;

```

where `X` is `row_vector[D] X[N,N]`

I will say, with `D = 64` and `N = 1000` this is extremely slow. After profiling it, the forward pass is about 1/3 the time and the autodiff takes 2/3.

I would love to “vectorize” this more and take advantage of matrix derivatives if there is a way to do so.

Any ideas are MOST welcome!

Dan

---

<div class="post-metadata">

**Author:** ![spinkney](https://avatars.discourse-cdn.com/v4/letter/s/dec6dc/32.png) [@spinkney](https://discourse.mc-stan.org/u/spinkney)\
**Post date:** [April 20, 2021, 3:12pm UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469/9 "2021-04-20T15:12:08Z")

</div>

Does using a `dot_product` speed up the calculation?

```stan
...
Omega[i, j] = dot_product(X[i, j] * M, X[i, j]);
...

```

You can try removing the symmetrize allocation and call outside the loop with `symmetrize_from_lower_tri(Omega)` though I’m not sure you’ll get much from it.

---

<div class="post-metadata">

**Author:** ![mathDR](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mathdr/32/1997_2.png) [@mathDR](https://discourse.mc-stan.org/u/mathDR)\
**Post date:** [April 20, 2021, 6:00pm UTC](https://discourse.mc-stan.org/t/more-efficient-anisotropic-calculation-for-gaussian-process-kernel/21469/10 "2021-04-20T18:00:40Z")

</div>

That helped a bit. Apologies, but re-looking at the different versions I have tried, the fastest (by a considerable margin) was

```
for (i in 1:N-1) {
      Omega[i,i] = 0.0;
      for (j in (i+1):N) {
        Omega[i,j] = quad_form(M, X[i,j]);
        Omega[j,i] = Omega[i,j];
      }
    }
    Omega[N,N] = 0.0;

```

where `X` is now `vector[D] X[N,N]`

I guess I could construct a very large block matrix out of the `X[i,j]` somehow?