Solvers improvement roadmap

[Sayandip Halder]

Hello! I hope to discuss the possible action plan for improving solvers and solveset here. The GitHub wiki is old and needs to be updated.

The following are what I can think of right now. Also, I want to know what people in the Sympy community think of this.

1. Solve Heaviside, Piecewise and floor functions correctly with solve()

2. Improve solving trigonometric functions with solveset()

3. Add linsolve, nonlinsolve and solve_decomposition to solveset()

4. Handler for union of ImageSets (PR 18489)

[Oscar Benjamin]

On Mon, 15 Feb 2021 at 21:01, Sayandip Halder <sayan...@gmail.com> wrote:
>
> Hello! I hope to discuss the possible action plan for improving solvers and solveset here. The GitHub wiki is old and needs to be updated.

Please do update it!

> The following are what I can think of right now. Also, I want to know what people in the Sympy community think of this.
>
> 1. Solve Heaviside, Piecewise and floor functions correctly with solve()

It might make more sense to fix these things in solveset and then make
use of solveset in solve rather than trying to fix things in solve
itself.

> 2. Improve solving trigonometric functions with solveset()

> 4. Handler for union of ImageSets (PR 18489)

I've put the above two points together. It would be good to fix these.
I think it should be straight-forward to make some progress.

> 3. Add linsolve, nonlinsolve and solve_decomposition to solveset()

I have recently made significant performance improvements in linsolve.
It would be good to make more use of linsolve in the other solvers.
I'm not sure what you mean for the other points because solveset is
for univariate equations and the others are multivariate solvers.

I think that a major point that none of the current solvers can
address is the conditional existence of solutions:
https://github.com/sympy/sympy/issues/16861
That's a significant problem for many users who want sympy to find
solutions to generically inconsistent systems e.g.:
https://github.com/sympy/sympy/issues/19859

Another point is that many users just want numeric solutions and want
to be able to get them more conveniently e.g.:
https://github.com/sympy/sympy/issues/19164

Another problem is how to make use of the outputs from solveset:
https://github.com/sympy/sympy/issues/10006

The diophantine solvers are awkward to use but could easily be made better e.g.:
https://github.com/sympy/sympy/issues/20682

We also need solvers for handling (at least linear) inequalities. The
current inequality solvers are severely lacking.

A major problem is that we have a mess of multiple solvers that have
overlapping functionality:

1. solve is actually the interface that most users want but the solve
function has an unreliable API and legacy behaviour
2. solveset is nice and formally correct but it's hard to use the
output and only for univariate equations
3. solvify tries to wrap solveset to give the output that solve does
but it doesn't cover as many cases
4. linsolve/nonlinsolve are reasonable but not feature complete and
also have an awkward return type

None of these functions really replaces any other so we just have a
range of different overlapping implementations and inconstent APIs. It
would be good to unify as much as possible of the implementations and
make the APIs consistent and more useful. Making the different
routines share a common clean codebase is needed but is a lot of work.

The most immediate thing that is needed is good documentation that
explains why there are different solvers and when you might want to
use each of them and discusses the caveats and limitations etc.

--
Oscar

[Sayandip Halder]

On Tuesday, February 16, 2021 at 5:21:21 AM UTC+5:30 Oscar wrote:

It might make more sense to fix these things in solveset and then make
use of solveset in solve rather than trying to fix things in solve
itself.

There are many instances where solveset() will return a range. Consider this: p = Piecewise((0, x < -1), (x**2, x <= 1), (log(x), True)).

Solveset will return an Interval for the region x<-1. If we use solveset in solve, then the behaviour of solve will change.

 

 

[Oscar Benjamin]

Currently solve gives an error for that case:
NotImplementedError: solve cannot represent interval solutions

If solve were to use solveset then it could just check for an interval
and raise the same error.

--
Oscar

[Sayandip Halder]

Your PR #18814 made solve_linear_system() a wrapper for linsolve(). But that was later replaced in some PR (which I am unable to find) by solve_lin_sys(). Why was it changed? Ref: here

-

Sayandip Halder

[Oscar Benjamin]

On Sat, 13 Mar 2021 at 22:23, Sayandip Halder <sayan...@gmail.com> wrote:
>
>
> Your PR #18814 made solve_linear_system() a wrapper for linsolve(). But that was later replaced in some PR (which I am unable to find) by solve_lin_sys(). Why was it changed? Ref: here

At the time it turned out that linsolve was slow for some cases where
solve_linear_system was faster. Since then I introduced a new matrix
class (DomainMatrix) that is used internally for some operations and
is much faster than the current matrix class in many common cases.
That is now used by linsolve and solve (and also solve_lin_sys).

There were a bunch of changes beginning with:
https://github.com/sympy/sympy/pull/18814
https://github.com/sympy/sympy/pull/18844

Most recent change:
https://github.com/sympy/sympy/pull/20780

See here for the current status:
https://github.com/sympy/sympy/issues/20987

I think that building out the implementation of DomainMatrix and using
it to make sympy much faster is a major priority for sympy right now.

Oscar

[Sayandip Halder]

I have recently made significant performance improvements in linsolve.
It would be good to make more use of linsolve in the other solvers.

Do you mean this PR?

I'm not sure what you mean for the other points because solveset is
for univariate equations and the others are multivariate solvers.

I was referring to this.

 

4. linsolve/nonlinsolve are reasonable but not feature complete and
also have an awkward return type

I am not sure what you mean by awkward return type. The FiniteSet and ImageSet returned by linsolve() and nonlinsolve() respectively

are the same as the return type of solveset(). As for feature complete, what more features do you think should be added to them?

