Solving systems of polynomials in SymPy

[Maaz Muhammad]

Hi all! I'm a PhD student at the University of Toronto, working on polynomials (optimization, algorithms, and applications). 

As part of my work, I needed to solve systems of polynomials (inequalities and equalities), which SymPy currently lacks. I have therefore been working on such a solver for SymPy. The standard algorithm for this is cylindrical algebraic decomposition, which decomposes R^n into "cells" over which each polynomial is sign-invariant. You can then evaluate the system in each cell, which allows you to either get a satisfying point (or, a satisfying point from each true cell), or the algorithm returns that the system is inconsistent. Here is a demo of my code:

The implementation here returns a single satisfying point if it exists, and otherwise returns an empty list if no satisfying point exists. 

The algorithm relies heavily on the Hong projection operator. This was really the main thing I had to code. The Hong operator relies on taking resultants, reducta, and derivatives -- all of which can already be done in SymPy, so it was just a matter of putting it together (after understanding Hong's paper). 

I previously posted this on my Twitter and was asked by the SymPy Twitter account to post here. I would love to make this part of SymPy. Note that thus far, the only two implementations of CAD are in Mathematica or a software called QEPCAD. 

For now, the algorithm works with sample points from each cell (these sample points are generated during the running of the algorithm). I'm thinking that the poly_solve() function could have a parameter like 'return_type' which can take values either 'one' or 'all' -- the former would return a satisfying point from a single cell, the latter would return a satisfying point from all valid cells. 

Another consideration is that CAD can also be used to construct algebraic formulas that describe all satisfying points. This is called the solution formula. Mathematica and QEPCAD implement this. Unfortunately, this step is quite challenging, and I have not fully understood it yet. Pretty much all the information about this step is given in a PhD thesis from the 90s by Christopher Brown (who actually wrote QEPCAD). I'm studying it right now.

Overall, I have the following questions:

1. What would be the process to integrate this algorithm into SymPy?

2. Should I wait until I figure out the solution formula step to try to integrate into SymPy?

3. Are there any rules about me eventually writing a paper about my implementation in, for example, the INFORMS Journal of Computing (https://pubsonline.informs.org/page/ijoc/editorial-statement#Software%20Tools)? 

Thanks!

[Aaron Meurer]

On Tue, Jul 2, 2024 at 4:22 PM Maaz Muhammad <mmaa...@gmail.com> wrote:

Hi all! I'm a PhD student at the University of Toronto, working on polynomials (optimization, algorithms, and applications). 

As part of my work, I needed to solve systems of polynomials (inequalities and equalities), which SymPy currently lacks. I have therefore been working on such a solver for SymPy. The standard algorithm for this is cylindrical algebraic decomposition, which decomposes R^n into "cells" over which each polynomial is sign-invariant. You can then evaluate the system in each cell, which allows you to either get a satisfying point (or, a satisfying point from each true cell), or the algorithm returns that the system is inconsistent. Here is a demo of my code:

The implementation here returns a single satisfying point if it exists, and otherwise returns an empty list if no satisfying point exists. 

The algorithm relies heavily on the Hong projection operator. This was really the main thing I had to code. The Hong operator relies on taking resultants, reducta, and derivatives -- all of which can already be done in SymPy, so it was just a matter of putting it together (after understanding Hong's paper). 

I previously posted this on my Twitter and was asked by the SymPy Twitter account to post here. I would love to make this part of SymPy. Note that thus far, the only two implementations of CAD are in Mathematica or a software called QEPCAD. 

Yes, that was me who talked to you on Twitter. SymPy has needed CAD for a while (you can find several references to it if you search the issue tracker and mailing list archives), so it's good to finally see an implementation.  

For now, the algorithm works with sample points from each cell (these sample points are generated during the running of the algorithm). I'm thinking that the poly_solve() function could have a parameter like 'return_type' which can take values either 'one' or 'all' -- the former would return a satisfying point from a single cell, the latter would return a satisfying point from all valid cells. 

Another consideration is that CAD can also be used to construct algebraic formulas that describe all satisfying points. This is called the solution formula. Mathematica and QEPCAD implement this. Unfortunately, this step is quite challenging, and I have not fully understood it yet. Pretty much all the information about this step is given in a PhD thesis from the 90s by Christopher Brown (who actually wrote QEPCAD). I'm studying it right now.

Overall, I have the following questions:

1. What would be the process to integrate this algorithm into SymPy?

The best place to start is to make a pull request on GitHub. No doubt there will be several things about your implementation that will need to be cleaned up to fit within the existing framework of SymPy, but there is no way to know about that until we see the code, so a pull request is the best place to start. For example, I don't know what sorts of tests you've written, but you will need to have extensive tests for your code so that we can be sure it is correct. And we'll need to discuss what all the public APIs for this should look like. 

Let us know if you are new to git/GitHub and we can help you with that.

 

2. Should I wait until I figure out the solution formula step to try to integrate into SymPy?

No, if what you have already works, you should add it. It's better to make incremental progress than to try to do everything all at once. I'm sure there will already be enough things that need to be done in your existing code without also having to worry about the more complicated things. We can always add new functionality later.

3. Are there any rules about me eventually writing a paper about my implementation in, for example, the INFORMS Journal of Computing (https://pubsonline.informs.org/page/ijoc/editorial-statement#Software%20Tools)? 

We don't have any rules about that. The journal might have rules. You'd have to check with them. I expect the main thing is that if a SymPy contributor ends up helping significantly with the implementation they may ask to be added as a co-author.

Aaron Meurer

 

Thanks!

--
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/d6688cc4-1e09-49fe-800e-a16e3f9e8ba4n%40googlegroups.com.

[Oscar Benjamin]

Hi Maaz,

This sounds like an excellent addition to SymPy.

There was also a recent GitHub issue that discussed a more limited form of CAD:

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

The implementations there are based on this paper by Strzebonski:

https://core.ac.uk/download/pdf/82649664.pdf

There are many different things that a CAD (or GCAD) would be useful for so there are likely a number of different interfaces that it would be good to provide.

If you send your code as a pull request then it will be easier to discuss it. Alternatively if the method is similar to the one discussed in gh-26177 then you could discuss it there or open a new issue to discuss it.

One thing that SymPy currently lacks is a way to get the roots of a polynomial with irrational algebraic coefficients in RootOf form:

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

In general that is needed to be able to represent points in the non-full-dimensional cells.

--

Oscar