Question: Harmonic Expressions in SymEngine and significant slowdowns when compared to SymPy 1.11.1

[michaelstepniczka]

Hello,

I was wondering about the difference in standards with outputs in SymEngine and SymPy when it comes to evaluating gamma functions and their derivatives. In particular, I am interested in looking at the harmonic form of the gamma/polygamma functions, which I can access in SymPy via e.g. the [.rewrite(sp.harmonic)](https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.gamma_functions.polygamma) function.

I note that I used to need to do this in SymPy version=1.11.1, though it seems to be automatic in the current version 1.13.3. As a minimal example, I am interested in the following:

import sympy as sp

n = sp.symbols('n')
f = sp.gamma(5*n+1)/sp.gamma(n+1)**5
df = sp.diff(f,n)
ddf = sp.diff(df,n)
dddf = sp.diff(ddf,n)

print([sp.expand(dddf.subs({n:x})) for x in range(5)])

In particular, in the current version of SymPy (1.13.3), I do not need to call sp.expand(dddf.subs({n:x}).rewrite(sp.harmonic)). As an aside, as I am looking at SymEngine for optimization purposes, I was curious to see that when calling this style of computation repeatedly many times, it was faster to run in SymPy 1.11.1 and manually call sp.expand(dddf.subs({n:x}).rewrite(sp.harmonic)) as opposed to work in SymPy 1.13.3 and call sp.expand(dddf.subs({n:x})). I found this to be a surprise.

This means that in the current version of SymPy, running:

sp.expand(dddf.subs({n:1}))

outputs

-28800*zeta(3) - 6900 + 7700*pi**2

In the old version (1.11.1), the same code outputs:

-1485715/36 + 7700*pi**2 - 600*polygamma(2, 2) + 15000*polygamma(2, 6)

and it was faster for me to manipulate this expression by directly calling .rewrite(sp.harmonic).

This code can be mapped directly to SymEngine, where in the first code block we simply replace each instance of sp.foo() with se.foo(). I was surprised to note that calling the same command

se.expand(dddf.subs({n:1}))

gave the old version (SymPy 1.11.1) of the output

-1485715/36 + 7700*pi**2 - 600*polygamma(2, 2) + 15000*polygamma(2, 6)

and I was wondering how I should go about manipulating this in SymEngine to find the harmonic form, as I did not directly see a .rewrite() function in the API.

I'll now say that I am interested in this harmonic representation of this expression as I am interested in the rational component of it. We see above that if I don't simplify the polygammas, we would naively say that the rational coefficient is 1485715/36, whereas the number I am interested in is the - 6900 as seen above in the SymPy 13.3.3 evaluation.

