Gradient function evaluation

Below is some code I wrote to evaluate the position of a point moving towards the minimum of a gradient/3d function (defined at the beginning as "eq"). roll.roll() does this by repeatedly evaluating the equation at point (x,y), moving it in the direction of the gradient, then repeating with the new point.

It is very very slow to run though. I think this is because either calculate() is inefficient, or sympy's symbolic equation manipulation in roll.roll is really slow. Does anyone have any ideas on how to speed this up? Is there another libray other than SymPy that is faster?

import sympy as smp

x, y = smp.symbols('x y')
eq = 1*smp.exp(-((x-5)/5)**2 - ((y-1)/2)**2) + \
    2*smp.exp(-((x+3)/2)**2 - ((y-3)/2)**2) + \
    3*smp.exp(-((x-4)/2)**2 - ((y-7)/2)**2)

# Evaluates the 2 input sympy symbolic function "expression" at points (x1,y1)
def calculate(expression,x1,y1):
    EQ = smp.lambdify((x,y), expression, 'numpy')
    return EQ(x1,y1)

class roll:
    xDiff = smp.diff(eq,x)
    yDiff = smp.diff(eq,y)
    normalize = eq/smp.sqrt(xDiff**2 + yDiff**2)

    def roll(x,y,duration):
        (x,y) = (x,y)
        for i in range(0,duration):
            (x,y) = (
                x-calculate((roll.normalize*roll.xDiff),x,y),
                y-calculate((roll.normalize*roll.yDiff),x,y)
                )
        return (x,y)

print(roll.roll(1,2,10))

Here is a visual to help see what this program is doing; the bigger the colored dots are, the greater the function f(x) is evaluated at that point. The draggable point represents what the program is attempting to find. https://www.desmos.com/calculator/c8mq2rijqn

I've tried to figure out if it's possible to pre-calculate normalize*xDiff not inside of roll.roll, but idk if thats possible.

Also, I believe that it is actually pretty easy to do this if the step size isn't dependent on the value of the function at the current point. I do need it to move faster when it's at a high point on the graph though (not just a point with a steep slope) so that has really been hard to figure out too.

You are calling lambdify inside the loop. The point of lambdify is that it returns a fast function but lambdify itself is a lot slower than the function that it returns. You should call lambdify once and then in the loop repeatedly use the function that was returned by it.

This code is equivalent to yours and returns the exact same result but the loop is 500x times faster:

import sympy as smp

x, y = smp.symbols('x y')
eq = 1*smp.exp(-((x-5)/5)**2 - ((y-1)/2)**2) + \
    2*smp.exp(-((x+3)/2)**2 - ((y-3)/2)**2) + \
    3*smp.exp(-((x-4)/2)**2 - ((y-7)/2)**2)

# Evaluates the 2 input sympy symbolic function "expression" at points (x1,y1)
def calculate(expression,x1,y1):
    EQ = smp.lambdify((x,y), expression, 'numpy')
    return EQ(x1,y1)

class roll:
    xDiff = smp.diff(eq,x)
    yDiff = smp.diff(eq,y)
    normalize = eq/smp.sqrt(xDiff**2 + yDiff**2)

    # call lambdify once
    fxy = smp.lambdify((x, y), (x - normalize*xDiff, y - normalize*yDiff))

    def roll(x,y,duration):
        (x,y) = (x,y)
        for i in range(0,duration):
            # in the loop call the function that was returned by lambdify
            x, y = roll.fxy(x, y)
        return (x,y)

print(roll.roll(1,2,10))

The question is you want to do gradient-descent.

If you define a function about grad value, it's the most efficient way.

def eq(x,y): 
    return y*x**2+2*x*y
def grad_x(x,y):
    return 2*x*y+2*y
def grad_y(x,y):
    return x**2+2*x

But when you confused about how to calculate the grad of your function, you can use other package which can autograd (i.e. numpy, pytorch).

Here is an example of pytorch:

import torch

def eq(x, y):
    return (
        1 * torch.exp(-(((x - 5) / 5) ** 2) - ((y - 1) / 2) ** 2)
        + 2 * torch.exp(-(((x + 3) / 2) ** 2) - ((y - 3) / 2) ** 2)
        + 3 * torch.exp(-(((x - 4) / 2) ** 2) - ((y - 7) / 2) ** 2)
    )


def roll(x, y, duration):
    x, y = (
        torch.tensor(x).requires_grad_(True),
        torch.tensor(y).requires_grad_(True),
    )
    for _ in range(0, duration):
        func_value = eq(x, y)
        xDiff = torch.autograd.grad(func_value, x, retain_graph=True)[0]
        yDiff = torch.autograd.grad(func_value, y)[0]
        normalize = func_value / torch.sqrt(xDiff ** 2 + yDiff ** 2)
        x = x - normalize * xDiff
        y = y - normalize * yDiff
    return (x, y)


print(roll(1.0, 2.0, 10))