hyint: hybrid (symbolic-numeric) integration

[Shahriar Iravanian]

Hi All,

I have uploaded a new package, hyint, to github and PyPi. It is a hybrid (symbolic-numeric) integration package on top of SymPy and numpy/scipy. I would appreciate it if it is added to the list of projects using SymPy.

You can find it at https://github.com/siravan/hyint or install it as 'pip install hyint'.

hyint uses an ansatz generation algorithm similar to the Risch-Bronstein poor man's integrator combined with a sparse regression algorithm 
adopted from the Sparse identification of nonlinear dynamics (SINDy) to solve indefinite integration problems to univariate expressions. It can solve a large subset of standard  elementary integrals despite a very small size (a few hundred lines of code).

I'm one of the principal contributors to SymbolicNumericIntegartion.jl (https://github.com/SciML/SymbolicNumericIntegration.jl), which is a Julia package for symbolic-numeric integration and is the basis of hyint.

Thanks,

Shahriar Iravanian

[Oscar Benjamin]

On Tue, 29 Aug 2023 at 07:09, Shahriar Iravanian <irvani...@gmail.com> wrote:
>
> Hi All,
>

Hi Shahriar,

This sounds excellent. Yes, you can add it to the list of projects
using SymPy using a pull request, I think to the website repo.

How exactly is it different from Bronstein's poor man's integrator?

SymPy has an integration algorithm called heurisch which is based on
the poor man's integrator but uses exact rather than approximate
solutions to the linear equations. Is exact vs approximate the main
distinction here between what hyint does and what heurisch does?

--
Oscar

[Aaron Meurer]

Judging from the examples, it looks like it

In: hyint.integrate(1 / (x**3 - 2*x + 1), x)
Out: 2.17082039324994*log(x - 1) + 1.34164078649988*log(x + 1/2 +
sqrt(5)/2) - 1.17082039324993*log(x**3 - 2*x + 1)

