# Use vector as variable for diff but as values for calculation

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.

edit retag close merge delete

Sort by » oldest newest most voted

One possibility:

def f(a, b, sigma=1, l=1): return sigma^2*exp(-1/(l^2)*sum(map(lambda u,v:(u-v)^2, a, b)))


Example of use :

sage: q=4
sage: X=vector([var("x_{}".format(u)) for u in range(q)])
sage: Y=vector([var("y_{}".format(u)) for u in range(q)])
sage: f(X,Y)
e^(-(x_0 - y_0)^2 - (x_1 - y_1)^2 - (x_2 - y_2)^2 - (x_3 - y_3)^2)
sage: f(X,Y).diff(x_1)
-2*(x_1 - y_1)*e^(-(x_0 - y_0)^2 - (x_1 - y_1)^2 - (x_2 - y_2)^2 - (x_3 - y_3)^2)
# Gradient of f(X,Y) as a function of X
sage: vector([f(X,Y).diff(u) for u in X])
(-2*(x_0 - y_0)*e^(-(x_0 - y_0)^2 - (x_1 - y_1)^2 - (x_2 - y_2)^2 - (x_3 - y_3)^2), -2*(x_1 - y_1)*e^(-(x_0 - y_0)^2 - (x_1 - y_1)^2 - (x_2 - y_2)^2 - (x_3 - y_3)^2), -2*(x_2 - y_2)*e^(-(x_0 - y_0)^2 - (x_1 - y_1)^2 - (x_2 - y_2)^2 - (x_3 - y_3)^2), -2*(x_3 - y_3)*e^(-(x_0 - y_0)^2 - (x_1 - y_1)^2 - (x_2 - y_2)^2 - (x_3 - y_3)^2))


Is that what you seek ?

A "real-world" use might check for equality of lengths of a and b and raise a convenient exception if not satisfied...

more