nonlinsolve() and linsolve() fail on a under-determined system

[sgsokol]

Hi,
I am new to SymPy and trying to solve a parametric linear system with 4 variables x1, x2, x3, x4 and 4 parameters a, b, i, and j. The system is under-determined, it has 1 free variable despite the fact that there are 5 equations. However, nonlinsolve() and linsolve() does not detect that the system is under-determined and provided meaningless result. Here is a minimal reproducible example:

import sympy
x1,x2,x3,x4=sympy.symbols(", ".join("x"+str(i+1) for i in range(4)))
i,j,a,b=sympy.symbols("i, j, a, b")
sy=[x1+x2+x3+x4-1, x1+x3-i, x2+x4-j, x3-a*(x3+x4), x4-b*(x3+x4)]
sympy.nonlinsolve(sy, [x1,x2,x3,x4])
# {(i, j, 0, 0)}

sympy.linsolve(sy, [x1,x2,x3,x4])
# EmptySet()

And here is an example of a solution:

dsn={x1: -(a*x4 - b*i)/b, x2: j - x4, x3: a*x4/b, x4: x4}
# dict of substitutions which follow from the original system
dsub={(i+j):1, (a+b): 1}

# check that dsn is indeed a solution
[sympy.simplify(v.subs(dsn)).subs(dsub) for v in sy]
# [0, 0, 0, 0, 0]

I am using Python 3.7.6 and sympy 1.3 on a rpm based Linux.

[sgsokol]

I just tested with sympy 1.6.1 (the latest available version for now) and the result is the same -- no solution found.

[oscarbenjamin]

In your example you substitute values for i and j but they are not given as unknowns. The symbols in the equations that are not unknowns are parameters of the system and in this case solutions only exist for particular values of those parameters:

In [27]: A, b = linear_eq_to_matrix(sy, [x1, x2, x3, x4])                                                                         

In [28]: Matrix.hstack(A, b).echelon_form().expand()                                                                              
Out[28]: 
Matrix([
[1,  1,     1,          1,          1],
[0, -1,     0,         -1,      i - 1],
[0,  0, 1 - a,         -a,          0],
[0,  0,     0, -a - b + 1,          0],
[0,  0,     0,          0, -i - j + 1]])

From the lower-right of the echelon form of the augmented matrix we can see that a solution only exists if i+j=1 which is generically untrue for arbitrary i and j. In the case of conditionally existent solutions like this sympy's solvers will usually return the solution for the generic case (no solution in this example) rather than the degenerate case (solution that exists only for precise parameter values).

A simpler example to demonstrate what you are showing is this:

In [45]: x, y = symbols('x, y')

In [46]: solve(eqs, [x])
Out[46]: []

In [47]: linsolve(eqs, [x])
Out[47]: ∅

In [48]: nonlinsolve(eqs, [x])
Out[48]: {(-1,)}

I said that sympy usually gives the generic solution because this is not uniformly applied and we can see e.g. that nonlinsolve has returned a non-generic solution that only exists if y=1.

If instead one of the conditional parameters is given as an unknown then all solvers can find the value of that unknown:

In [50]: solve(eqs, [x, y])
Out[50]: {x: -1, y: 1}

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

In [52]: nonlinsolve(eqs, [x, y])
Out[52]: {(-1, 1)}

That will also work in the original example:

In [53]: linsolve(sy, [x1,x2,x3,x4,i])
Out[53]: {(1 - j, j, 0, 0, 1 - j)}

In [54]: linsolve(sy, [x1,x2,x3,x4,i,j])
Out[54]: {(1 - j, j, 0, 0, 1 - j, j)}

I don't think the example you have showed is incorrect. Rather I think that Out[48] above is incorrect because it returns a non-generic solution. The current solvers should return generic solutions.

We should support conditionally existent solutions as described in #16861. The solver from there fails for this example though...

[oscarbenjamin]

To be clear the system with only x1, x2, x3 and x4 as unknowns is not underdetermined. It is in fact overdetermined. There are 5 equations for 4 unknowns and those 5 equations are only consistent if i and j (not given as unknowns) take particular values. For most possible values of i and j no solution exists.

