SymPy for symbolic matrix inversion
Johann
It was nice meeting you at EuroScipy in Leipzig and to chat with you
about SymPy. As we discussed, I would be interested in finding out
about possible use of SymPy for the symbolic matrix inversions we are
implementing in our PySCeS software (at the moment using Maxima).
I have uploaded 3 examples of such symbolic matrices to the "Files"
section. They are
pickled python objects, the matrix is represented as a list of lists.
Each element is currently a string. You'll see that many of them are
1 or 0, the other symbols are frequently preceded by a minus (this is
intentional, they are negative), and a symbol such as 'ecR1_S1' is a
single symbol (i.e. entity). In some cases matrix elements consist of
two terms that are being divided, e.g. 'JR1/JR2'. Here JR1 and JR2
are separate symbols for the sake of algebraic analysis.
What we are requiring is a symbolic inverse of these matrices, with
the determinant factored out of the terms (i.e. we need the
determinant separately).
There are 3 examples:
- matrix1 is a 10x10 matrix which is pretty sparse
- matrix2 is also 10x10 but has more elements
- matrix3 is 11x11 but has more terms
The analysis of the first 2 takes less than 5s using our system, the
3rd does take a couple of minutes.
This is to give you an idea of the types of matrices in our analysis.
The systems can get bigger (we've tried up to 20x20) but the
computational time required increases exponentially! These 3 examples
should be sufficient for general proof of concept. Also I don't know
whether one could use algos for sparse matrices (maybe the first one,
I don't think the 3rd one would classify as sparse?).
I'd be interested if you could take a look at these matrices and let
me know your experiences with SymPy. Also let me know if the format
of the matrices I supplied is OK for you to use, or if you need
additional info.
I look forward to hearing from you!
Best wishes
Johann Rohwer
Ondrej Certik
2008/8/7 Johann <j.m.r...@gmail.com>:
>
> Hi Ondřej and other SymPy developers
>
> It was nice meeting you at EuroScipy in Leipzig and to chat with you
> about SymPy. As we discussed, I would be interested in finding out
> about possible use of SymPy for the symbolic matrix inversions we are
> implementing in our PySCeS software (at the moment using Maxima).
>
> I have uploaded 3 examples of such symbolic matrices to the "Files"
> section. They are
> pickled python objects, the matrix is represented as a list of lists.
> Each element is currently a string. You'll see that many of them are
> 1 or 0, the other symbols are frequently preceded by a minus (this is
> intentional, they are negative), and a symbol such as 'ecR1_S1' is a
> single symbol (i.e. entity). In some cases matrix elements consist of
> two terms that are being divided, e.g. 'JR1/JR2'. Here JR1 and JR2
> are separate symbols for the sake of algebraic analysis.
>
> What we are requiring is a symbolic inverse of these matrices, with
> the determinant factored out of the terms (i.e. we need the
> determinant separately).
>
> There are 3 examples:
> - matrix1 is a 10x10 matrix which is pretty sparse
> - matrix2 is also 10x10 but has more elements
> - matrix3 is 11x11 but has more terms
>
> The analysis of the first 2 takes less than 5s using our system, the
> 3rd does take a couple of minutes.
Those numbers are using Maxima?
>
> This is to give you an idea of the types of matrices in our analysis.
> The systems can get bigger (we've tried up to 20x20) but the
> computational time required increases exponentially! These 3 examples
> should be sufficient for general proof of concept. Also I don't know
> whether one could use algos for sparse matrices (maybe the first one,
> I don't think the 3rd one would classify as sparse?).
>
> I'd be interested if you could take a look at these matrices and let
> me know your experiences with SymPy. Also let me know if the format
> of the matrices I supplied is OK for you to use, or if you need
> additional info.
Ok, I'll try to find time this afternoon and code it.
But you can try this yourself right now. Construct symbols:
In [1]: a, b, c, d = var("sym1 sym2 c something")
Construct a zero matrix:
In [2]: M = zeros((10, 10))
In the example below I'll construct a unit matrix:
In [10]: M = eye(10)
Fill it in:
In [3]: M[0, 1] = a
In [4]: M[0, 3] = b
In [5]: M[1, 5] = c
In [6]: M[2, 7] = d
print it:
In [16]: M
Out[16]:
⎡1 sym₁ 0 sym₂ 0 0 0 0 0 0⎤
⎢ ⎥
⎢0 1 0 0 0 c 0 0 0 0⎥
⎢ ⎥
⎢0 0 1 0 0 0 0 something 0 0⎥
⎢ ⎥
⎢0 0 0 1 0 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 1 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 1 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0 1 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0 0 1 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0 0 0 1 0⎥
⎢ ⎥
⎣0 0 0 0 0 0 0 0 0 1⎦
And invert it:
In [17]: M.inv()
Out[17]:
⎡1 -sym₁ 0 -sym₂ 0 c⋅sym₁ 0 0 0 0⎤
⎢ ⎥
⎢0 1 0 0 0 -c 0 0 0 0⎥
⎢ ⎥
⎢0 0 1 0 0 0 0 -something 0 0⎥
⎢ ⎥
⎢0 0 0 1 0 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 1 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 1 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0 1 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0 0 1 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0 0 0 1 0⎥
⎢ ⎥
⎣0 0 0 0 0 0 0 0 0 1⎦
And determinant:
In [19]: M.det()
Out[19]: 1
Then we can play with the inversion algorithm to get optimal results.
If I remember right, currently we have implemented 3 different ones,
as you can find out by:
In [18]: M.inv?
Since you need the determinant too, you probably want to use:
In [20]: M.inv("ADJ")
But for the example above, the "ADJ" method is slower than "GE" which
is the default. Well, we'll figure out for the matrices that you
need.
I am curious how fast it will be. And thanks a lot for providing us
the matrices, because this is a real life problem, so we'll use it as
one of our benchmarks.
How do you use maxima? Did you wrote your own wrappers? My plan is to
be able to optionally call Maxima from SymPy, so that you just code it
in SymPy and in case sympy was not fast enough for some reason, it'd
just call Maxima in the background to do the matrix inversion. As a
temporary crutch, until we make sympy faster. :)
Ondrej
Johann
> 2008/8/7 Johann <j.m.roh...@gmail.com>:
> > - matrix1 is a 10x10 matrix which is pretty sparse
> > - matrix2 is also 10x10 but has more elements
> > - matrix3 is 11x11 but has more terms
>
> > The analysis of the first 2 takes less than 5s using our system, the
> > 3rd does take a couple of minutes.
>
> Those numbers are using Maxima?
> I am curious how fast it will be. And thanks a lot for providing us
> the matrices, because this is a real life problem, so we'll use it as
> one of our benchmarks.
problems, let me know.
> How do you use maxima? Did you wrote your own wrappers?
module from Python and pipes to connect to Maxima. If you'd like more
info on the exact workings, contact Tim Akhurst (PhD student in our
group who developed the code)
takhurst at sun dot ac dot za
> My plan is to
> be able to optionally call Maxima from SymPy, so that you just code it
> in SymPy and in case sympy was not fast enough for some reason, it'd
> just call Maxima in the background to do the matrix inversion. As a
> temporary crutch, until we make sympy faster. :)
CAS, but a solution tailored to a very specific problem. The symbolic
matrix is automatically generated from a model of a cellular pathway,
this is inverted using Maxima, then the results are parsed and
symbolic expressions for the coefficients are extracted/generated. I'm
not sure in how far our interface to Maxima would be appropriate for a
generic CAS. Hasn't SAGE done this though?
Johann
Ondrej Certik
>
> On Aug 7, 11:59 am, "Ondrej Certik" <ond...@certik.cz> wrote:
>> 2008/8/7 Johann <j.m.roh...@gmail.com>:
>> > There are 3 examples:
>> > - matrix1 is a 10x10 matrix which is pretty sparse
>> > - matrix2 is also 10x10 but has more elements
>> > - matrix3 is 11x11 but has more terms
>>
>> > The analysis of the first 2 takes less than 5s using our system, the
>> > 3rd does take a couple of minutes.
>>
>> Those numbers are using Maxima?
>
> Yes.
>
>> I am curious how fast it will be. And thanks a lot for providing us
>> the matrices, because this is a real life problem, so we'll use it as
>> one of our benchmarks.
>
> OK, you're welcome. If you need larger matrices for real-life
> problems, let me know.
Thanks. We'll see after we try the small ones.
>
>
>> How do you use maxima? Did you wrote your own wrappers?
>
> Yes, we wrote our own wrappers. Basically, we're using the subprocess
> module from Python and pipes to connect to Maxima. If you'd like more
> info on the exact workings, contact Tim Akhurst (PhD student in our
> group who developed the code)
> takhurst at sun dot ac dot za
>
>> My plan is to
>> be able to optionally call Maxima from SymPy, so that you just code it
>> in SymPy and in case sympy was not fast enough for some reason, it'd
>> just call Maxima in the background to do the matrix inversion. As a
>> temporary crutch, until we make sympy faster. :)
>
> This sounds cool. Note however, that we have not implemented a generic
> CAS, but a solution tailored to a very specific problem. The symbolic
> matrix is automatically generated from a model of a cellular pathway,
> this is inverted using Maxima, then the results are parsed and
> symbolic expressions for the coefficients are extracted/generated. I'm
> not sure in how far our interface to Maxima would be appropriate for a
> generic CAS. Hasn't SAGE done this though?
Yes, I plan to use Sage wrappers for this. So that SymPy can do
basically everything that one needs if one is willing to install
maxima. And then improve.
Ondrej
Ondrej Certik
>
> Hi Ondřej and other SymPy developers
>
> It was nice meeting you at EuroScipy in Leipzig and to chat with you
> about SymPy. As we discussed, I would be interested in finding out
> about possible use of SymPy for the symbolic matrix inversions we are
> implementing in our PySCeS software (at the moment using Maxima).
>
> I have uploaded 3 examples of such symbolic matrices to the "Files"
> section. They are
> pickled python objects, the matrix is represented as a list of lists.
I download matrix1 fromt he google groups to my desktop and did:
In [7]: cPickle.load(open("/home/ondra/Desktop/matrix1"))
---------------------------------------------------------------------------
UnpicklingError Traceback (most recent call last)
/home/ondra/sympy/<ipython console> in <module>()
UnpicklingError: invalid load key, ' '.
But the matrices you sent to my private email worked. It turned out
it's extremely easy to play with this in SymPy, you just take the
object and do
Matrix(a)
and that's it. It will automatically create the correct symbols and
numbers. I was myself surprised. So here is the full script:
In [1]: import cPickle
In [2]: a = cPickle.load(open("/home/ondra/Desktop/Downloads/matrix1"))
In [3]: A=Matrix(a)
In [4]: %time f= A.inv()
CPU times: user 0.53 s, sys: 0.00 s, total: 0.53 s
Wall time: 0.53 s
In [6]: m1 = cPickle.load(open("/home/ondra/Desktop/Downloads/matrix1"))
In [7]: m2 = cPickle.load(open("/home/ondra/Desktop/Downloads/matrix2"))
In [8]: m3 = cPickle.load(open("/home/ondra/Desktop/Downloads/matrix3"))
In [9]: %time i1 = Matrix(m1).inv()
CPU times: user 0.28 s, sys: 0.00 s, total: 0.28 s
Wall time: 0.28 s
In [11]: %time i2 = Matrix(m2).inv()
CPU times: user 2.84 s, sys: 0.01 s, total: 2.85 s
Wall time: 2.88 s
In [13]: %time i3 = Matrix(m3).inv()
CPU times: user 164.29 s, sys: 0.19 s, total: 164.48 s
Wall time: 165.80 s
In [26]: %time d1 = Matrix(m1).det()
CPU times: user 1.43 s, sys: 0.00 s, total: 1.43 s
Wall time: 1.45 s
In [31]: %time d2 = Matrix(m1).det()
CPU times: user 0.20 s, sys: 0.00 s, total: 0.20 s
Wall time: 0.20 s
In [35]: %time d3 = Matrix(m3).det()
CPU times: user 111.41 s, sys: 0.01 s, total: 111.43 s
Wall time: 111.86 s
So I think that's not that bad as I feared. However now it depends on
what exactly you have and need. E.g. you said you need the determinant
factored out. Currently the inversion is in the form the gaussian
elimination returned it. E.g. try i1[0,0] and you'll see the mess. To
simplify it, just type simplify(i1[0,0]). Still a mess. However, when
you multiply with the determiant, you get something smaller:
In [42]: simplify(i1[0,0])*d1
Out[42]: -ecR1_S1⋅ecR2_S2⋅ecR3_S3⋅ecR4_S4⋅ecR5_S5⋅ecR6_S6⋅ecR7_S7⋅ecR8_S8⋅ecR9_S9
Is this correct?
Unfortunately i2[0,0] is even a bigger mess, which trim(i2[0,0]) just
hangs when simplifying it. So one should probably use some other
method when solving it.
What are you doing with the data after that? E.g. in which form do you
need it? We definitely need to employ maxima as one of our algorithms
for solving the matrix, so that we can compare. The LU method also
works fine, but ADJ is just too slow.
Ondrej
Vinzent Steinberg
On Aug 7, 5:05 pm, "Ondrej Certik" <ond...@certik.cz> wrote:
[...]
> need it? We definitely need to employ maxima as one of our algorithms
> for solving the matrix, so that we can compare. The LU method also
> works fine, but ADJ is just too slow.
for dense matrices. It's also more accurate than some other methods
(only matters when using floating point arithmetic). But it's not
really intended for symbolic purposes, so it probably produces a mess.
Vinzent
Ondrej Certik
I think so too. It's good to have it around though.
Ondrej