I have found a way to access this by using the `.args_as_sympy()' function, and this successfully gave me each argument in terms of the updated version of SymPy where it automatically outputs values in their harmonic representation. Namely, running this on the expression above, we find:

[-1485715/36, -1200 + 1200*zeta(3), 7700*pi**2, 1280515/36 - 30000*zeta(3)]

where it seems the order of the expressions was not preserved. (The order is not important to me). With this, I can then find the rational component of each expression, finding the desired result of:

-1485715/36 -1200 + 0 + 1280515/36 = -6900

My questions are threefold:

• Is there a more natural way to find the harmonic form of the expression directly in SymEngine so I avoid passing to SymPy? I am hopeful that such an implementation will still retain the ~3x advantage in speed that SymPy 1.11.1 had when compared against running the same specific example in SymPy 1.13.3. I do not know how the relative timings change when I push the computation further.
• As I am calling .subs() on a given expression for many different values, is there any simplification I should be making on the expression itself? In SymPy 1.11.1 I was able to directly rewrite the expression using rewrite(sp.harmonic) before evaluation, rather than calling it in every individual instance of evaluation.
  • I did need to be careful to only do this after the derivatives were taken of the gamma functions themselves; otherwise, if I e.g. rewrote f_harm = expr.rewrite(sp.harmonic) and took derivatives of f_harm, evaluation would be left in terms of unevaluated derivatives of harmonic, e.g. \frac{d}{d \xi_{1}} \operatorname{harmonic}{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=5 }}. This is no longer a problem in SymPy 1.13.3, but the code appeared to be slower. The speed results quoted above were without rewriting the expressions themselves, and rather calling either .subs() in 1.13.3 or .subs().rewrite(sp.harmonic) in 1.11.1.
• Is there a best practice around finding the rational value of an expression? I was naively taking the expression, and initially checking if (a) the entire expression is rational or not, and if not, (b) checking if the expression is of type Add. If the expression is of type Add, I went through each argument of the expression running the .is_rational check and took the sum of these to be my result.
  • Therefore, explicitly in the above, -1200 + 1200*zeta(3) is of Add type and so had an is_rational call on both -1200 and 1200*zeta(3), and returned -1200. I'm being sloppy with whether these are SymPy or SymEngine expressions/functions being called, but this doesn't affect the story. (I'll note that for a different project, I was surprised to find that if I declared a variable to be rational, that without evaluating the variable, the .is_rational command returned True).

I'll also note that I haven't actually tried running the analogous code in SymEngine as of yet, and so I can't speak as to what the speedups are. I'm currently rewriting my code to be compatible with both the newer version of SymPy as well as SymEngine. Edit: see the additional comment below. SymEngine seems to be faster than SymPy 1.13.3 in its evaluations, but significantly slower than 1.11.1, which makes me worried that I am not using SymEngine/SymPy as intended in this new version....

Thank you for all the advice! I hope everything is properly documented and the questions are clear and well formulated, but if not, please ask for clarification!

[michaelstepniczka]

Edit:

When I try to evaluate the above expressions with SymEngine or SymPy 1.13.3, I am obtaining significantly slower evaluation times, and I am not sure why this is happening. Explaining a bit again what I said above, I am hoping to find the rational component of the evaluation of these gamma and polygamma functions. Let me express the minimally working examples directly:

import symengine as se
import sympy as sp
import time

def _get_rational_value_sp(expr):
    working_sum = sp.Rational(0,1)
    if expr.is_rational: 
        working_sum=expr
    elif type(expr) == sp.core.add.Add: 
        components = expr.args
        for component in components:
            if component.is_rational:
                working_sum += component
    return working_sum

def _get_rational_value_se(expr):       
    working_sum = se.Rational(0,1)
    if expr.is_Rational: 
        working_sum=expr
    elif type(expr) == se.Add: 
        components = expr.args_as_sympy()
        for component in components:
            working_sum += _get_rational_value_sp(component)
    return working_sum

Then, upon running the following block:

n=se.symbols('n')
f0=se.gamma(5*n+1)/se.gamma(n+1)**5
f1=se.diff(f0,n)
f2=se.diff(f1,n)
f3=se.diff(f2,n)

order = 1000
start_time = time.time()
res0 = [_get_rational_value_se(se.expand(f0.subs({n:x}))) for x in range(order)]
res1 = [_get_rational_value_se(se.expand(f1.subs({n:x}))) for x in range(order)]
res2 = [_get_rational_value_se(se.expand(f2.subs({n:x}))) for x in range(order)]
res3 = [_get_rational_value_se(se.expand(f3.subs({n:x}))) for x in range(order)]
end_time = time.time()
print(f'Total time: {end_time-start_time} seconds to go to order={order}') 

the evaluation time on my machine is Total time: 481.2087595462799 seconds to go to order=1000.

Similarly, if I clear the kernel of my Jupyter notebook and instead run the code in SymPy 1.13.3:

n=sp.symbols('n')
f0=sp.gamma(5*n+1)/sp.gamma(n+1)**5
f1=sp.diff(f0,n)
f2=sp.diff(f1,n)
f3=sp.diff(f2,n)

order = 1000
start_time = time.time()
res0 = [_get_rational_value_sp(sp.expand(f0.subs({n:x}))) for x in range(order)]
res1 = [_get_rational_value_sp(sp.expand(f1.subs({n:x}))) for x in range(order)]
res2 = [_get_rational_value_sp(sp.expand(f2.subs({n:x}))) for x in range(order)]
res3 = [_get_rational_value_sp(sp.expand(f3.subs({n:x}))) for x in range(order)]
end_time = time.time()
print(f'Total time: {end_time-start_time} seconds to go to order={order}') 

the evaluation time on my machine is at least 480 seconds, as I cut off the evaluation manually after I noticed it was slower than the 481 seconds above. A good few minutes into the computation, I started a timer manually and cut if off at 480 seconds; I wouldn't be surprised if this was already a similar amount of time into the computation. In this time, after I aborted the evaluation, I saw that res0, res1, and res2 had completed computation. I wasn't using tqdm and so I don't know the relative evaluation times. Despite the long evaluation times, it seems that the results that were found were correct.

Towards the relative evaluation times, when I commented out res3 = [_get_rational_value_sp(sp.expand(f3.subs({n:x}))) for x in range(order)], I found Total time: 723.2369108200073 seconds to go to order=1000. Similarly, when I also commented out res2 = ..., I found Total time: 168.6326892375946 seconds to go to order=1000. Lastly, when I also commented out res1 = ..., I was left with Total time: 1.3772249221801758 seconds to go to order=1000. Somehow the complexity is in these polygamma functions and their simplifications. I note that for res1, all of the EulerGamma factors cancel out; this is by design. Despite this, however, it seems that it still takes a long time to compute.

However, I remain skeptical of my implementation because when I run the same code in SymPy 1.11.1:

n=sp.symbols('n')
f0=sp.gamma(5*n+1)/sp.gamma(n+1)**5
f1=sp.diff(f0,n)
f2=sp.diff(f1,n)
f3=sp.diff(f2,n)

order = 1000
start_time = time.time()
res0 = [_get_rational_value_sp(sp.expand(f0.subs({n:x}).rewrite(sp.harmonic))) for x in range(order)]
res1 = [_get_rational_value_sp(sp.expand(f1.subs({n:x}).rewrite(sp.harmonic))) for x in range(order)]
res2 = [_get_rational_value_sp(sp.expand(f2.subs({n:x}).rewrite(sp.harmonic))) for x in range(order)]
res3 = [_get_rational_value_sp(sp.expand(f3.subs({n:x}).rewrite(sp.harmonic))) for x in range(order)]
end_time = time.time()
print(f'Total time: {end_time-start_time} seconds to go to order={order}') 

the evaluation time on my machine is: Total time: 27.064862728118896 seconds to go to order=1000.

Question: I believe I must be doing something wrong; there is no reason for SymPy and SymEngine to be slower at evaluating these gamma and polygamma functions than an old version of SymPy. Does anyone have insight as to where SymPy/SymEngine has changed to make the evaluation of these examples worse?

For completeness, when running SymEngine and SymPy 1.13.3, pip list returns:

Package                  Version
------------------------ -----------
absl-py                  1.4.0
aiofiles                 22.1.0
aiosqlite                0.19.0
anyio                    3.6.2
argon2-cffi              21.3.0
argon2-cffi-bindings     21.2.0
arrow                    1.2.3
asttokens                2.2.1
attrs                    23.1.0
Babel                    2.12.1
backcall                 0.2.0
beautifulsoup4           4.12.2
bleach                   6.0.0
certifi                  2022.12.7
cffi                     1.15.1
charset-normalizer       3.1.0
comm                     0.1.3
contourpy                1.0.7
cvxopt                   1.3.0
cycler                   0.11.0
cysignals                1.11.2
Cython                   0.29.34
cytools                  1.2.3
daqp                     0.5.1
debugpy                  1.6.7
decorator                5.1.1
defusedxml               0.7.1
dnspython                2.3.0
ecos                     2.0.12
executing                1.2.0
fastjsonschema           2.16.3
fonttools                4.39.3
fqdn                     1.5.1
gekko                    1.0.6
gmpy2                    2.1.5
h5py                     3.8.0
idna                     3.4
importlib-metadata       6.5.0
importlib-resources      5.12.0
ipykernel                6.22.0
ipython                  8.12.0
ipython-genutils         0.2.0
isoduration              20.11.0
jedi                     0.18.2
Jinja2                   3.1.2
json5                    0.9.11
jsonpointer              2.3
jsonschema               4.17.3
jupyter_client           8.2.0
jupyter_core             5.3.0
jupyter-events           0.6.3
jupyter_server           2.5.0
jupyter_server_fileid    0.9.0
jupyter_server_terminals 0.4.4
jupyter_server_ydoc      0.8.0
jupyter-ydoc             0.2.4
jupyterlab               3.6.3
jupyterlab-pygments      0.2.2
jupyterlab_server        2.22.1
kiwisolver               1.4.4
MarkupSafe               2.1.2
matplotlib               3.7.1
matplotlib-inline        0.1.6
mistune                  2.0.5
Mosek                    10.1.11
mpmath                   1.3.0
nbclassic                0.5.5
nbclient                 0.7.3
nbconvert                7.3.1
nbformat                 5.8.0
nest-asyncio             1.5.6
notebook                 6.5.4
notebook_shim            0.2.2
numpy                    1.24.2
ortools                  9.6.2534
osqp                     0.6.2.post9
packaging                23.1
pandocfilters            1.5.0
parso                    0.8.3
pexpect                  4.8.0
pickleshare              0.7.5
Pillow                   9.5.0
pip                      23.0.1
platformdirs             3.2.0
pplpy                    0.8.7
prometheus-client        0.16.0
prompt-toolkit           3.0.38
protobuf                 4.22.3
psutil                   5.9.5
ptyprocess               0.7.0
pure-eval                0.2.2
pycparser                2.21
Pygments                 2.15.1
pymongo                  4.3.3
pyparsing                3.0.9
pyrsistent               0.19.3
python-dateutil          2.8.2
python-flint             0.6.0
python-json-logger       2.0.7
PyYAML                   6.0
pyzmq                    25.0.2
qdldl                    0.1.7
qpsolvers                3.4.0
requests                 2.28.2
rfc3339-validator        0.1.4
rfc3986-validator        0.1.1
scikit-sparse            0.4.12
scipy                    1.10.1
scs                      3.2.3
Send2Trash               1.8.0
setuptools               65.5.0
six                      1.16.0
sniffio                  1.3.0
soupsieve                2.4.1
stack-data               0.6.2
symengine                0.13.0
sympy                    1.13.3
terminado                0.17.1
tinycss2                 1.2.1
tomli                    2.0.1
tornado                  6.3
tqdm                     4.65.0
traitlets                5.9.0
typing_extensions        4.5.0
uri-template             1.2.0
urllib3                  1.26.15
wcwidth                  0.2.6
webcolors                1.13
webencodings             0.5.1
websocket-client         1.5.1
y-py                     0.5.9
ypy-websocket            0.8.2
zipp                     3.15.0
WARNING: There was an error checking the latest version of pip.
Note: you may need to restart the kernel to use updated packages.

When I run with the version SymPy 1.11.1:

Package                  Version
------------------------ -----------
absl-py                  1.4.0
aiofiles                 22.1.0
aiosqlite                0.19.0
anyio                    3.6.2
argon2-cffi              21.3.0
argon2-cffi-bindings     21.2.0
arrow                    1.2.3
asttokens                2.2.1
attention_grabber        0.1.3
attrs                    23.1.0
Babel                    2.12.1
backcall                 0.2.0
beautifulsoup4           4.12.2
bleach                   6.0.0
certifi                  2022.12.7
cffi                     1.15.1
charset-normalizer       3.1.0
comm                     0.1.3
contourpy                1.0.7
cvxopt                   1.3.0
cycler                   0.11.0
cysignals                1.11.2
Cython                   0.29.34
cytools                  1.2.3
daqp                     0.5.1
debugpy                  1.6.7
decorator                5.1.1
defusedxml               0.7.1
dnspython                2.3.0
ecos                     2.0.12
executing                1.2.0
fastjsonschema           2.16.3
fonttools                4.39.3
fqdn                     1.5.1
galois                   0.3.7
gekko                    1.0.6
gmpy2                    2.1.5
gurobipy                 11.0.2
h5py                     3.8.0
idna                     3.4
importlib-metadata       6.5.0
importlib-resources      5.12.0
ipykernel                6.22.0
ipython                  8.12.0
ipython-genutils         0.2.0
ipywidgets               7.7.1
isoduration              20.11.0
jedi                     0.18.2
Jinja2                   3.1.2
json5                    0.9.11
jsonpointer              2.3
jsonschema               4.17.3
jupyter_client           8.2.0
jupyter_core             5.3.0
jupyter-events           0.6.3
jupyter_server           2.5.0
jupyter_server_fileid    0.9.0
jupyter_server_terminals 0.4.4
jupyter_server_ydoc      0.8.0
jupyter-ydoc             0.2.4
jupyterlab               3.6.3
jupyterlab-pygments      0.2.2
jupyterlab_server        2.22.1
jupyterlab-widgets       1.1.1
kiwisolver               1.4.4
llvmlite                 0.41.1
MarkupSafe               2.1.2
matplotlib               3.7.1
matplotlib-inline        0.1.6
mistune                  2.0.5
Mosek                    10.1.11
mpmath                   1.3.0
multiset                 3.2.0
nbclassic                0.5.5
nbclient                 0.7.3
nbconvert                7.3.1
nbformat                 5.8.0
nest-asyncio             1.5.6
networkx                 3.0
notebook                 6.5.4
notebook_shim            0.2.2
numba                    0.58.1
numpy                    1.24.2
ortools                  9.6.2534
osqp                     0.6.2.post9
packaging                23.1
pandas                   2.2.2
pandocfilters            1.5.0
parso                    0.8.3
patsy                    0.5.6
pexpect                  4.8.0
pickleshare              0.7.5
Pillow                   9.5.0
pip                      23.0.1
platformdirs             3.2.0
plotly                   5.13.1
plotly-express           0.4.1
pplpy                    0.8.7
pprofile                 2.1.0
prometheus-client        0.16.0
prompt-toolkit           3.0.38
protobuf                 4.22.3
psutil                   5.9.5
ptyprocess               0.7.0
pure-eval                0.2.2
pycparser                2.21
Pygments                 2.15.1
pymongo                  4.3.3
Pympler                  1.0.1
pyparsing                3.0.9
pyrsistent               0.19.3
python-dateutil          2.8.2
python-flint             0.3.0
python-json-logger       2.0.7
pytz                     2024.1
PyYAML                   6.0
pyzmq                    25.0.2
qdldl                    0.1.7
qpsolvers                3.4.0
requests                 2.28.2
rfc3339-validator        0.1.4
rfc3986-validator        0.1.1
scikit-sparse            0.4.12
scipy                    1.10.1
scs                      3.2.3
Send2Trash               1.8.0
setuptools               65.5.0
six                      1.16.0
sniffio                  1.3.0
soupsieve                2.4.1
stack-data               0.6.2
statsmodels              0.14.2
symengine                0.13.0
sympy                    1.11.1
tenacity                 9.0.0
terminado                0.17.1
tinycss2                 1.2.1
tomli                    2.0.1
tornado                  6.3
tqdm                     4.65.0
traitlets                5.9.0
typing_extensions        4.5.0
tzdata                   2024.1
uri-template             1.2.0
urllib3                  1.26.15
wcwidth                  0.2.6
webcolors                1.13
webencodings             0.5.1
websocket-client         1.5.1
widgetsnbextension       3.6.8
y-py                     0.5.9
ypy-websocket            0.8.2
zipp                     3.15.0
WARNING: There was an error checking the latest version of pip.
Note: you may need to restart the kernel to use updated packages.

For reasons unrelated to this project, I am running this all within two different Docker images, with the installed packages as above. I don't see any reason why such significant slowdowns would come from the minimal changes here. I hope I have included enough information; I apologize any is extraneous.