There are a number of open issues about proper solutions of under/over-determine…d systems of equations. I think that the fundamental problem is that we don't have a way to return conditional solution sets (from solve or solveset). I've written a linear equation solver that (ab)uses Piecewise to do this and I'll show the code below but examples first.
Consider the equation `a*x=b`. We can solve this easily:
```julia
In [2]: a, b, c, d, e = symbols('a, b, c, d, e')
In [3]: x, y, z = symbols('x, y, z')
In [4]: solve(a*x-b, x)
Out[4]:
⎡b⎤
⎢─⎥
⎣a⎦
In [5]: solveset(a*x-b, x)
Out[5]:
⎧b⎫
⎨─⎬
⎩a⎭
```
The problem with this is that these solutions are not valid for `a=0`. The solver I have written gives
```julia
In [7]: s = linsolve_cond([a*x-b], [x])
In [8]: s
Out[8]:
⎧ ℂ for a = 0 ∧ b = 0
⎪
⎪ ∅ for a = 0 ∧ b ≠ 0
⎨
⎪⎧⎛b ⎞⎫
⎪⎨⎜─,⎟⎬ for a ≠ 0
⎩⎩⎝a ⎠⎭
```
This correctly handles the three possible cases for the solution set. We can use subs with the result:
```julia
In [9]: s.subs(a, 1)
Out[9]: {(b,)}
In [10]: s.subs(a, 0)
Out[10]:
⎧ℂ for b = 0
⎨
⎩∅ otherwise
In [11]: s.subs(a, 0).subs(b, 1)
Out[11]: ∅
```
That latter example fails with the output of solve/solveset:
```julia
In [19]: solveset(a*x-b, x)
Out[19]:
⎧b⎫
⎨─⎬
⎩a⎭
In [20]: solveset(a*x-b, x).subs(a, 0)
Out[20]: {zoo⋅b}
In [21]: solveset(a*x-b, x).subs(a, 0).subs(b, 1)
Out[21]: {zoo}
In [22]: solve(a*x-b, x)[0].subs(a, 0).subs(b, 1)
Out[22]: zoo
```
The solver I have written uses split-recursion. It performs row-reduction and splits recursively each time it encounters a pivot that isn't known to be zero. This means it *should* give a correct but possibly complicated answer for any system of linear equations although I haven't tested it extensively. Here's some examples:
```julia
In [23]: linsolve_cond([a*x+y, x+y], [x, y])
Out[23]:
⎧ {(0, 0)} for a ≠ 1
⎪
⎨⎧⎛-τ₀ ⎞ ⎫
⎪⎨⎜────, τ₀⎟ | τ₀ ∊ ℂ⎬ otherwise
⎩⎩⎝ a ⎠ ⎭
In [24]: linsolve_cond([a*x+y, x+y+a-1], [x, y])
Out[24]:
⎧ {(1, 0)} for a = 0
⎪
⎪⎧⎛-τ₀ ⎞ ⎫
⎪⎨⎜────, τ₀⎟ | τ₀ ∊ ℂ⎬ for a = 1
⎨⎩⎝ a ⎠ ⎭
⎪
⎪⎧⎛-(1 - a) 1 - a⎞⎫
⎪⎨⎜─────────, ─────⎟⎬ for (a > 0 ∧ a < 1) ∨ a > 1 ∨ a < 0
⎩⎩⎝a⋅(a - 1) a - 1⎠⎭
In [25]: linsolve_cond([x+y, x+y+a], [x, y])
Out[25]:
⎧{(-τ₀, τ₀) | τ₀ ∊ ℂ} for a = 0
⎨
⎩ ∅ otherwise
In [26]: linsolve_cond([x+y, x+y+a, x-y], [x, y])
Out[26]:
⎧{(0, 0)} for a = 0
⎨
⎩ ∅ otherwise
In [27]: linsolve_cond([a*x+b*y, c*x+d*y], [x, y])
Out[27]:
⎧ 2
⎪ ℂ for a = 0 ∧ b = 0 ∧ c = 0 ∧ d = 0
⎪
⎪ {(τ₀, 0) | τ₀ ∊ ℂ} for a = 0 ∧ c = 0 ∧ (b ≠ 0 ∨ d ≠ 0)
⎪
⎪⎧⎛-d⋅τ₀ ⎞ ⎫
⎪⎨⎜──────, τ₀⎟ | τ₀ ∊ ℂ⎬ for a = 0 ∧ b = 0 ∧ c ≠ 0
⎨⎩⎝ c ⎠ ⎭
⎪
⎪ {(0, 0)} for (a = 0 ∨ a⋅d - b⋅c ≠ 0) ∧ (a ≠ 0 ∨ b ≠ 0) ∧ (a ≠ 0 ∨ c ≠ 0) ∧ (b ≠ 0 ∨ a⋅d - b⋅c ≠ 0) ∧ (c ≠ 0 ∨ a⋅d - b⋅c ≠ 0)
⎪
⎪⎧⎛-b⋅τ₀ ⎞ ⎫
⎪⎨⎜──────, τ₀⎟ | τ₀ ∊ ℂ⎬ for a⋅d - b⋅c = 0 ∧ a ≠ 0
⎪⎩⎝ a ⎠ ⎭
⎩
```
I think that SymPy really needs something like this although I'm not sure exactly how it should be integrated:
1. Should this be what linsolve (or solveset) does?
2. Should solveset generally return solutions of this form?
3. Piecewise is a subset of Expr whereas Set is a subset of Basic so I don't think Piecewise should be used here (but a similar class could be written)
4. The conditions in the output could be simpler in some cases. I'm not sure exactly how to simplify them though.
The code for the solver is here:
```python
from sympy import *
def linsolve_cond(eqs, unknowns, unique=False):
if not eqs:
return S.Complexes**len(unknowns)
# Preprocessing
A, b = linear_eq_to_matrix(eqs, unknowns)
Aaug = Matrix.hstack(A, b).tolist()
# Main workhorse:
sols_conds = _linsolve_cond(Aaug)
# sols_conds is a list of 3-tuples:
# [(solset, pivot_conds, consistency_conds),...]
#
# solset: solution set as a FiniteSet or ImageSet
# pivot_conds: list of conditions (e.g. a!=0) assumed in pivoting
# consistency_conds: list of conditions needed for existence of solutions
# Build all the separate cases into a Piecewise:
sets_conds = []
for solset, pivot_conds, consistency_conds in sols_conds:
pivot_cond = And(*pivot_conds)
consistency_cond = And(*consistency_conds)
if consistency_cond is not S.false:
sets_conds.append((solset, pivot_cond & consistency_cond))
if consistency_cond is not S.true:
sets_conds.append((S.EmptySet, pivot_cond & Not(consistency_cond)))
sets_conds_d = {}
for ss, conds in sets_conds:
if ss not in sets_conds_d:
sets_conds_d[ss] = conds
else:
sets_conds_d[ss] = Or(sets_conds_d[ss], conds)
if unique:
sets_conds_d = {s: c for s, c in sets_conds_d.items() if isinstance(s, FiniteSet)}
return Piecewise(*sets_conds_d.items())
def _linsolve_cond(Aaug, _recurse=None):
Nr, Nc = len(Aaug), len(Aaug[0])
Aorig = Matrix(Aaug)
if _recurse is None:
row, col, pivots, pivot_conds = 0, 0, [], []
else:
row, col, pivots, pivot_conds = _recurse
if pivots:
row, col = pivots[-1]
row += 1
col += 1
else:
row, col = 0, 0
sols_conds = []
# Call self recursively for alternate pivots
def recurse_zero_pivot(r, c):
pivot = Aaug[r][c]
Aaugr = [[Arc.subs(pivot, 0) for Arc in Arow] for Arow in Aaug]
pivot_condsr = pivot_conds[:] + [Eq(pivot, 0)]
_recurse = (r, c, pivots[:], pivot_condsr)
sols_conds.extend(_linsolve_cond(Aaugr, _recurse=_recurse))
while row < Nr and col < Nc-1:
# Find pivot row and swap into position
for r in range(row, Nr):
is_zero = Aaug[r][col].is_zero
if not is_zero:
if is_zero is None:
# Recurse for the case that the pivot is zero
recurse_zero_pivot(r, col)
pivot_conds.append(Ne(Aaug[r][col], 0))
if r != row:
Aaug[r], Aaug[row] = Aaug[row], Aaug[r]
break
else:
# All zeros, next column
col += 1
continue
if pivots:
assert pivots[-1][0] != row
pivots.append((row, col))
pivot_row = Aaug[row]
pivot_div = Aaug[row][col]
for r in range(row+1, Nr):
pivot_mul = Aaug[r][col]
if pivot_mul.is_zero:
continue
Aaug[r][col] = S.Zero
for c in range(col+1, Nc):
Aaug[r][c] = Aaug[r][c]*pivot_div - pivot_row[c]*pivot_mul
# Next row/column...
row += 1
col += 1
# Back substitute and list of possibilities
sol_set, consistency_conds = _back_substitute(Aaug, pivots)
sols_conds.append((sol_set, pivot_conds, consistency_conds))
return sols_conds
def _back_substitute(Aaug, pivots):
Nc = len(Aaug[0])
# Check conditions for existence of solutions then find solutions by
# back-substitution below
consistency_conds = []
for row in reversed(range(len(Aaug))):
is_zero = [e.is_zero for e in Aaug[row]]
if not all(x is True for x in is_zero[:-1]):
break
elif is_zero[-1] is False:
consistency_conds.append(S.false)
elif is_zero[-1] is None:
consistency_conds.append(Eq(Aaug[row][-1], 0))
assert (row == 0 and not pivots) or row == pivots[-1][0]
# Matrix of all zeros?
if not pivots:
solset = S.Complexes**(Nc-1)
return solset, consistency_conds
gen = numbered_symbols('tau')
params = []
sol = [None] * (Nc-1)
pivots_cols = {c:r for r, c in pivots}
for col in reversed(range(Nc-1)):
if col in pivots_cols:
r = pivots_cols[col]
lhsterms = (Aaug[r][c]*sol[c] for c in range(col+1, Nc-1))
sol[col] = (Aaug[r][-1] - Add(*lhsterms)) / Aaug[r][col]
else:
# Non-pivot gets a free symbol
sym = next(gen)
params.append(sym)
sol[col] = sym
if params:
solset = ImageSet(Lambda(tuple(params), tuple(sol)), *[S.Complexes]*len(params))
else:
solset = FiniteSet(tuple(sol))
return solset, consistency_conds
x, y, z = symbols('x, y, z')
a, b, c, d, e = symbols('a, b, c, d, e', finite=True)
unknowns = (x, y, z)
eqs = [sqrt(3)*x+y, sqrt(2)*z]
sol = linsolve_cond(eqs, unknowns)
pprint(sol)
M = Matrix(symbols('M:9')).reshape(3, 3)
xs = Matrix(symbols('x:3'))
b = Matrix(symbols('b:3'))
sol3 = linsolve_cond(list(M*xs - b), list(xs), unique=True)
```