Matrix inversion bug
henjo
I have been using Sympy for quite some time now in a symbolic electric
circuit solver (which I will release soon).
To the problem, I have a matrix that gives an incorrect inverse:
>>> A = Matrix([[x+y, -x, 0], [-x-y, x, 1], [0,1,0]])
>>> Ainv = A.inv()
>>> Ainv.simplify()
>>> print Ainv
⎡0 0 0⎤
⎢ ⎥
⎢0 0 1⎥
⎢ ⎥
⎣0 1 0⎦
>>> print A*Ainv
⎡ 0 0 0⎤
⎢ ⎥
⎢ 0 1 0⎥
⎢ ⎥
⎣-x - y x 1⎦
Regards,
Henrik
Ondrej Certik
Thanks for the bug report. I made it:
http://code.google.com/p/sympy/issues/detail?id=1168
In the meantime, use a different method for calculating the inverse, e.g.:
In [1]: A = Matrix([[x+y, -x, 0], [-x-y, x, 1], [0,1,0]])
In [2]: Ainv = A.inv("ADJ")
In [3]: Ainv
Out[3]:
⎡ -1 -x ⎤
⎢────── 0 ──────⎥
⎢-x - y -x - y⎥
⎢ ⎥
⎢ 0 0 1 ⎥
⎢ ⎥
⎣ 1 1 0 ⎦
In [4]: Ainv*A
Out[4]:
⎡ x y ⎤
⎢- ────── - ────── 0 0⎥
⎢ -x - y -x - y ⎥
⎢ ⎥
⎢ 0 1 0⎥
⎢ ⎥
⎣ 0 0 1⎦
In [6]: (Ainv*A).applyfunc(simplify)
Out[6]:
⎡1 0 0⎤
⎢ ⎥
⎢0 1 0⎥
⎢ ⎥
⎣0 0 1⎦
If you have time, it'd be cool if you could figure out where the
problem is exactly, see the issue. Looking forward for your symbolic
electric
circuit solver. When you release it, please send also an email to this
list, I am sure people will like to play with it. Well, to speak for
myself, I will. :)
Ondrej
henrik johansson
Thanks for the help. I dug a little deeper and found that the problem was that both the LU and GaussElimination solver
could not detect zeros in the pivot elements. The problems disappeared if the pivot elements were simplified first.
For the GE solver I could just set the simplifed argument to True in the call to the "rref" method.
But there are probably cases where a simplify is not enough. I do not know if there is a more robust way of detecting if an expression is zero.
For the circuit solver I will use the ADJ version from now on. I will keep you informed when it is ready for release.
Best regards,
Henrik
Ondrej Certik
> Hi Ondrej,
>
> Thanks for the help. I dug a little deeper and found that the problem was
> that both the LU and GaussElimination solver
> could not detect zeros in the pivot elements. The problems disappeared if
> the pivot elements were simplified first.
> For the GE solver I could just set the simplifed argument to True in the
> call to the "rref" method.
>
> But there are probably cases where a simplify is not enough. I do not know
> if there is a more robust way of detecting if an expression is zero.
Mathematically speaking, no. But in practise, I am sure it is. Try to
research some literature, and also please create a new issue for each
expression, that you cannot simplify. Also we should be able to pass a
custom function for simplifying (zero finding) to our matrix inversion
method, so that the user can pass there his own super effective zero
finder, suited for his problems. Could you please create new issues
for that?
>
> For the circuit solver I will use the ADJ version from now on. I will keep
> you informed when it is ready for release.
Thanks, see also my latest email to this list.
Ondrej
henrik johansson
> But there are probably cases where a simplify is not enough. I do not knowMathematically speaking, no. But in practise, I am sure it is. Try to
> if there is a more robust way of detecting if an expression is zero.
research some literature, and also please create a new issue for each
expression, that you cannot simplify. Also we should be able to pass a
custom function for simplifying (zero finding) to our matrix inversion
method, so that the user can pass there his own super effective zero
finder, suited for his problems. Could you please create new issues
for that?
Maybe it is better to just add a warning in the docstring to the matrix inversion and solve_linear_system functions
that the solver might have problems if the matrix not simplified. Otherwise there will be unnecessary overhead if the
expressions are already simple enough. And file all occasions where simplify cannot detect a zero as a bug.
/Henrik
Ondrej Certik
>> > But there are probably cases where a simplify is not enough. I do not
>> > know
>> > if there is a more robust way of detecting if an expression is zero.
>>
>> Mathematically speaking, no. But in practise, I am sure it is. Try to
>> research some literature, and also please create a new issue for each
>> expression, that you cannot simplify. Also we should be able to pass a
>> custom function for simplifying (zero finding) to our matrix inversion
>> method, so that the user can pass there his own super effective zero
>> finder, suited for his problems. Could you please create new issues
>> for that?
>
> Maybe it is better to just add a warning in the docstring to the matrix
> inversion and solve_linear_system functions
> that the solver might have problems if the matrix not simplified.
Indeed. Could you please prepare a patch? We'll merge it.
> Otherwise
> there will be unnecessary overhead if the
> expressions are already simple enough. And file all occasions where simplify
> cannot detect a zero as a bug.
Right. That's why we should allow the user to pass his own simplifying
functions. I think there will always be cases when simplify cannot do
it, so that when the user can pass it's own zero finder, he can fix
that easily. And let's report all failing cases as bugs to simplify.
We'll then judge on case by case basis if simplify should handle it or
not and how it should be approached.
Thanks,
Ondrej
henrik johansson
Indeed. Could you please prepare a patch? We'll merge it.
Right. That's why we should allow the user to pass his own simplifying
functions. I think there will always be cases when simplify cannot do
it, so that when the user can pass it's own zero finder, he can fix
that easily. And let's report all failing cases as bugs to simplify.
We'll then judge on case by case basis if simplify should handle it or
not and how it should be approached.
Here is a patch with an optional zero finder function argument in the inv method plus a note about the pivoting problem in the docstring. What do you think?
And I will add issues for the cases where simplify fails to detect a zero that I found so far.
Regards,
Henrik
Ondrej Certik
> Hi,
>
>>
>> Indeed. Could you please prepare a patch? We'll merge it.
>>
>> Right. That's why we should allow the user to pass his own simplifying
>> functions. I think there will always be cases when simplify cannot do
>> it, so that when the user can pass it's own zero finder, he can fix
>> that easily. And let's report all failing cases as bugs to simplify.
>> We'll then judge on case by case basis if simplify should handle it or
>> not and how it should be approached.
>
>
> Here is a patch with an optional zero finder function argument in the inv
> method plus a note about the pivoting problem in the docstring. What do you
> think?
Very nice patch, thanks! It's in:
http://hg.sympy.org/sympy/rev/5fc85c34d2ae
Welcome to our team:
http://hg.sympy.org/sympy/rev/dbb9f31ff099
> And I will add issues for the cases where simplify fails to detect a zero
> that I found so far.
Excellent, thanks a lot!
Ondrej