#19038, #17566, #16861, PR 12424 are the most important issues and PRs that should be tackled at first, in my opinion. Do you have any more in mind that

should be added to this list?

 

None of these functions really replaces any other so we just have a
range of different overlapping implementations and inconstent APIs. It
would be good to unify as much as possible of the implementations and
make the APIs consistent and more useful.

 

I hope to take up this as my GSOC project. Will there be any mentor for this?

The most immediate thing that is needed is good documentation that
explains why there are different solvers and when you might want to
use each of them and discusses the caveats and limitations etc.

Agreed. We should start with solve(), which is the most lacking when it comes to documentation.

 

--
Oscar

[Oscar Benjamin]

On Mon, 15 Mar 2021 at 08:54, Sayandip Halder <sayan...@gmail.com> wrote:
>
>
>> I have recently made significant performance improvements in linsolve.
>> It would be good to make more use of linsolve in the other solvers.
>
> Do you mean this PR?

https://github.com/sympy/sympy/pull/20780

Yes, and the ones it references.

>> 4. linsolve/nonlinsolve are reasonable but not feature complete and
>> also have an awkward return type
>
>
> I am not sure what you mean by awkward return type. The FiniteSet and ImageSet returned by linsolve() and nonlinsolve() respectively
> are the same as the return type of solveset().

No, solveset returns a proper set that you can use in set operations
etc. With linsolve you get either:

a) empty set
b) a FiniteSet containing one tuple representing the unique solution
c) A FiniteSet with a tuple that has the unknowns in as free symbols

The last case c) is not a proper set:

In [1]: linsolve([x + y], [x, y])
Out[1]: {(-y, y)}

It should be something like:

In [2]: ImageSet(Lambda(y, (-y, y)), Complexes)
Out[2]: {(-y, y) │ y ∊ ℂ}

I'm not sure that users actually want an ImageSet like this but the
point is that it means the return value is not a proper set. There is
no value in returning a Set if it can not actually be used correctly
in intersections etc. With that in mind actually the FiniteSet is
quite awkward because you can't index into it like you can with the
list returned by solve. It only ever has one element so it would be
better just to return the tuple instead of wrapping it awkwardly in a
FiniteSet.

There is the same problem with nonlinsolve for positive dimensional
solution sets:

In [3]: nonlinsolve([x + y], [x, y])
Out[3]: {(-y, y)}

Again that's not a proper Set. With nonlinsolve it gets worse:

In [5]: nonlinsolve([sin(y)], [x, y])
Out[5]: {(x, {2⋅n⋅π + π │ n ∊ ℤ}), (x, {2⋅n⋅π │ n ∊ ℤ})}

This is a FiniteSet containing two Tuples of two elements one of which
is an Expr and the other is an ImageSet. There are so many things
jumbled up in there but again it is not a proper Set which would be:

In [6]: ImageSet(Lambda((x, n), (x, (2*n + 1)*pi)), Complexes, Integers)
Out[6]: {(x, π⋅(2⋅n + 1)) │ x ∊ ℂ, n ∊ ℤ}

Again I'm not sure if that is what users really want but since the
output from nonlinsolve does not represent an actual Set it can't be
used as a Set. Worse it is mixing up Expr and Set in the same Tuple
which just looks like gibberish to me. Mixing up Expr and Set also
breaks nonlinsolve's own internal logic:

In [7]: nonlinsolve([x + y, sin(y)], [x, y])
sympy/core/operations.py:64: SymPyDeprecationWarning:
Add/Mul with non-Expr args has been deprecated since SymPy 1.7. Use
Expr args instead. See https://github.com/sympy/sympy/issues/19445 for
more info.
SymPyDeprecationWarning(
---------------------------------------------------------------------------
TypeError: unsupported operand type(s) for /: 'ImageSet' and 'One'

I really think that we need to go back to the drawing board as far as
solver return types are concerned. None of them can handle
conditionally existent solutions and none of them has a satisfactory
way of handling infinite solution sets. Only solveset uses sets in a
way that is actually consistent but its return value is actually
really awkward to use.

https://github.com/sympy/sympy/issues/16861

> As for feature complete, what more features do you think should be added to them?
> #19038, #17566, #16861, PR 12424 are the most important issues and PRs that should be tackled at first, in my opinion. Do you have any more in mind that
> should be added to this list?

I think that the biggest problem is just how inefficient many of these
solve functions are. They should be much faster. Work is definitely
needed on polynomial systems as they are particularly slow.

>> None of these functions really replaces any other so we just have a
>> range of different overlapping implementations and inconstent APIs. It
>> would be good to unify as much as possible of the implementations and
>> make the APIs consistent and more useful.
>
>
> I hope to take up this as my GSOC project. Will there be any mentor for this?

There are people who could mentor this.

>
>> The most immediate thing that is needed is good documentation that
>> explains why there are different solvers and when you might want to
>> use each of them and discusses the caveats and limitations etc.
>
>
> Agreed. We should start with solve(), which is the most lacking when it comes to documentation.

Actually what I mean is that we need a page called "How to solve
equations in sympy" that explains all the possible solvers and when
you might use one or when you would use another. It should discuss in
detail what the different approaches are for and what the strengths
and weaknesses of each are.


Oscar

[Sayandip Halder]

Hello! The docs say that solve() does not return single solution, multiple solutions, no solution, etc consistently. However, solve() does return all

of them consistently, except floor() and infinite solution(which is actually handled by solveset(). Shouldn't this part of the docs be updated?

Also, what exactly is meant by "we don't know"?

---

Sayandip