It seems that we have two different rref (reduced row echelon form) algorithms, one in Matrix.rref and one in solve_linear_system in solvers.py. I don't know too much about these algorithms, so I don't know if they really are the same algorithm, but in either case, I think solve should just be calling Matrix.rref.

This might speed things up too (examples taken from Lay's Linear Algebra, Third Edition):

In [15]: c = Matrix([[1, -2, 1, 0],[0, 2, -8, 8],[-4, 5, 9, -9]])

In [16]: c Out[16]: ⎡1 -2 1 0 ⎤ ⎢ ⎥ ⎢0 2 -8 8 ⎥ ⎢ ⎥ ⎣-4 5 9 -9⎦

In [17]: c.rref() Out[17]: ⎛⎡1 0 0 29⎤, [0, 1, 2]⎞ ⎜⎢ ⎥ ⎟ ⎜⎢0 1 0 16⎥ ⎟ ⎜⎢ ⎥ ⎟ ⎝⎣0 0 1 3 ⎦ ⎠

In [19]: solve_linear_system(c, x, y, z) Out[19]: {x: 29, y: 16, z: 3}

In [25]: %timeit solve_linear_system(c, x, y, z) 100 loops, best of 3: 11.2 ms per loop

In [26]: %timeit c.rref() 1000 loops, best of 3: 3.59 ms per loop

In [30]: e = Matrix([[4, -1, -1, 0, 0, 0, 0, 0, 5],[-1, 4, 0, -1, 0, 0, 0, 0, 15],[-1, 0, 4, -1, -1, 0, 0, 0, 0],[0, -1, -1, 4, 0, -1, 0, 0, 10],[0, 0, -1, 0, 4, -1, -1, 0, 0],[0, 0, 0, -1, -1, 4, 0, -1, 10],[0, 0, 0, 0, -1, 0, 4, -1, 20],[0, 0, 0, 0, 0, -1, -1, 4, 30]])

In [32]: e.rref() Out[32]: <output omitted>

In [33]: var('x1 x2 x3 x4 x5 x6 x6 x7 x8') Out[33]: (x₁, x₂, x₃, x₄, x₅, x₆, x₆, x₇, x₈)

In [34]: solve_linear_system(e, x1, x2, x3, x4, x5, x6, x7, x8) Out[34]: ⎧ 827 1377 886 1546 1171 1831 1967 2517⎫ ⎨x₁: ───, x₂: ────, x₃: ───, x₄: ────, x₅: ────, x₆: ────, x₇: ────, x₈: ────⎬ ⎩ 209 209 209 209 209 209 209 209 ⎭

In [35]: %timeit e.rref() 10 loops, best of 3: 65 ms per loop

In [36]: %timeit solve_linear_system(e, x1, x2, x3, x4, x5, x6, x7, x8) 10 loops, best of 3: 84.5 ms per loop

Sorry, you need to run

from sympy.solvers.solvers import solve_linear_system

first for the above to work.

One thing that will need to be added to Matrix.rref is the ability to exit early if the system is inconsistent, as solve_linear_system does, so that we do not waste time reducing an inconsistent matrix when don't care.

I'm looking at the two functions and, as I understand them so far, they are using almost identical algorithms. The largest difference that I see is that solve_linear_system has the early exits (as it should) while rref continues to reduce as much as it can. I can try to fold these into a single rref call in the matrix module.

solve_linear_system() is not user-friendly, and has un-conventional behavior. It and the related functions for solving linear equations that should be replaced.

In many elementary algebraic settings, one can obtain a system of linear equations that needs solving. It would be great if we could just call solve_linear_system(equations, variables) and get our solution. However, solve_linear_system() assumes our system is in an implicit matrix form. We can obtain solutions using solve(equations,variables) and solve_poly_system(equations,variables). However, this relies on both functions being smart enough to recognize the special linear properties of the system, which we are often already aware of. We could also write code to transform our system to matrix form for solve_linear_system(), but seems like much un-necessary work.

It seems natural to me that solve_linear_system() or some equivalent function should be able to recognize and pre-process all natural forms of linear systems.

A related issue is that the behaviors of solve_linear_system and solve_poly_system are inconsistent.

In [25]: from sympy import *
In [26]: from sympy.abc import x,y
In [27]: m = Matrix([[1,2,0],[3,-4,1]])
In [28]: solve_linear_system(m,x,y)
Out[28]: {x: 1/5, y: -1/10}

In [29]: solve_linear_system(m,[x,y])
Out[29]: []

In [32]: e = [x+2*y,3*x-4*y-1]

In [33]: solve_poly_system(e,x,y)
Out[33]: [(1/5, -1/10)]

In [34]: solve_poly_system(e,[x,y])
Out[34]: [(1/5, -1/10)]

Both approaches have merit, and return the correct answer. solve_linear_system's returning of a dictionary makes potential assignments explicit. However, it fails when the variables are given inside an iterable. The output of solve_poly_system() is consistent. And solve_poly_system also tolerates different input formats.

We further note that solve_linear_system() assumes that our matrix is in the form of an augmented matrix, where the last column is on the right of the equals sign when writen in matrix-vector form. However, there is no way to tell that this is an augmented matrix rather than a regular matrix, just from inspection of the input. Users might naturally assume that solve_linear_system calculates the nullspace of the given matrix, which would return the opposite sign of the correct answer. In fact, solve_linear_system() can be replaced with a 1-liner transforming the matrix to reduced row-eschelon form.

In [35]: dict(zip((x,y),list(m.rref()[0][:,-1])))
Out[35]: {x: 1/5, y: -1/10}

There is a 3rd function, solve_linear_system_LU,

In [37]: solve_linear_system_LU(m,[x,y])
Out[37]: {x: 1/5, y: -1/10}

solve_linear_system_LU() has it's own idiosyncratic behavior, as it only accepts iterables of variables.

In [38]: solve_linear_system_LU(m,x,y)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/home/junk/src/sympy/bin/<ipython-input-38-4683da0e9657> in <module>()
----> 1 solve_linear_system_LU(m,x,y)

TypeError: solve_linear_system_LU() takes exactly 2 arguments (3 given)

The 4th function, solve_linear(), seems to be intended for use in the solution of systems of polynomial equations. However, it also has non-generic behavior, sometimes returning numerator/denominator of the argument, and sometimes returning a solution of a linear equation. The exact reason for this behavior remains unclear to me, and I'd appreciate anybody who can help me out in understanding this. To me, the name atleast appears misleading.

Currently, all 4 functions are used sparesly within sympy's code base. This is a surprise, in itself, given the frequency with which linear systems are employed in mathematics, and suggests that they are not fulfilling the niches they are named for. A consequence is that a project-wide change would be a managable undertaking at this point in time. The potential impact on existing code bases outside the project is unclear.

There are several potential options.

1) Standardize on one particular form of argument style for solvers, and reimplement. Deprecate obsolete versions.

2) Add code so that the existing solvers allow variety of argument styles.

3) Leave as-is. Pros: no changes needed. Cons: New users may be confused by different styles, and may be less likely to adopt for regular usage. Users who do adopt will be channeled toward "solve()" which may be in-efficient in many simple cases. As usage spreads, design changes breaking backward compatibility will get more difficult.

The code below is my proposal for a replacement, and I am looking for feedback. The function is just a wrapper around rref(), which does all the hard work.

from sympy import *

def solve_lin_sys(eqs,xs): """ eqs is an iterable of linear equations maybe Eq(,) or expr if expr, we assume it is implicitly set to zero. xs is an iterable of variables we are solving for """ if isinstance(xs,Symbol): xs = (xs,) else: xs = tuple(xs,)

if isinstance(eqs,Matrix):
    assert eqs.shape[1] <= len(xs)+1, "too many matrix columns for variables"
    assert eqs.shape[1] >= len(xs)+1, "too few matrix columns for variables"

else:
    # move all equation arguments to one side
    if not hasattr(eqs, "__iter__"):
        eqs = ( eqs, )
    new_eqs = []
    for e in eqs:
        if isinstance(e,Equality):
            e = e.lhs-e.rhs
        for x in xs:
            assert diff(e,x,2) == 0, "not a linear system"
        new_eqs.append(e)
    eqs = tuple(new_eqs)

    # transform from equations to matrix form
    S = [ (x,0)  for x in xs ]
    M = zeros(len(eqs), len(xs)+1)
    for j,e_j in enumerate(eqs):
        for i,x_i in enumerate(xs):
            M[j,i] = diff(e_j,x_i)
        M[j,-1] = -1*e_j.subs(S)
    eqs = M

# solve by row-reduction
eschelon, pivots = eqs.rref()

# construct the returnable form of the solutions
p = len(pivots)
if p > len(xs):
    return None
sols = eschelon[0:p,p:]*Matrix([ [x] for x in xs[p:] ]+[ [1] ])
return dict(zip(xs,sols))

def example(): from sympy.abc import x,y,z eqs = [x + y - 3, x - y - 4] vs = [x,y]

# works for a list of expressions and a list of variables
print solve_lin_sys(eqs,vs)

# works for a list of equations
eqs = ( Eq(x + y, 3), Eq(x, y + 4) )
print solve_lin_sys(eqs,vs)

# works for a matrix and a list of variables
eqs = Matrix([[1,1,3],[1,-1,4]])
print solve_lin_sys(eqs,vs)

# works for simple equations
print solve_lin_sys(x+y**2 + 3*x,x)

# for over-determined systems
eqs = Matrix([[1,1,3],[1,-1,4],[2,-1,-5],[4,4,-2]])
print solve_lin_sys(eqs,vs)

# for under-determined systems, leave the last variables free
eqs = Matrix([[1,1,0,3],[1,-1,1,4]])
print solve_lin_sys(eqs,[x,y,z])

example()

Thanks for writing this little manifesto. You clearly care about the problem!

I for the most part agree with what you're saying here. This is a problem. I actually don't understand why we even need all these separate solve functions at the user level. Shouldn't solve() be able to figure out what is meant from the input? Matching an equation or system of equations to linear or polynomial or whatever is fast and easy. And if the user wants something different, we could employ a hints system. Actually, according to my original comment, solve_linear_system didn't use to be part of the user namespace. But for some reason, it was added.

Regarding changing the API, we will have to deprecate the old behavior, if possible. This will make things easiest on our users.

Regarding uses in SymPy, I personally just use Matrix directly. One very important (perhaps the most important) use of solve_linear_system in SymPy is the one in heurisch. See issue 1441, in particular comment 5. So that can provide a useful benchmark for the linear system stuff. That's the old Risch. The new Risch (which I wrote) also has to solve some linear systems, but as I said, it just uses Matrix and rref directly. Even so, the ability to exit early when we know the system is inconsistant would be nice.

Sorry for the length -- it grew over a few days.

We could go with just having solve() do everything. But solve() has to introduce overhead to classify systems, and as we see in dsolve(), this overhead can be significant. Often, we already know the system we're solving is linear, so it seems natural to have such a function.

A second reason would be that there is a very natural way to return solutions for under-determined linear systems, but not in general. Having a linear system solver can avoid the need for post-processing. (I'm particularly thinking of the calculation of series solutions of linear ODE's) Your need to exit early for over-determined systems is a similar benefit.

But these are not strong arguments, and I agree that dropping the current functions in favor of solve + hints would be a good improvement.

And now there is also the solve_lin_sys in the polys.

Cool extension being able to specify the ring. When solving underdetermined systems, you can specify the free variables by re-ordering the variables in ring.

Can we extend this by allowing ring to be an iterator or sequence of unknowns also? (a natural way would be to generate a ring based on the given sequence, but I don't see how to get ring() to do this)

If this gets fixed, then perhaps we start depricating the other methods. To summarize current state, to solve linear systems, we have

solve(equations,unknowns)
solve_linear_system(aug_matrix)
solve_linear_system_LU(aug_matrix,unknowns)
solve_poly_system(equations,unknowns)
solve_linear(single_equation)
polys.solvers.solve_lin_sys(equations,ring)

and

dict(zip(unknowns,list(aug_matrix.rref()[0][:,-1])))

We have moved issues to GitHub https://github.com/sympy/sympy/issues.

Migrated to http://github.com/sympy/sympy/issues/4949