How to optimize evaluation of hessian in sympy?

Question

I'm using statsmodels library for generic likelihood models. As I have a quite complicated likelihood function, I used sympy to calculate gradient and hessian for me. This works fine, but it is too slow for my needs, because likelihood function contains term b0*x0 + b1*x1 + ... + bn*xn. That way hessian size increases by N^2 and so does the complexity.

Elements of the hessian are often pretty similar like expensive_operation * x0 and expensive_operation * x1, etc. It means that if I could pre-calculate expensive_operation and use it in functions in hessian, I would drastically increase performance.

So the question is - is there a tool which would take list of functions, optimize them and then evaluate them effectively? Something like numexpr which would take list of functions?

Answer 1

SymPy has cse, which stands for common subexpression elimination. See the docs.

A simple example:

>>> print(cse(sin(x**2)*cos(x**2) + 2*sin(x**2) - cos(x**2)))
([(x0, x**2), (x1, sin(x0)), (x2, cos(x0))], [x1*x2 + 2*x1 - x2])