[sgsokol]

In fact, i and j are implicitly linked by the first three equations. If we combine them by adding eq2 and eq3, we get immediately i+j=1. But I admit that it could be too much to ask to sympy to take it into account by itself.

We should support conditionally existent solutions as described in #16861.

It could also make sense to add an optional parameter to linsolve() and co. say psubs which (if given) could contain substitute expression on parameters like [(i+j, 1), (a+b, 1)]. These substitutions could be applied to both lsh and rhs after gauss-jordan reduction or similar stages.

[oscarbenjamin]

You can include the equation Eq(i + j, 1) in the system of equations

In [12]: linsolve(sy + [i + j - 1], [x1, x2, x3, x4, i])                                                                                       
Out[12]: {(1 - j, j, 0, 0, 1 - j)}

In [13]: linsolve(sy + [i + j - 1], [x1, x2, x3, x4, i, j])                                                                                    
Out[13]: {(1 - j, j, 0, 0, 1 - j, j)}

In [14]: linsolve(sy + [i + j - 1], [x1, x2, x3, x4, j])                                                                                       
Out[14]: {(i, 1 - i, 0, 0, 1 - i)}

If you want linsolve to discover what the relationship between i and j should be then just make i and j unknowns without including the extra equation:

In [7]: linsolve(sy, [x1,x2,x3,x4,i])                                                                                                                         
Out[7]: {(1 - j, j, 0, 0, 1 - j)}

But I admit that it could be too much to ask to sympy to take it into account by itself.

It wouldn't be hard for sympy to do this automatically but the implication would be that linsolve returns solution sets that are not generically correct. If i and j are parameters of the system rather than unknowns then for almost all possible values of i and j the system has no solution.

That is why linsolve returns the empty set. That might not be the ideal behaviour in your usecase but it is a consistent definition.

[sgsokol]

You can include the equation Eq(i + j, 1) in the system of equations

In [12]: linsolve(sy + [i + j - 1], [x1, x2, x3, x4, i])                                                                                       
Out[12]: {(1 - j, j, 0, 0, 1 - j)}

This workaround is not as generic as suggested psubs solution. For example, we cannot do the same with condition a+b=1 as it makes appear nonlinear term a*x3. We have to use nonlinsolve() instead of staying with linsolve():

>>> sympy.nonlinsolve(sy + [i + j - 1] + [a + b - 1], [x1, x2, x3, x4, i, a])
FiniteSet((-j + x4 + 1 - x4/b, j - x4, (-b*x4 + x4)/b, x4, 1 - j, 1 - b), (-j + x4 + 1 - x4/b, j - x4, -(b*x4 - x4)/b, x4, 1 - j, 1 - b))

looks close to the searched solution.

[oscarbenjamin]

Rather than psubs it would probably make more sense to have a parameter called something like conditions that could express equalities, unequalities and inequalities as well e.g.:

linsolve(..., conditions=Eq(i + j, 1) & Eq(a + b, 1))

On the other hand at some point we probably want to be able to support an assumptions argument here so maybe that could be used as the mechanism like

linsolve(..., assumptions=Q(Eq(i + j, 1) & Eq(a + b, 1)))

For now you don't need to use nonlinsolve if you eliminate the redundant symbols:

In [6]: conditions = [Eq(i + j, 1), Eq(a + b, 1)]                                                                                                             

In [7]: rep = solve(conditions, [j, b])                                                                                                                       

In [8]: rep                                                                                                                                                   
Out[8]: {b: 1 - a, j: 1 - i}

In [9]: linsolve([s.subs(rep) for s in sy], [x1,x2,x3,x4])                                                                                                    
Out[9]: 
⎧⎛ a⋅x₄                   -a⋅x₄     ⎞⎫
⎨⎜───── + i, -i - x₄ + 1, ──────, x₄⎟⎬
⎩⎝a - 1                   a - 1     ⎠⎭