I am implementing the Squared Exponential kernel in SageMath, for now let's say it's defined as follows:

$f(x_i,x_k)=σ^2 \exp\left(−\frac{1}{2 \ell^2} \sum_{j=1}^q (x_{i,j} − x_{k,j})^2 \right)$

With $x_i$ and $x_k$ vectors of variable, but equal, length, $\sigma$ and $l$ constant. (In the future $l$ might be vector valued as well, I hope this can be handled then).

The function $f$ must be differentiable in $x_i$ and $x_k$ (not its entries!).

An example implementation might look as follows:

```
n = 8
x1 = vector(list(var('v1_%d' % i) for i in range(1, n+1)))
x2 = vector(list(var('v2_%d' % i) for i in range(1, n+1)))
sigma,l = var('sigma, l')
f = sigma*e^(
sum(vector([(n - m)^2 for n, m in zip(x1.coefficients(), x2.coefficients())]))
*(-1/(2*length_scale)))
```

But what I actually want is this:

```
f(x1, x2) = sigma*e^(sum((x1[i]-x2[i])^2, i, 0, n))*(-1/(2*length_scale)))
```

So that, hopefully, I can do `f.diff(x1)`

and it considers the sum etc. properly.
The problem is that `sum`

needs a variable, but `x1[i]`

demands that `i`

is an integer, so it can grab the elements from the vector.

Is there any obvious method I have missed, or is it not implemented in sage yet? -> If there is no such method, what would be good starting points for me to implement that myself?

**Notes:**
For later processing I need the function to be a differentiable SageMath expression instead of e.g. a Python-function call.
Further, it would be great (but not 100% required) to pass variable length vectors to $f$ instead of specifying the length beforehand.