I am currently solving the following problem (Gaussian tail approximation) using solve() :
x, y = symbols('x y')
l_bound = np.mean(self.data[j]) - 3 * np.std(self.data[j]) # e.g. 4
u_bound = np.mean(self.data[j]) + 3 * np.std(self.data[j]) # e.g. 13
sol_dict = solve([y/(x-1) - 3 * y/sqrt(((x-1)**2)*(x-2))-l_bound, y/(x-1) + 3 * y/sqrt(((x-1)**2) * (x-2))-u_bound], [x, y], dict=True)
I am trying to find a way to speed up the process, because it takes forever and I have quite many rows of a matrix of data to process (each j is a row). The examples I found in the documentaries show how to evaluate a function - Is there a possibility to use solve() on GPU?
What I've done so far: I've read through documentaries, unable to find the right expression. - As I am new to the topic, I might have missed some essential part (sorry if my question might be trivial). Thanks in advance for any pointer in a more suited direction, appreciate any help.
Sympy is a CAS based purely on Python, so I don't believe that complex computations such a solvers can easily be ported to work on a GPU.
However, by looking at your example it appears that you are trying to get numerical solutions. Hence, solve is not the right tool for the job. nsolve (numerical solve) might be better suited for the job. Generally, you will need to provide appropriate initial guesses. For example:
l_bound = 4
u_bound = 13
expr1 = y/(x-1) - 3 * y/sqrt(((x-1)**2)*(x-2))-l_bound
expr2 = y/(x-1) + 3 * y/sqrt(((x-1)**2) * (x-2))-u_bound
nsolve([expr1, expr2], [x, y], [10, 10])
# out: (34.1111111111111, 281.444444444444)
You'll notice that for this particular case, nsolve is much faster than solve. However, for each iteration of your loop, nsolve will convert the system of symbolic equations to a system of numerical functions using lambdify: this operation is expensive. If the form of your equations is not changing with iterations, than we can convert the system of symbolic equations to a system of numerical functions once, at the beginning, and reuse that inside the loop with updated values. For example:
from mpmath import findroot
l_bound, u_bound = symbols("l u")
expr1 = y/(x-1) - 3 * y/sqrt(((x-1)**2)*(x-2))-l_bound
expr2 = y/(x-1) + 3 * y/sqrt(((x-1)**2) * (x-2))-u_bound
system = Matrix([expr1, expr2]).T
J = system.jacobian([x, y])
# convert symbolic objects to numerical functions
system = lambdify([x, y, l_bound, u_bound], system.T, modules="mpmath")
J = lambdify([x, y, l_bound, u_bound], J, modules="mpmath")
l_bounds, u_bounds = [4], [13]
IG = [10, 10]
for l, u in zip(l_bounds, u_bounds):
# provide new values of lower and upper bounds with each iteration
sol = findroot(lambda x, y: system(x, y, l, u), IG, J=lambda x, y: J(x, y, l, u))
# continue your task
Notes:
nsolve is going to do exactly what we have done in the previous code block: convert the system of equations to a symbolic matrix, compute a symbolic jacobian, converts both system and Jacobian to numerical functions and use mpmath's findroot to find the solution.findroot you use some other root finding algorithm from Scipy. That's for you to find out :)