CODEWITHSUNDEEP
pythonsympy
slaw
slaw

Reputation: 6869

Speeding Up NumPy Array Generation from SymPy Poly Expression

I have a multivariable polynomial expression that I would like to convert into a square numpy matrix:

import sympy

n = 225
x, C = sympy.symbols(f'x:{n+1}'), sympy.symbols('C')

expr = -1
for i in range(1, n+1):
    expression += x[i]
    
penalty = sympy.Poly(C*(expr)**2)

So, C is just some constant (e.g., C=0.5) and the penalty would be the expansion of:

Ultimately, what I am looking for is the corresponding numpy matrix that represents this expanded expression (minus/dropping the constant term).

Also, note that we assume:

In other words, all linear terms are equal to their squared equivalent (we'll see more of this below).

To provide a more concrete example, let's say n=3 then the penalty will be:

And this expands to:

Due to the relationship/equality of linear term mentioned above, this further simplifies to:

And dropping the constant term, C, would allow us to represent the coefficients of this polynomial as a 3x3 matrix with the squared terms positioned along the diagonal and the cross terms along the off diagonals:

or its equivalent acceptable symmetric form:

And this is the 2D numpy array that I would like to generate but, unfortunately, when n=225 this takes a long time and eventually causes Python to crash. Is there a more efficient approach that I can take for any n that is large (less than 1000)?

Upvotes: 0

Views: 97

Answers (1)

Davide_sd
Davide_sd

Reputation: 13135

SymPy offers many ways to work with polynomials:

  1. symbolic expressions.
  2. the Poly class, which can use two different representations (but that's a long story...).
  3. the ring domains.

Depending on the problem we are going to solve, one approach may be substantially better than the others. Here, with ring domains I'm able to generate the polynomial up to n=1000 in less than two minutes (on my machine). Here is how:

import sympy
from sympy import *
import timeit

n = 1000
x, C = sympy.symbols(f'x:{n+1}'), sympy.symbols('C')

t1 = timeit.default_timer()
expr = -1
for i in range(1, n+1):
    expr += x[i]
poly_expr = C*expr**2

t2 = timeit.default_timer()
print("Step 1 [s]: ", t2 - t1)

# polynomial ring of symbolic coefficients
K = EX[x]
# polynomial built with a polynomial ring
p = K.from_sympy(poly_expr)
t3 = timeit.default_timer()
print("Step 2 [s]: ", t3 - t2)

# new "symbolic" polynomial
new_poly = K.to_sympy(p)

t4 = timeit.default_timer()
print("Step 3 [s]: ", t4 - t3)

# apply assumption xi=xi**2
pows = new_poly.find(Pow)
sd = {p: p.base for p in pows}
new_poly = new_poly.xreplace(sd)
t5 = timeit.default_timer()
print("Step 4 [s]: ", t5 - t4)

# out:
# Step 1 [s]:  1.1785079769997537
# Step 2 [s]:  78.12489726600052
# Step 3 [s]:  26.54897919799987
# Step 4 [s]:  6.294931798999642

As for the generation of the matrix, I believe you have already asked that question a few weeks ago and received the answer...

NOTE: it shouldn't matter in this case because the coefficients of the polynomial are small, but if you regularly works with polynomials you will benefit a lot from installing gmpy2.

Upvotes: 0

Related Questions