And looking at the code, I'd say a difference is that the Bronstein
algorithm is more rigorous and general in how it generates its ansatz.
But actually one of the reasons heurisch() is so slow is that it can
generate ansatz that are huge (especially when it's given an integrand
that doesn't have an answer to begin with). Employing some heuristics
like the ones used here in heurisch could help it out.

Actually it wouldn't be much work to generalize this to something that
gives an exact answer at least some of the time. For example, you
could use nsimplify() to guess at exact values for numeric
coefficients. As long as diff(answer) - integrand symbolically
simplifies to 0 you know you have a correct answer.

Aaron Meurer

>
> --
> Oscar
>
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxRyJ3FUZ1pN6GACMoHZBYGbUJAavKq5Y1RW5__1nFAT8w%40mail.gmail.com.

[Shahriar Iravanian]

Thank you very much for your comments. 

While hyint ansatz generation is similar in spirit to heurisch, it is less rigorous and more 

heuristic. One reason is that the numerical portion of hyint is much more forgiving and 

can filter out many undesirable ansatzes (by using QR decomposition to remove linearly 

dependent ones and sparse regression to solve the integral). Therefore, we can generate 

many different forms of ansatz, which allows hyint to solve integrals beyond what 

heurisch can do. 

For example, I use a set of 170 test cases (listed in test.py) to test hyint. The current version 

can solve 154 of them; heurisch solves 126 (of course, this list is curated for hyint). 

Some integrals that hyint can solve and heurisch cannot:

x/sqrt(x - 1) => 2*x*sqrt(x - 1)/3 + 4*sqrt(x - 1)/3
1/sqrt(x**2 + 4) => asinh(x/2)
sqrt(x)*log(x) => 2*x**(3/2)*log(x)/3 - 4*x**(3/2)/9
log(log(x))/x => log(x)*log(log(x)) - log(x)
log(cos(x))*tan(x) => -log(cos(x))**2/2
sin(x + 1)/(x + 1) => Si(x + 1)

sqrt(1 - sin(x)) => 2*cos(x)/sqrt(1 - sin(x))

Some integrals that heurisch can solve and hyint cannot:

exp(x)/(exp(2*x) - 1) => log(exp(x) - 1)/2 - log(exp(x) + 1)/2
tan(x)**3 => -log(tan(x)**2 + 1)/2 + tan(x)**2/2

> Actually it wouldn't be much work to generalize this to something that

> gives an exact answer at least some of the time.

I completely agree. I believe the main application of hyint is as an ansatz discovery

engine rather than a standalone integrator. This is a line of research I'm actively

pursuing and would appreciate help, comments, and collaboration!

-- Shahriar 

[Aaron Meurer]

On Thu, Aug 31, 2023 at 8:22 PM Shahriar Iravanian
<irvani...@gmail.com> wrote:
>
> Thank you very much for your comments.
>
> While hyint ansatz generation is similar in spirit to heurisch, it is less rigorous and more
> heuristic. One reason is that the numerical portion of hyint is much more forgiving and
> can filter out many undesirable ansatzes (by using QR decomposition to remove linearly
> dependent ones and sparse regression to solve the integral). Therefore, we can generate
> many different forms of ansatz, which allows hyint to solve integrals beyond what
> heurisch can do.
>
> For example, I use a set of 170 test cases (listed in test.py) to test hyint. The current version
> can solve 154 of them; heurisch solves 126 (of course, this list is curated for hyint).
>
> Some integrals that hyint can solve and heurisch cannot:
>
> x/sqrt(x - 1) => 2*x*sqrt(x - 1)/3 + 4*sqrt(x - 1)/3
> 1/sqrt(x**2 + 4) => asinh(x/2)
> sqrt(x)*log(x) => 2*x**(3/2)*log(x)/3 - 4*x**(3/2)/9
> log(log(x))/x => log(x)*log(log(x)) - log(x)
> log(cos(x))*tan(x) => -log(cos(x))**2/2
> sin(x + 1)/(x + 1) => Si(x + 1)
> sqrt(1 - sin(x)) => 2*cos(x)/sqrt(1 - sin(x))

By the way, a big problem with heurisch is that it's very sensitive to
the way the derivative of the expression is represented, since that's
how it generates the ansatz. This is fundamental to the approach, but
an issue is that it just directly uses whatever sympy.diff() returns.
This is a big problem for algebraic functions because they can be
represented in multiple different ways, and that representation can
affect the result of integration (a big piece of the "real" Risch
algorithm for algebraic functions is normalizing things so that this
isn't a problem). For example, the derivative of s = sqrt(x) can be
represented as 1/(2*s) or s/(2*x). SymPy happens to use the former
form, because of the way the core combines exponents, but there are
integrals with square roots where heurisch can't find the answer but
would if it used the latter form (not so much because the latter form
is better, although in some ways it is since it's a polynomial in s,
but usually just by dumb luck of the way it generates the ansatz).

>
> Some integrals that heurisch can solve and hyint cannot:
>
> exp(x)/(exp(2*x) - 1) => log(exp(x) - 1)/2 - log(exp(x) + 1)/2
> tan(x)**3 => -log(tan(x)**2 + 1)/2 + tan(x)**2/2
>
> > Actually it wouldn't be much work to generalize this to something that
> > gives an exact answer at least some of the time.
>
> I completely agree. I believe the main application of hyint is as an ansatz discovery
> engine rather than a standalone integrator. This is a line of research I'm actively
> pursuing and would appreciate help, comments, and collaboration!

I like the idea of using some faster linear algebra as a preprocessor
to reduce the size of the problem to be solved. You could then take
the reduced problem and solve it with exact symbolic arithmetic. The
worst that can happen is you might drop an ansatz that shouldn't have
been dropped, but you would still always give a symbolically correct
answer when you do.

If heurisch could support much larger, possibly linearly dependent
ansatz, then that would open up a lot of possibilities, like for
instance, trying to use ansatz coming from both sqrt(x) and x/sqrt(x)
simultaneously. But right now it's way too slow even with just the
linear system it generates.

I'm also wondering if you've thought about symbolic constants and if
there's any tricks you could do to support them. I think there might
be, especially when still thinking about things in terms of just using
faster linear algebra as a pre-processor to a symbolic solve.

Aaron Meurer

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/c0d397a1-640f-477e-9897-e1157bb6ca67n%40googlegroups.com.

[Oscar Benjamin]

On Fri, 1 Sept 2023 at 06:45, Aaron Meurer <asme...@gmail.com> wrote:
>
> I like the idea of using some faster linear algebra as a preprocessor
> to reduce the size of the problem to be solved. You could then take
> the reduced problem and solve it with exact symbolic arithmetic. The
> worst that can happen is you might drop an ansatz that shouldn't have
> been dropped, but you would still always give a symbolically correct
> answer when you do.
>
> If heurisch could support much larger, possibly linearly dependent
> ansatz, then that would open up a lot of possibilities, like for
> instance, trying to use ansatz coming from both sqrt(x) and x/sqrt(x)
> simultaneously. But right now it's way too slow even with just the
> linear system it generates.

I'm not sure if I've interpreted this point correctly but when
heurisch is slow it is not usually because of the linear algebra part.
Probably that is the easiest part to speed up as well.

The slow parts of heurisch are converting back and forth between Expr
and the domains, differentiation, and cancel.

If heurisch was being rewritten then I think it could be written
without the back and forth conversions. The differentiation could be
handled in the domain representation. The cancellation could be
handled there as well.

I suspect that heurisch was implemented at a time when SymPy did not
have as many of the pieces that it now has that you would want to use
to implement the algorithm within the algebra subsystem rather than
the symbolic subsystem.

Possibly the limiting factor for SymPy implementing heurisch in the
ideal way right now is support for symbolic algebraic extensions like
Q(x)[sqrt(x)] within the domain system although polys.agca has
FiniteMonogenicExtension that could represent that.

--
Oscar

[Aaron Meurer]

On Fri, Sep 1, 2023 at 4:12 AM Oscar Benjamin
<oscar.j....@gmail.com> wrote:
>
> On Fri, 1 Sept 2023 at 06:45, Aaron Meurer <asme...@gmail.com> wrote:
> >
> > I like the idea of using some faster linear algebra as a preprocessor
> > to reduce the size of the problem to be solved. You could then take
> > the reduced problem and solve it with exact symbolic arithmetic. The
> > worst that can happen is you might drop an ansatz that shouldn't have
> > been dropped, but you would still always give a symbolically correct
> > answer when you do.
> >
> > If heurisch could support much larger, possibly linearly dependent
> > ansatz, then that would open up a lot of possibilities, like for
> > instance, trying to use ansatz coming from both sqrt(x) and x/sqrt(x)
> > simultaneously. But right now it's way too slow even with just the
> > linear system it generates.
>
> I'm not sure if I've interpreted this point correctly but when
> heurisch is slow it is not usually because of the linear algebra part.
> Probably that is the easiest part to speed up as well.

Maybe things have changed since the introduction of things like
DomainMatrix, but when I benchmarked it years ago this was the case.
That and the fact that it was using the dense poly representation with
hundreds of variables.

>
> The slow parts of heurisch are converting back and forth between Expr
> and the domains, differentiation, and cancel.
>
> If heurisch was being rewritten then I think it could be written
> without the back and forth conversions. The differentiation could be
> handled in the domain representation. The cancellation could be
> handled there as well.

We could probably merge a lot of the code in risch and heurisch, but
loosening the restrictive rigorous code from risch. Actually Risch
itself can solve some algebraic integrals if the restrictions are
removed (just the nonelementary integral errors would be incorrect).

>
> I suspect that heurisch was implemented at a time when SymPy did not
> have as many of the pieces that it now has that you would want to use
> to implement the algorithm within the algebra subsystem rather than
> the symbolic subsystem.

Yes, heurisch predates the polys module. I integrated it with the
polys back in 2012 or 2013, but I didn't really refactor it in any
major way.

Aaron Meurer

>
> Possibly the limiting factor for SymPy implementing heurisch in the
> ideal way right now is support for symbolic algebraic extensions like
> Q(x)[sqrt(x)] within the domain system although polys.agca has
> FiniteMonogenicExtension that could represent that.
>
> --
> Oscar
>

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxSN%2BVBubjX6i3mq%3DwE7nXZSorg7drnSP_X3-2daxbzE0Q%40mail.gmail.com.

[Shahriar Iravanian]

> This is a big problem for algebraic functions because they can be

> represented in multiple different ways, and that representation can
> affect the result of integration (a big piece of the "real" Risch
> algorithm for algebraic functions is normalizing things so that this
> isn't a problem). For example, the derivative of s = sqrt(x) can be
> represented as 1/(2*s) or s/(2*x).

The strength of the numerical method is that the distinction between 

1/(2*s) and s/(2*x) is lost. We can have both expressions in the ansatz 

list, and one will be dropped because they are linearly dependent.

> I like the idea of using some faster linear algebra as a preprocessor
> to reduce the size of the problem to be solved. You could then take
> the reduced problem and solve it with exact symbolic arithmetic. The
> worst that can happen is you might drop an ansatz that shouldn't have
> been dropped, but you would still always give a symbolically correct
> answer when you do.

This is the plan. We can use hyint to prune the set of ansatzes to those 

actually in the integral and then solve symbolically for the coefficients.

> I'm also wondering if you've thought about symbolic constants and if
> there's any tricks you could do to support them. I think there might
> be, especially when still thinking about things in terms of just using
> faster linear algebra as a pre-processor to a symbolic solve.

The basic idea is to generate ansatzes with intact symbolic constants, but
then substitute them with random numbers (same as we do for x), solve 

the integral with sparse regression, and then, starting with the 

pruned set of symbolic ansatzes, solve for the coefficients as above.

Regarding improving the speed of calculations, I cannot comment on 

sympy linear solvers and polys module but should mention that 

symbolic-numeric methods can also improve the speed. Currently, 

hyint separates the symbolic part (ansatz generation) from the numerical 

part. But it is possible to mix numerical methods in ansatz generation. 

We can calculate only the diff of factors of the integrands and not the 

whole expressions. Afterward, ansatz generation can proceed mainly 

numerically (while keeping track of the path for each ansatz; polys may 

be helpful here). Finally, we can reconstruct the relevant ansatz 

symbolically. This way, we can skip most of the symbolic diffs and 

cancellations.

Thanks,

Shahriar Iravanian  

[Shahriar Iravanian]

> I like the idea of using some faster linear algebra as a preprocessor
> to reduce the size of the problem to be solved. You could then take
> the reduced problem and solve it with exact symbolic arithmetic. The
> worst that can happen is you might drop an ansatz that shouldn't have
> been dropped, but you would still always give a symbolically correct
> answer when you do.

I'm testing this idea and it seems to be working fine. I just uploaded a new 

version of hyint that exports a `hints(eq, x)` function to return a filtered list 

of ansatz for eq. This version accepts expressions with symbolic 

constants. 

The output of hyint.hints() can then be fed to the `hints` argument of 

heurisch as a preprocessor. Example:

>>> from sympy import log
>>> from sympy.integrals.heurisch import heurisch
>>> from sympy.abc import a, x
>>> import hyint  

>>> y = log(log(x) + a) / x

>>> heurisch(y, x)
# returns None, no solution

>>> hints = hyint.hints(y, x)
>>> print(hints)
[log(a + log(x)), log(x)*log(a + log(x)), log(x)]
   
>>> heurisch(y, x, hints=hints)
a*log(a + log(x)) + log(x)*log(a + log(x)) - log(x)
    

Thanks,

Shahriar Iravanian

On Friday, September 1, 2023 at 3:58:27 PM UTC-4 asme...@gmail.com wrote: