New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Linsolve: Linear System solver #9438

Merged
merged 12 commits into from Jun 29, 2015

Conversation

Projects
None yet
7 participants
@aktech
Copy link
Member

aktech commented May 29, 2015

PR for implementing Linear system Solver.

TODO

  • linear_eq_to_matrix method
  • Tests for linear_eq_to_matrix
  • gauss_jordan_solve
  • linsolve
  • Tests for linsolve
  • Tests for gauss_jordan_solve
  • Fix Merge Conflicts
  • 100% Coverage

@hargup @flacjacket

@jcrist jcrist added the GSoC label May 29, 2015

[ 1 2 3 1 ]
[ 3 1 1 -6 ]
[ 2 4 9 2 ]

This comment has been minimized.

@moorepants

moorepants May 29, 2015

Member

Why does it return a single matrix? I would expect it to return A and B of Ax=B.

This comment has been minimized.

@aktech

aktech May 29, 2015

Author Member

It basically returns an augmented Matrix of A & b.
It can be easily split into A and b by slicing.

This comment has been minimized.

@moorepants

moorepants May 29, 2015

Member

When would you ever use the augmented matrix [A|b] as a whole? If there is no use case for keeping it together, I'd have this return A and b separately.

This comment has been minimized.

@aktech

aktech May 29, 2015

Author Member

Since, I was using the augmented Matrix in linsolve for converting it to reduced row echleon form M.rref() , I made it return the Augmented matrix.
Though, as you said there are more use cases for the former (A, b) I will make it return A & b.


m, n = aug[:, :-1].shape

# solve by row reduced echleon form using Gauss Elimination

This comment has been minimized.

@moorepants

moorepants May 29, 2015

Member

Isn't this stuff already implemented in Matrix.LUSolve?

This comment has been minimized.

@aktech

aktech May 29, 2015

Author Member

I am not implementing it again, I am using the rref from Matrices.

This comment has been minimized.

@aktech

aktech May 29, 2015

Author Member

and also none of the solve methods in matrices, works as a General linear system solver.
Linsolve aims to implement the General Linear system solver, supporting all kinds of general forms of linear systems & input formats.

This comment has been minimized.

@moorepants

moorepants May 29, 2015

Member

Yes but isn't this all that needs to be done?:

eqs = [...]
A, b = to_matrix(eqs)
sol = A.LUSolve(b)

This comment has been minimized.

@aktech

aktech May 29, 2015

Author Member

Since, LUsolve() works for square matrices, we can't use this method for undetermined systems.

For example:

In [14]: M
Out[14]: 
Matrix([
[1, 2, 1,  1],
[1, 2, 2, -1],
[2, 4, 0,  6]])

In [15]: b
Out[15]: 
Matrix([
[ 7],
[12],
[ 4]])

In [17]: M.LUsolve(b)

Out[17]: .....
NonSquareMatrixError: A Matrix must be square to apply LUdecomposition_Simple().

This comment has been minimized.

@aktech

aktech May 29, 2015

Author Member

Though I don't have a written plan, only abstract is written in Proposal.
As solveset is the future replacement for solve .
linsolve will replace solve_linear_system method of solve.

I am not reusing solve_linear_system because it's duplicates rref method of Matrices.
So, basically I intend to support the functionality of the above function.

This comment has been minimized.

@moorepants

moorepants May 29, 2015

Member

I just read your GSoC proposal about this. Nice overview.

I have two main suggestions here that I think will help us review your work:

  1. Write all of the unit tests for linsolve before you write any implementation code (this could also include class/function prototypes if you want)
  2. Let the Matrix code be the core source for all linear solve operations (i.e. remove all algorithmic solve code from everywhere else in sympy). Then for any linear solve stuff that can't be handled by the Matrix code base, implement that here (not sure what that would be, but 1. will show us). Ideally, linsolve will only have functionality to manipulate the different styles of linear system specifications and put it into a common form (like Ax=b).

This comment has been minimized.

@aktech

aktech May 29, 2015

Author Member

Thanks for you suggestions, I will work upon it.

This comment has been minimized.

@hargup

hargup May 29, 2015

Member

I totally agree with moorepants suggestions, in additions to that you should only one canonical type to do the computations and write a wrapper to convert the other input format to the canonical type. To keep things even cleaner, you should write this wrapper as a decorator.

This comment has been minimized.

@aktech

aktech May 31, 2015

Author Member

@moorepants @hargup

I have pushed the tests for linsolve, I am working on the implementation (almost done).
In the implementation, I am only doing the pre-processing of inputs & post processing of the solutions (in linsolve function) obtained from .rref() matrix method (that's all Matrix code base can handle), for well behaved, under & over-determined systems.

###############################################################################


def eq_to_matrix(equations, *symbols):

This comment has been minimized.

@moorepants

moorepants May 29, 2015

Member

This is a very generic name. Wouldn't it be more clear if you call it linear_eq_to_matrix?

This comment has been minimized.

@aktech

aktech May 29, 2015

Author Member

Indeed. +1

This comment has been minimized.

@aktech

aktech May 29, 2015

Author Member

Though the name was inspired from Matlab.

This comment has been minimized.

@moorepants

moorepants May 29, 2015

Member

Matlab has some of the worst function names of all programming languages. Definitely don't follow their lead on that.

@aktech aktech force-pushed the aktech:linsolve branch from b4e75dd to 8adafa6 May 30, 2015

# Test for different input forms
x1, x2, x3, x4 = symbols('x1, x2, x3, x4')

M = Matrix([[1, 2, 1, 1, 7], [1, 2, 2, -1, 12], [2, 4, 0, 6, 4]])

This comment has been minimized.

@moorepants

moorepants May 31, 2015

Member

I would like to see either most of the tests with all symbols instead of integers or at least one nice fully symbolic test.

2*x1 + 4*x2 + 6*x4 - 4]

assert linsolve(M, (x1, x2, x3, x4)) == \
FiniteSet((-2*x2 - 3*x4 + 2, x2, 2*x4 + 5, x4))

This comment has been minimized.

@moorepants

moorepants May 31, 2015

Member

Do FiniteSets have a fixed order? Is that order guaranteed to match the order of (x1, x2, x3, x4)?

This comment has been minimized.

@moorepants

moorepants May 31, 2015

Member

Actually I'm confused what this is providing. Shouldn't it just provide four numbers?

This comment has been minimized.

@moorepants

moorepants May 31, 2015

Member

M doesn't even have symbols, so this doesn't make any sense to me...

This comment has been minimized.

@aktech

aktech May 31, 2015

Author Member

Actually FiniteSet is unordered, that's why we should be using FiniteSet of Ordered tuple, which is guaranteed to match the order of the variables.

This comment has been minimized.

@moorepants

moorepants May 31, 2015

Member

Ah, A is non-square so it provides a reduced form of the equations?

This comment has been minimized.

@aktech

aktech May 31, 2015

Author Member

I am almost done with linsolve implementation it would give a little more clear picture. I would push that very soon.

This comment has been minimized.

@moorepants

moorepants May 31, 2015

Member

I didn't see that A was non-square immediately that's why I was confused.

So FiniteSet is unordered? You say that you should be using an ordered tuple. It seems to me that you must preserve order somehow with respect the provided unknowns. What is the reason for returning a FiniteSet in the first place?

This comment has been minimized.

@aktech

aktech May 31, 2015

Author Member

Actually I'm confused what this is providing. Shouldn't it just provide four numbers?

FiniteSet is used to maintain consistency in the return type in the solveset code base, that's basically one of the main goal of redesigning the new solve 'solveset'.

This comment has been minimized.

@moorepants

moorepants May 31, 2015

Member

So if it is unordered, does it need to be?

This comment has been minimized.

@aktech

aktech May 31, 2015

Author Member

When we were dealing with univariate then FiniteSet was introduced in solveset as a general return type, Now that we are dealing with multivariate we should have a consensus over one return type. I think @hargup can answer this better, as it was suggested here 4th post.

@moorepants

This comment has been minimized.

Copy link
Member

moorepants commented May 31, 2015

Tests look great. Now it is very clear how you want things to work. Thanks.

@aktech aktech force-pushed the aktech:linsolve branch 2 times, most recently from 6359199 to b1f1036 May 31, 2015

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 1, 2015

@moorepants @flacjacket @hargup Please have a look.

@aktech aktech changed the title [WIP] Linsolve: Linear System solver Linsolve: Linear System solver Jun 2, 2015

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 3, 2015

@aktech aktech force-pushed the aktech:linsolve branch from b1f1036 to 177673f Jun 4, 2015

@hargup hargup referenced this pull request Jun 4, 2015

Merged

Implemented ComplexPlane Class: Complex Sets #9463

6 of 6 tasks complete
@hargup

This comment has been minimized.

Copy link
Member

hargup commented Jun 4, 2015

@aktech I've restarted the build.

@hargup

This comment has been minimized.

Copy link
Member

hargup commented Jun 4, 2015

The travis build has been completed successfully. I've verified the maths used in this code and it looks good to me and I'm +1 for a merge. @flacjacket can you add your comments. I would also like @asmeurer to look at the API and see if it is a good fit.

@moorepants

This comment has been minimized.

Copy link
Member

moorepants commented Jun 4, 2015

Two things from me:

  1. I still don't understand how an unordered set result is a good idea. How do you know which solution goes with which unknown symbol?
  2. It also isn't clear to me (probably because I haven't dug too deep) if this comment of mine was addressed completely:

Let the Matrix code be the core source for all linear solve operations (i.e. remove all algorithmic solve code from everywhere else in sympy). Then for any linear solve stuff that can't be handled by the Matrix code base, implement that here (not sure what that would be, but 1. will show us). Ideally, linsolve will only have functionality to manipulate the different styles of linear system specifications and put it into a common form (like Ax=b).

In general I would have expected more calls to high level methods in the Matrix codebase.

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 5, 2015

@moorepants

  1. I still don't understand how an unordered set result is a good idea. How do you know which solution goes with which unknown symbol?

The result here is not an unordered set, rather it is FiniteSet of Ordered tuple. See this.
The FiniteSet around the tuple is just a wrapper, so that solution looks like a set, which is consistent with rest of solveset.
If you don't like { } around the tuple we can also just return a tuple, It doesn't have much significance, except it looks like a set, (Since we don't have a SymPy ordered set).

  1. It also isn't clear to me (probably because I haven't dug too deep) if this comment of mine was addressed completely:

I am only doing the pre-processing of inputs & post processing of the solutions here obtained from .rref() matrix method that's all Matrix code base can handle, but if you think this also should be handled by Matrix code base, then all we can be done is to create a method solve_general_linear_system somewhere in matrix.py & all the code in linsolve will go into it except first few lines (of preprocessing input). Is this should I do?

@moorepants

This comment has been minimized.

Copy link
Member

moorepants commented Jun 5, 2015

The result here is not an unordered set, rather it is FiniteSet of Ordered tuple. See this.

In that link you say that FiniteSet is unordered. Is unordered or isn't it?

If you don't like { } around the tuple we can also just return a tuple, It doesn't have much significance, except it looks like a set, (Since we don't have a SymPy ordered set).

I don't care whether it is a set or a tuple. I want to ensure is that if I solve for x, y and z, I know which of the results are for x, y, and z.

2 I guess that I was imagining that once you expose A and b, then you'd examine A and b and call the appropriate method on A: A.necessary_solve_method(b) instead of doing some other processing in your function. Can you explain why it can't be this way? Or why you think it shouldn't?

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 5, 2015

Sorry for not being able to make it clear.

  • FiniteSet is unordered, So FiniteSet(1, 2, 3) = {1, 2, 3} is equivalent to FiniteSet(3, 1, 2) = {3, 1, 2}. Here the arguments to FiniteSet are integers.
  • A tuple is ordered, so a = (1, 2, 3) is not equivalent to (3, 1, 2).

Now, If I pass tuple as an argument to FiniteSet:
like arg = (1, 2, 3) toFiniteSet( arg ) which results in FiniteSet( (1, 2, 3) ) (Note the extra pair of parenthesis).
It is returned as { (1, 2, 3) }
Now this is ordered as no other arguments are passed to FiniteSet except a tuple, which is ordered!

Hence { (1, 2, 3) } is not Equivalent to { (3, 1, 2) } or { (3, 2, 1) }.
That's the reason I call the output as FiniteSet of ordered tuple, (since we just pass an ordered tuple as an argument to FiniteSet).


I guess that I was imagining that once you expose A and b, then you'd examine A and b and call the appropriate method on A: A.necessary_solve_method(b) instead of doing some other processing in your function. Can you explain why it can't be this way? Or why you think it shouldn't?

Here is an brief audit of all the methods to solve linear system contained in Matrices.py :

  • upper_triangular_solve: Matrix must be square.
  • cholesky_solve: Matrix must be square.
  • diagonal_solve: Solves only diagonal matrix
  • LDLsolve: Matrix must be square
  • LUsolve: Matrix must be square
  • QRsolve: The number of rows must be greater than columns
  • pinv_solve: Rank-deficient matrices are not yet supported

All the method defined is matrices.py are good for solving well-behaved (n x n) systems only, none of them can be used to consistently represent solutions for overdetermined and underdetermined systems, except a few extra cases solved by QR & pinv solve.

That's the prime reason linsolve is a wrapper around the reduced row echleon form method rref() .
It solves a general system by converting it into RREF form and then returns solution by determining if it's effectively well-behaved, or effectively over or under-determined system.

@moorepants

This comment has been minimized.

Copy link
Member

moorepants commented Jun 5, 2015

FiniteSet

Ok, I get it now. But I'm not sure why this is a good solution. Should you create an OrderedFiniteSet instead of doing this? It seems to me that in all of your test cases you return a set of a single tuple. If there is no case for a set of multiple tuples, I'm not sure why this is needed. It seems to unnecessarily over-complicate things. I don't know anything about the new solver's @hargup has created, so I don't really know what the reason if for returning unordered sets. Returning a dictionary that maps unknowns to their solutions seems much clearer to me, especially since it can immediately be used in methods like subs(). This design detaches the rhs of the solution equations. If you are going to have a FiniteSet of a tuple why not have a FiniteSet of a dictionary?

Using the matrix code.

I understand this now, thanks for the clear explanation. So it seems like you should use the methods you list if your A and b fit the right criteria (and you may want to offer the user the option to choose the method when you can use more than one for the given A and b matrices). Now for the over and under determined systems, would it be better to add new methods to Matrix so that you can solve those types of systems with Matrices? If not, why?

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 5, 2015

Should you create an OrderedFiniteSet instead of doing this?

That would be good, I agree, but I don't see a reason for why FiniteSet was made unordered in the first place, I would probably like to see FiniteSet being modified to be ordered.

If there is no case for a set of multiple tuples, I'm not sure why this is needed. It seems to unnecessarily over-complicate things.

I think @hargup can answer this better as it was suggested here (4th post).

Returning a dictionary that maps unknowns to their solutions seems much clearer to me, especially since it can immediately be used in methods like subs(). This design detaches the rhs of the solution equations. If you are going to have a FiniteSet of a tuple why not have a FiniteSet of a dictionary?

There has been a long discussion in the past regarding a consistent output for solvers, and the conclusion was to use sets as a common output. See this.

Now for the over and under determined systems, would it be better to add new methods to Matrix so that you can solve those types of systems with Matrices? If not, why?

The matrices code base have most of the general methods to handle linear systems, adding new method doesn't seems a good thing to me, since rather than adding new methods, we can improve the already present methods, to accept more general forms of linear systems. As there already are lots of NotImplementedError . See this.

@moorepants

This comment has been minimized.

Copy link
Member

moorepants commented Jun 5, 2015

That would be good, I agree, but I don't see a reason for why FiniteSet was made unordered in the first place, I would probably like to see FiniteSet being modified to be ordered.

That sounds fine to me.

There has been a long discussion in the past regarding a consistent output for solvers, and the conclusion was to use sets as a common output. See this.

I'm not going to read all that, sorry. Maybe someone can summarize why sets are the best choice.

adding new method doesn't seems a good thing to me

The point I'm trying to make is not about whether new methods are needed in matrices, but that all of the solve algorithm code should belong with the matrices, not inside this function you have just created.

we can improve the already present methods

Yes, that is perfectly fine.

If you are going to spend time to fix this, it'd be best to fix up the Matrix based solve code so that all of the necessary algorithms reside there. Then your function solve_linear_system is simply a very light wrapper.

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 5, 2015

If you are going to spend time to fix this, it'd be best to fix up the Matrix based solve code so that all of the necessary algorithms reside there.

OK, we can do this, by moving the algorithmic solve code to matrices.py . Since linsolve is using reduced row echleon form to solve a system, which is basically gauss jordan elimination.
We can create a method gauss_jordan_solve() and move all algorithmic code there and linsolve will simply preprocess the input and call gauss_jordan_solve() and return solutions.

@hargup

This comment has been minimized.

Copy link
Member

hargup commented Jun 5, 2015

Sets Module
The main reason for using sets as output to solvers is that it can consistently
represent many types of solutions. For single variable case it can represent:

  • No solution (by null set)
  • Finitely many solutions (by FiniteSet)
  • Infinitely many solutions, both countably and uncountably infinite solutions
    (using the ImageSet module)
  • Interval
  • There can also be bizarre solutions to equations like set of all rational number.

No other programmer's object (list, dictionary, python sets) provides the
flexibility of mathematical sets which our sets module try to emulate. The
second reason to use sets is that they are close to the entities which
mathematician's deals with and it makes it easier to reason about them.
Another advantage of using objects closer to mathematical entities is that the
user won't have to "learn" our representation and she can have her expectations
transferred from her mathematical experience.

For multivariate case we are representing solutions as a set of points in a n
dimensional space and a point is represented by a ordered tuple (maybe we should
use point from the geometry module). Just like single variable case there can be
one solution, finitely many solutions, countably infinite number of solutions
and uncountable infinite number of solutions and I feel the sets module is the
best fit for this purpose.

To give an example of two points solution a user can give a system of equations
of two circles:

>>> solveset([x**2 + y**2 -1, x**2 + (y-1/2)**2 - 1], (x, y))

Dictionary as Output

Dictionary are easy to deal with programatically but mathematically they are not
very precise and use of them can quickly lead to inconsistency and a lot of
confusion. For example addition of two solutions is a fairly valid operations
but dictionaries by themselves have no add operation. Since we are representing
solutions as dictionary we can define a custom add operation where value of x
variable is added with value of other x variable and value of y variable
with the value of other y variable. But what if the user tries to add solution
of two equations which have different sets of variables? Say one has x and y
as variables and other has u and v? We cannot allow this addition because
the keys of dictionaries are unordered but this behavior will be inconsistent
with the one variable case because it will be absurd to not to allow addition of
two numbers, say one obtained by solution of x = 1 and other by solving y = 2.

@moorepants

This comment has been minimized.

Copy link
Member

moorepants commented Jun 5, 2015

OK, we can do this, by moving the algorithmic solve code to matrices.py

That seems reasonable, but isn't LUsolve Gaussian Elimination too? It may be redundant.

Also, only do this if others agree. @hargup Do you agree that these algorithms should be in the Matrix codebase? Or should it remain here as @aktech has it.

@moorepants

This comment has been minimized.

Copy link
Member

moorepants commented Jun 5, 2015

@hargup Thank you for this explanation. Is it reasonable to have these sets be ordered? In the linsolve case I'm not seeing any example that has more than one set of tuples as a solution. So it is odd to me that you have a FiniteSet wrapped around a tuple. Maybe this looks fine in the grand scheme of solvers, but for this particular solver it just seems like an extra layer that makes it harder to get at what you want.

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 5, 2015

It may be redundant.

Why do you think so?
Gauss Jordan elimination is a different algorithm from LUsolve.

BTW, you can have a look at the last commit, I made the changes.
@hargup is fine with this, I talked to him recently (on gitter just now).

@moorepants

This comment has been minimized.

Copy link
Member

moorepants commented Jun 5, 2015

Why do you think so?

Not completely sure but from http://en.wikipedia.org/wiki/LU_decomposition, it says "The LU decomposition can be viewed as the matrix form of Gaussian elimination." Also "The LU decomposition is basically a modified form of Gaussian elimination."

@debugger22

This comment has been minimized.

Copy link
Member

debugger22 commented Jun 9, 2015

It's very less likely to work even after the restart. See this #9489

@debugger22

This comment has been minimized.

Copy link
Member

debugger22 commented Jun 9, 2015

I've restarted it, though.

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 9, 2015

@debugger22

This comment has been minimized.

Copy link
Member

debugger22 commented Jun 9, 2015

And that too with a good margin. I've no idea about that.

@aktech aktech force-pushed the aktech:linsolve branch from df9a32c to 471ff19 Jun 18, 2015

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 18, 2015

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 20, 2015

Coverage Report

  • [1] gauss_jordan_solve 100 %
  • [2] linsolve : 100: %

@hargup @flacjacket

Ref:
[1] http://iamit.in/sympy/coverage-report/matrices/
[2] http://iamit.in/sympy/coverage-report/solveset/

@aktech aktech added solvers.solveset and removed solvers labels Jun 29, 2015

@aktech aktech self-assigned this Jun 29, 2015

hargup added a commit that referenced this pull request Jun 29, 2015

Merge pull request #9438 from aktech/linsolve
Linsolve: Linear System solver

@hargup hargup merged commit 36f2be8 into sympy:master Jun 29, 2015

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
@hargup

This comment has been minimized.

Copy link
Member

hargup commented Jun 29, 2015

@aktech Thanks! It is in.

@aktech

This comment has been minimized.

Copy link
Member Author

aktech commented Jun 29, 2015

@aktech aktech referenced this pull request Jun 30, 2015

Merged

Add linsolve sphinx docs #9587

1 of 1 task complete

@hargup hargup referenced this pull request Jul 18, 2015

Closed

General linear solve #2580

1 of 5 tasks complete

@aktech aktech deleted the aktech:linsolve branch Jul 25, 2015

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment