solveset incorrectly returns EmptySet, should return ConditionSet

[agryman]

While converting some nuclear physics code from Maple to SymPy, I ran into a situation where solveset returned EmptySet when there was in fact a solution. The expression depended on a parameter B. The behaviour of solveset changed for B=14 versus B=15.

I confirmed that solutions exist in both cases using nsolve.

solveset found the solution for B=14 but not for B=15 in which case it returned EmptySet. According to the docs for solveset, it should return ConditionSet when it can't find a solution.

I also used simplify on the expression, but this actually made matters worse. solveset returned EmptySet in both cases.

I've attached a script (debug_rwc_alam.py.txt ) to reproduce the bug. Here is the output it generated:

Python 3.10.0 (default, Oct 13 2021, 06:44:31) [Clang 12.0.0 (clang-1200.0.32.29)]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.0.1 -- An enhanced Interactive Python. Type '?' for help.
PyDev console: using IPython 8.0.1
Python 3.10.0 (default, Oct 13 2021, 06:44:31) [Clang 12.0.0 (clang-1200.0.32.29)] on darwin
runfile('/Users/arthurryman/Documents/repositories/agryman/acmpy/src/acmpy/debug/debug_rwc_alam.py', wdir='/Users/arthurryman/Documents/repositories/agryman/acmpy/src/acmpy/debug')
==========Solving equation for B = 14 ==========
Success: nsolve solution = 19.676355436230185
Success: first solveset solution = 19.676355436230185
Error: solveset returns EmptySet after simplification
==========Solving equation for B = 15 ==========
Success: nsolve solution = 20.911495067823516
Error: solveset returns EmptySet
Error: solveset returns EmptySet after simplification

I'm using the latest SymPy from Pip. I've attached my requirements.txt file
debug_rwc_alam.py.txt
requirements.txt
.

[eagleoflqj]

Note that 3/2 results in a float 1.5, which may lead to further inaccuracy and solver failure.
If you use Rational(3, 2) or S(3)/2 your example will work fine.

==========Solving equation for B = 14 ==========
Success: nsolve solution = 19.676355436230185
Success: first solveset solution = 19.676355436230185
Success: first solveset solution after simplification = 19.676355436230185
==========Solving equation for B = 15 ==========
Success: nsolve solution = 20.911495067823516
Success: first solveset solution = 20.911495067823516
Success: first solveset solution after simplification = 20.911495067823516

[agryman]

@eagleoflqj Thanks for the suggestion. I appreciate it. I verified it worked for the example. However, the example was a simplified version of the actual code which takes some float parameters. In the example that 3 was actually the value of a float input c1.

The nuclear physics problem behind the code is to find the parameter value that minimizes the expectation value of a Hamiltonian operator. This is accomplished by finding the zeros of the first derivative. However, there could be more than one zero, so the code picks the one that gives the minimum expectation value. So solveset is preferable to nsolve in this case.

As I mentioned above, the current Python implementation is a very direct translation from some Maple code, so that why it uses solveset. However, given that this is a numerical task, it might be a better approach to use the SciPy minimize function.

That being said, I'd still think solveset should be fixed so that it returns ConditionSet when it fails to solve the equations.

[oscarbenjamin]

This is the code:

from sympy import *

A = Symbol('A', real=True)
B = symbols('B')
mu = sqrt(9 + (A * 3 / 2) ** 2)
E = (-3 / 2) ** 2 * (-9 * A ** 5 / mu ** 2
                            - A ** 3 * B ** 2 * 3
                            + A ** 2 * B ** 2 * 2 * (mu + 3)) \
          + A ** 3 * (2 * mu + 9) \
          - B ** 2 * mu * (mu + 2) * (-A * 3 + 2 * (mu + 4))

E15 = E.subs(B, 15)
print(solveset(E15, A, Reals))

The correct solution is found internally by solveset but then rejected here:

sympy/sympy/solvers/solveset.py

Lines 883 to 890 in 918eb81

	f_set = [] # solutions for FiniteSet
	c_set = [] # solutions for ConditionSet
	for s in result:
	if checksol(f, symbol, s):
	f_set.append(s)
	else:
	c_set.append(s)
	solution_set = FiniteSet(*f_set) + ConditionSet(symbol, Eq(f, 0), FiniteSet(*c_set))

Firstly checksol fails because substituting the float value into the equation give 1e-9 which I guess crosses the threshold for checksol:

(Pdb) p f
-20.25*_R**5/(9*_R**2/4 + 9) + _R**3*(2*sqrt(9*_R**2/4 + 9) + 9) - 1518.75*_R**3 + 1012.5*_R**2*(sqrt(9*_R**2/4 + 9) + 3) - 225*sqrt(9*_R**2/4 + 9)*(sqrt(9*_R**2/4 + 9) + 2)*(-3*_R + 2*sqrt(9*_R**2/4 + 9) + 8)
(Pdb) p s
20.9114950678235
(Pdb) p f.subs(symbol, s)
2.79396772384644e-9
(Pdb) p checksol(eq, symbol, s)
False

It should then go into the ConditionSet but that evaluates to the empty set:

(Pdb) p ConditionSet(symbol, Eq(f, 0), FiniteSet(*c_set))
EmptySet

That happens because of the check here:

sympy/sympy/sets/conditionset.py

Lines 140 to 147 in 918eb81

	if isinstance(base_set, FiniteSet):
	sifted = sift(
	base_set, lambda _: fuzzy_bool(condition.subs(sym, _)))
	if sifted[None]:
	know = FiniteSet(*sifted[True])
	base_set = FiniteSet(*sifted[None])
	else:
	return FiniteSet(*sifted[True])

Again that's because substituting the value into an Eq gives False:

(Pdb) p Eq(f, 0).subs(symbol, s)
False

I wonder if maybe evaluating an Eq with floats to false is something that should be handled more carefully. Certainly the checking in solve_radical could be smarter. If we have a good approximation of the root then we can make it more accurate:

In [27]: sol = 20.9114950678235

In [28]: sol_high_prec = nsolve(E15, A, sol, prec=50)

In [29]: sol_high_prec
Out[29]: 20.911495067823516743451965389930231658709938885148

In [30]: checksol(E15, A, sol)
Out[30]: False

In [31]: checksol(E15, A, sol_high_prec)
Out[31]: True

In [32]: E15.subs(A, sol)
Out[32]: -1.16415321826935e-9

In [33]: E15.subs(A, sol_high_prec)
Out[33]: -5.0446744715693414553254264998436980726089429907555e-44

If checksol used interval arithmetic then it could be more robust at verifying zeros.

[agryman]

@oscarbenjamin Thanks for digging into this. In my application I really don't need a very high-precision result. Perhaps solveset should allow the user to specify a tolerance along the lines of the math.isclose() function?

[oscarbenjamin]

You might not need a high precision result but checksol does. Because your original equation has radicals solveset will try to unradicalise it to get a polynomial equation:

In [12]: from sympy.solvers.solvers import unrad

In [13]: E15
Out[13]: 
                ⎛       __________    ⎞                          ⎛     __________    ⎞            _
         5      ⎜      ╱    2         ⎟                          ⎜    ╱    2         ⎟           ╱ 
  20.25⋅A     3 ⎜     ╱  9⋅A          ⎟            3           2 ⎜   ╱  9⋅A          ⎟          ╱  
- ──────── + A ⋅⎜2⋅  ╱   ──── + 9  + 9⎟ - 1518.75⋅A  + 1012.5⋅A ⋅⎜  ╱   ──── + 9  + 3⎟ - 225⋅  ╱   
     2          ⎝  ╲╱     4           ⎠                          ⎝╲╱     4           ⎠       ╲╱    
  9⋅A                                                                                              
  ──── + 9                                                                                         
   4                                                                                               

_________ ⎛     __________    ⎞ ⎛              __________    ⎞
   2      ⎜    ╱    2         ⎟ ⎜             ╱    2         ⎟
9⋅A       ⎜   ╱  9⋅A          ⎟ ⎜            ╱  9⋅A          ⎟
──── + 9 ⋅⎜  ╱   ──── + 9  + 2⎟⋅⎜-3⋅A + 2⋅  ╱   ──── + 9  + 8⎟
 4        ⎝╲╱     4           ⎠ ⎝         ╲╱     4           ⎠
                                                              
                                                              
                                                              

In [14]: unrad(E15)
Out[14]: 
⎛         12             10              9               8                7                 6      
⎝- 144.0⋅A   - 196128.0⋅A   + 1101600.0⋅A  + 79672788.0⋅A  + 162810000.0⋅A  + 1236978720.0⋅A  - 513

          5                 4                  3                  2                                
993600.0⋅A  + 7013952000.0⋅A  - 11588832000.0⋅A  + 17496000000.0⋅A  - 27993600000.0⋅A + 16329600000

      ⎞
.0, []⎠

The polynomial equation can be solved symbolically and its roots give the solutions of the original radical equation except that unradicalisation introduces spurious roots. For example these are the real roots of the polynomial:

In [18]: roots(unrad(E15)[0], filter='R')
Out[18]: {-15.7204549068037: 1, 20.9114950678235: 1}

Only one of those is a solution of the original equation. Here checksol is used to figure out which roots of the polynomial are genuine solutions of the original equation. Evaluating the equation at positive (correct) root gives a larger than expected residual which checksol rejects:

In [19]: E15.subs(A, 20.9114950678235)
Out[19]: -1.16415321826935e-9

In [20]: checksol(E15, A,-1.16415321826935e-9)
Out[20]: False

Then checksol decides that this isn't a valid root.

There are several ways that this can be improved:

1. Instead of using floats we can use a RootOf type object that can be evaluated to higher precision. For polynomial roots we already have RootOf (although it could be made faster) but it is also possible to make similar objects for non-polynomial equations by using the interval Newton method.
2. Rather than checking if something is a solution of the equation we can check if the signs match up i.e. does the sign under the radical match what would be needed. The issue here is just about which branch of the radical is chosen so we just need to verify that we have the correct branch and then we know that it is a solution.

[agryman]

@oscarbenjamin Is it correct to say that checksol only needs enough precision to be able to reject spurious solutions introduced by unradicalization?