[sympy] Re: Symbolic matrix inversion
On Sep 29, 9:49 pm, Luke wrote: > I was thinking that it might be nice to have pre-computed matrix inverses > for n x n matrices. Matrix inversion is O(n^3), so it would be nice to have > all this precomputed symbolically, and this would greatly speed up Sympy's > matrix capabilities. Inverses up to say 100x100 could be computed (or maybe > something smaller), and then when you need an inverse, everything would be > fast. This could also be used behind the scenes (by introduction of > symbolic substitution dictionaries) for inverting a matrix full of sympy > expressions. I think there won't be a huge speed advantage, because in the end you get a lazily evaluated version of the result of the inversion algorithm. If nothing simplifies somehow (and I think it doesn't for the general dense case), you only avoid some sympy overhead, if at all. Vinzent --~--~-~--~~~---~--~~ You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to sympy@googlegroups.com To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sympy?hl=en -~--~~~~--~~--~--~---
[sympy] Re: Symbolic matrix inversion
On Tue, Sep 29, 2009 at 4:03 PM, Luke wrote:
>
>
>
> On Sep 29, 1:09 pm, Ondrej Certik wrote:
>> On Tue, Sep 29, 2009 at 12:49 PM, Luke wrote:
>> > I'm using Sympy from within PyDy to generate the equations of motion for
>> > mechanical systems. At the end of the day, the equations can be most
>> > generally written as:
>> > M(x) * x'' = F(x, x', t)
>> > M(x) is what is known as the mass matrix, and will in general depend on the
>> > configuration of the system (positions and angles). This matrix needs to
>> > be
>> > inverted in order to solve for x'', which then will allow for numerical
>> > integration or stability analysis.
>> > I am generating M(x) symbolically, and in some case is it sparse, but in
>> > many cases, it isn't. Each entry of the matrix is a Sympy expression, but
>> > instead of trying to invert a matrix of Sympy expressions, I introduce a
>> > dummy symbol for each non-zero entries and then invert the matrix of dummy
>> > symbols. Humble systems might have 5 degrees of freedom (so a 5x5 needs to
>> > be inverted), so this inversion isn't so bad, but beyond this, the matrix
>> > inversions take a really long time, especially when the matrices are full
>> > (no zero entries).
>> > I was thinking that it might be nice to have pre-computed matrix inverses
>> > for n x n matrices. Matrix inversion is O(n^3), so it would be nice to
>> > have
>> > all this precomputed symbolically, and this would greatly speed up Sympy's
>> > matrix capabilities. Inverses up to say 100x100 could be computed (or
>> > maybe
>> > something smaller), and then when you need an inverse, everything would be
>> > fast. This could also be used behind the scenes (by introduction of
>> > symbolic substitution dictionaries) for inverting a matrix full of sympy
>> > expressions.
>>
>> The inversion of a 100x100 dense matrix would be quite a big mess, wouldn't
>> it?
>>
>
> Yes, it would be disaster.
>
>> But I see what you mean, and I think it makes sense to cache it, if it
>> speeds things up. But see below.
>>
>> > Does anybody know if this has been done by somebody somewhere, or have any
>> > other ideas on how it could be done better than the way I suggested?
>>
>> I would first try to profile the inversion code to see why it is slow.
>> Because for example the adjoint method is just a ratio of two
>> determinants, so it may be the determinants calculation that is slow
>> (and should be cached). However, there are 100^2 different
>> determinants (right?), but I guess caching 1 expressions should
>> still be doable. But if this is the case, we should just speed up the
>> determinant calculation, imho. Looking at the code, we have 2
>> algorithms implemented, bareis and berkowitz.
>>
>
> I did some timings of the three matrix inversion (GE, ADJ, LU) using
> the timeit module. I also timed the two determinant methods as well
> to see they the two perform side by side. It seems that bareis (the
> default one) is drastically slower than berkowitz, even when you
> call .expand() on the berkowitz determinant to make it end up with
> identical symbolic expressions. Maybe this should be the default
> method for .det()? I went up to 5x5 matrices, it started to get
> really slow after that. Gaussian elimination seems to be the slowest
> one, at least for dense matrices like the ones I used. Here is the
> code I used to do the timings:
>
> import timeit
> from numpy import zeros, max
> import matplotlib.pyplot as plt
>
> # Dimension of matrix to invert
> n = range(2, 6)
> # Number of times to invert
> number = 20
> # Store the results
> t = zeros((len(n), 5))
>
> for i in n:
> setup_code = "from sympy import Matrix, Symbol\nM = Matrix(%d,\
> %d"%(i,i)+",lambda i,j: Symbol('m%d%d'%(i,j)))"
> t[i-2, 0] = timeit.Timer('M.inverse_GE()', setup_code).timeit
> (number)
> t[i-2, 1] = timeit.Timer('M.inverse_ADJ()', setup_code).timeit
> (number)
> t[i-2, 2] = timeit.Timer('M.inverse_LU()', setup_code).timeit
> (number)
> t[i-2, 3] = timeit.Timer('M.det(method="bareis")',
> setup_code).timeit(number)
> t[i-2, 4] = timeit.Timer('M.det(method="berkowitz").expand()',
> setup_code).timeit(number)
>
>
> plt.plot(n, t[:,0]/number, label='GE')
> plt.plot(n, t[:,1]/number, label='ADJ')
> plt.plot(n, t[:,2]/number, label='LU')
> plt.plot(n, t[:,3]/number, label='bareis_det')
> plt.plot(n, t[:,4]/number, label='berkowitz_det.expand()')
> plt.legend(loc=0)
> plt.title('Average time to complete 1 matrix inversion/determinant')
> plt.xlabel('matrix dimension')
> plt.ylabel('Time [seconds]')
> plt.xticks(n)
> plt.axis([2, n[-1], 0, max(t/number)])
> plt.show()
>
>
> I'd be curious to know if others get similar results as I do. I
> posted the results of the above script here:
> http://i35.tinypic.com/so09hw.jpg
>
> It looks like Gaussian elimination is suffering from the bottleneck in
> the Bareis determinant since it has an assertion that calls .det() to
> make sure the determinant
[sympy] Re: Symbolic matrix inversion
On Sep 30, 8:39 am, Luke wrote: > The methods you suggest essentially takes care of the mass matrix > problem by solving a linear system numerically during numerical > integration. I am familiar with tools out there that do this, but > this isn't what I'm looking to do. I haven't seen one that is written > directly usable in Python -- do you know of one? The netlib packages > have this capability, but I'm no Fortran programmer. > > What I am interested in doing is solving the linear system > symbolically so that first order equations can be generated > symbolically and the most generic of ODE solvers will work. This also > eliminates the iteration that is being done by the ODE solver during > time integration. > By eliminating iteration here I mean eliminating the iterations of the root finding / linear system solving algorithm, not adaptive time stepping type iteration. ~Luke > Thanks, > ~Luke > > On Sep 29, 8:07 pm, Tim Lahey wrote: > > > On Sep 29, 2009, at 7:15 PM, Alan Bromborsky wrote: > > > > Are there differential equation solvers where you don't have to invert > > > the matrix? > > > A Newmark-Beta scheme will directly solve a second-order system of ODEs. > > The standard form uses iteration to solve the system so no inversion is > > necessary. For linear second-order problems you can rewrite things to > > use matrix algebra. > > > For more information, I recommend Bathe and Wilson, > > > Klaus-Jürgen Bathe and Edward L. Wilson. Numerical Methods in Finite > > Element Analysis. Prentice Hall, Englewood Cliffs, New Jersey, > > 1976. > > > There are other second order solvers out there too. > > > Cheers, > > > Tim. > > > --- > > Tim Lahey > > PhD Candidate, Systems Design Engineering > > University of Waterloohttp://www.linkedin.com/in/timlahey > > --~--~-~--~~~---~--~~ You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to sympy@googlegroups.com To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sympy?hl=en -~--~~~~--~~--~--~---
[sympy] Re: Symbolic matrix inversion
Alan, For the systems I have studied at the moment, the most complicated inverses I have need to compute are 3x3 dense matrices (for the nonlinear equations of motion of a benchmark bicycle model [0]) and 6x6 sparse inverses (for solving the kinematic equations of motion for the derivatives of the configuration variables in terms of the generalized speeds). For an example, you can look at my github in the bicycle_work branch and look at examples/rollingdisc/ rollingdisc_lib.py. All the functions in that file were generated using Sympy/PyDy and include the equations of motion in first order form. In that particular example, the mass matrix is diagonal, so inversion is trivial. In general though this is not the case and for larger dimensional systems, this symbolic inversion could literally take days to complete and this is what I want to avoid. Again, the reason I want to do this symbolically is so that the generated equations of motion are in first order form and no special ODE solver is needed. As far as I know the Scipy numerical integration routines and the GSL routines require equations in this form, so I would like to keep my equations in a form that can be immediately useful to those audiences with minimal hassle. ~Luke [0] -- J. P. Meijaard, Jim M. Papadopoulos, Andy Ruina, A. L. Schwab, 2007 ``Linearized dynamics equations for the balance and steer of a bicycle: a benchmark and review,'' Proceedings of the Royal Society A 463:1955-1982. On Sep 29, 7:34 pm, Alan Bromborsky wrote: > Luke wrote: > > I also check the GSL (GNU Scientific Library). They have a nice > > numerical integrator, but it doesn't allow for a mass matrix. > > ~Luke > > On Sep 29, 4:44 pm, Luke wrote: > > >> Yes, this is something I should look into. I am pretty sure that the > >> Netlib codes have this functionality, but it hasn't been wrapped into > >> the Python scientific packages that I know of, at least not yet. > >> scipy.integrate has odeint and ode, but both need everything in first > >> order form, no mass matrix is allowed. > > >> I figured that if I was going to work symbolically, I may as well go > >> as far as I can go, so if I can invert the matrices symbolically, I'd > >> prefer that. This would also be nice because then the equations could > >> be integrated with any numerical integrator out there, as the ode > >> formulation would be as generic as possible. > > >> ~Luke > > >> On Sep 29, 4:15 pm, Alan Bromborsky wrote: > > >>> Luke wrote: > > On Sep 29, 1:09 pm, Ondrej Certik wrote: > > > On Tue, Sep 29, 2009 at 12:49 PM, Luke wrote: > > >> I'm using Sympy from within PyDy to generate the equations of motion > >> for > >> mechanical systems. At the end of the day, the equations can be most > >> generally written as: > >> M(x) * x'' = F(x, x', t) > >> M(x) is what is known as the mass matrix, and will in general depend > >> on the > >> configuration of the system (positions and angles). This matrix needs > >> to be > >> inverted in order to solve for x'', which then will allow for numerical > >> integration or stability analysis. > >> I am generating M(x) symbolically, and in some case is it sparse, but > >> in > >> many cases, it isn't. Each entry of the matrix is a Sympy expression, > >> but > >> instead of trying to invert a matrix of Sympy expressions, I introduce > >> a > >> dummy symbol for each non-zero entries and then invert the matrix of > >> dummy > >> symbols. Humble systems might have 5 degrees of freedom (so a 5x5 > >> needs to > >> be inverted), so this inversion isn't so bad, but beyond this, the > >> matrix > >> inversions take a really long time, especially when the matrices are > >> full > >> (no zero entries). > >> I was thinking that it might be nice to have pre-computed matrix > >> inverses > >> for n x n matrices. Matrix inversion is O(n^3), so it would be nice > >> to have > >> all this precomputed symbolically, and this would greatly speed up > >> Sympy's > >> matrix capabilities. Inverses up to say 100x100 could be computed (or > >> maybe > >> something smaller), and then when you need an inverse, everything > >> would be > >> fast. This could also be used behind the scenes (by introduction of > >> symbolic substitution dictionaries) for inverting a matrix full of > >> sympy > >> expressions. > > > The inversion of a 100x100 dense matrix would be quite a big mess, > > wouldn't it? > > Yes, it would be disaster. > > > But I see what you mean, and I think it makes sense to cache it, if it > > speeds things up. But see below. > > >> Does anybody know if this has been done by somebody somewhere, or have > >> any > >> other ideas on how it could be done better than the way I suggested? > > > I would first try to profile the inversion code to se
[sympy] Re: Symbolic matrix inversion
The methods you suggest essentially takes care of the mass matrix problem by solving a linear system numerically during numerical integration. I am familiar with tools out there that do this, but this isn't what I'm looking to do. I haven't seen one that is written directly usable in Python -- do you know of one? The netlib packages have this capability, but I'm no Fortran programmer. What I am interested in doing is solving the linear system symbolically so that first order equations can be generated symbolically and the most generic of ODE solvers will work. This also eliminates the iteration that is being done by the ODE solver during time integration. Thanks, ~Luke On Sep 29, 8:07 pm, Tim Lahey wrote: > On Sep 29, 2009, at 7:15 PM, Alan Bromborsky wrote: > > > Are there differential equation solvers where you don't have to invert > > the matrix? > > A Newmark-Beta scheme will directly solve a second-order system of ODEs. > The standard form uses iteration to solve the system so no inversion is > necessary. For linear second-order problems you can rewrite things to > use matrix algebra. > > For more information, I recommend Bathe and Wilson, > > Klaus-Jürgen Bathe and Edward L. Wilson. Numerical Methods in Finite > Element Analysis. Prentice Hall, Englewood Cliffs, New Jersey, > 1976. > > There are other second order solvers out there too. > > Cheers, > > Tim. > > --- > Tim Lahey > PhD Candidate, Systems Design Engineering > University of Waterloohttp://www.linkedin.com/in/timlahey --~--~-~--~~~---~--~~ You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to sympy@googlegroups.com To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sympy?hl=en -~--~~~~--~~--~--~---
[sympy] Re: Symbolic matrix inversion
On Sep 29, 2009, at 7:15 PM, Alan Bromborsky wrote: > Are there differential equation solvers where you don't have to invert > the matrix? A Newmark-Beta scheme will directly solve a second-order system of ODEs. The standard form uses iteration to solve the system so no inversion is necessary. For linear second-order problems you can rewrite things to use matrix algebra. For more information, I recommend Bathe and Wilson, Klaus-Jürgen Bathe and Edward L. Wilson. Numerical Methods in Finite Element Analysis. Prentice Hall, Englewood Cliffs, New Jersey, 1976. There are other second order solvers out there too. Cheers, Tim. --- Tim Lahey PhD Candidate, Systems Design Engineering University of Waterloo http://www.linkedin.com/in/timlahey --~--~-~--~~~---~--~~ You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to sympy@googlegroups.com To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sympy?hl=en -~--~~~~--~~--~--~---
[sympy] Re: Symbolic matrix inversion
Luke wrote: > I also check the GSL (GNU Scientific Library). They have a nice > numerical integrator, but it doesn't allow for a mass matrix. > ~Luke > On Sep 29, 4:44 pm, Luke wrote: > >> Yes, this is something I should look into. I am pretty sure that the >> Netlib codes have this functionality, but it hasn't been wrapped into >> the Python scientific packages that I know of, at least not yet. >> scipy.integrate has odeint and ode, but both need everything in first >> order form, no mass matrix is allowed. >> >> I figured that if I was going to work symbolically, I may as well go >> as far as I can go, so if I can invert the matrices symbolically, I'd >> prefer that. This would also be nice because then the equations could >> be integrated with any numerical integrator out there, as the ode >> formulation would be as generic as possible. >> >> ~Luke >> >> On Sep 29, 4:15 pm, Alan Bromborsky wrote: >> >> >> >> >>> Luke wrote: >>> On Sep 29, 1:09 pm, Ondrej Certik wrote: > On Tue, Sep 29, 2009 at 12:49 PM, Luke wrote: > >> I'm using Sympy from within PyDy to generate the equations of motion for >> mechanical systems. At the end of the day, the equations can be most >> generally written as: >> M(x) * x'' = F(x, x', t) >> M(x) is what is known as the mass matrix, and will in general depend on >> the >> configuration of the system (positions and angles). This matrix needs >> to be >> inverted in order to solve for x'', which then will allow for numerical >> integration or stability analysis. >> I am generating M(x) symbolically, and in some case is it sparse, but in >> many cases, it isn't. Each entry of the matrix is a Sympy expression, >> but >> instead of trying to invert a matrix of Sympy expressions, I introduce a >> dummy symbol for each non-zero entries and then invert the matrix of >> dummy >> symbols. Humble systems might have 5 degrees of freedom (so a 5x5 needs >> to >> be inverted), so this inversion isn't so bad, but beyond this, the matrix >> inversions take a really long time, especially when the matrices are full >> (no zero entries). >> I was thinking that it might be nice to have pre-computed matrix inverses >> for n x n matrices. Matrix inversion is O(n^3), so it would be nice to >> have >> all this precomputed symbolically, and this would greatly speed up >> Sympy's >> matrix capabilities. Inverses up to say 100x100 could be computed (or >> maybe >> something smaller), and then when you need an inverse, everything would >> be >> fast. This could also be used behind the scenes (by introduction of >> symbolic substitution dictionaries) for inverting a matrix full of sympy >> expressions. >> > The inversion of a 100x100 dense matrix would be quite a big mess, > wouldn't it? > Yes, it would be disaster. > But I see what you mean, and I think it makes sense to cache it, if it > speeds things up. But see below. > >> Does anybody know if this has been done by somebody somewhere, or have >> any >> other ideas on how it could be done better than the way I suggested? >> > I would first try to profile the inversion code to see why it is slow. > Because for example the adjoint method is just a ratio of two > determinants, so it may be the determinants calculation that is slow > (and should be cached). However, there are 100^2 different > determinants (right?), but I guess caching 1 expressions should > still be doable. But if this is the case, we should just speed up the > determinant calculation, imho. Looking at the code, we have 2 > algorithms implemented, bareis and berkowitz. > I did some timings of the three matrix inversion (GE, ADJ, LU) using the timeit module. I also timed the two determinant methods as well to see they the two perform side by side. It seems that bareis (the default one) is drastically slower than berkowitz, even when you call .expand() on the berkowitz determinant to make it end up with identical symbolic expressions. Maybe this should be the default method for .det()? I went up to 5x5 matrices, it started to get really slow after that. Gaussian elimination seems to be the slowest one, at least for dense matrices like the ones I used. Here is the code I used to do the timings: import timeit from numpy import zeros, max import matplotlib.pyplot as plt # Dimension of matrix to invert n = range(2, 6) # Number of times to invert number = 20 # Store the results t = zeros((len(n), 5)) for i in n: setup_code = "from sympy import Matrix, Symbol\nM = Matrix(%d,\
[sympy] Re: Symbolic matrix inversion
I also check the GSL (GNU Scientific Library). They have a nice
numerical integrator, but it doesn't allow for a mass matrix.
~Luke
On Sep 29, 4:44 pm, Luke wrote:
> Yes, this is something I should look into. I am pretty sure that the
> Netlib codes have this functionality, but it hasn't been wrapped into
> the Python scientific packages that I know of, at least not yet.
> scipy.integrate has odeint and ode, but both need everything in first
> order form, no mass matrix is allowed.
>
> I figured that if I was going to work symbolically, I may as well go
> as far as I can go, so if I can invert the matrices symbolically, I'd
> prefer that. This would also be nice because then the equations could
> be integrated with any numerical integrator out there, as the ode
> formulation would be as generic as possible.
>
> ~Luke
>
> On Sep 29, 4:15 pm, Alan Bromborsky wrote:
>
>
>
> > Luke wrote:
>
> > > On Sep 29, 1:09 pm, Ondrej Certik wrote:
>
> > >> On Tue, Sep 29, 2009 at 12:49 PM, Luke wrote:
>
> > >>> I'm using Sympy from within PyDy to generate the equations of motion for
> > >>> mechanical systems. At the end of the day, the equations can be most
> > >>> generally written as:
> > >>> M(x) * x'' = F(x, x', t)
> > >>> M(x) is what is known as the mass matrix, and will in general depend on
> > >>> the
> > >>> configuration of the system (positions and angles). This matrix needs
> > >>> to be
> > >>> inverted in order to solve for x'', which then will allow for numerical
> > >>> integration or stability analysis.
> > >>> I am generating M(x) symbolically, and in some case is it sparse, but in
> > >>> many cases, it isn't. Each entry of the matrix is a Sympy expression,
> > >>> but
> > >>> instead of trying to invert a matrix of Sympy expressions, I introduce a
> > >>> dummy symbol for each non-zero entries and then invert the matrix of
> > >>> dummy
> > >>> symbols. Humble systems might have 5 degrees of freedom (so a 5x5
> > >>> needs to
> > >>> be inverted), so this inversion isn't so bad, but beyond this, the
> > >>> matrix
> > >>> inversions take a really long time, especially when the matrices are
> > >>> full
> > >>> (no zero entries).
> > >>> I was thinking that it might be nice to have pre-computed matrix
> > >>> inverses
> > >>> for n x n matrices. Matrix inversion is O(n^3), so it would be nice to
> > >>> have
> > >>> all this precomputed symbolically, and this would greatly speed up
> > >>> Sympy's
> > >>> matrix capabilities. Inverses up to say 100x100 could be computed (or
> > >>> maybe
> > >>> something smaller), and then when you need an inverse, everything would
> > >>> be
> > >>> fast. This could also be used behind the scenes (by introduction of
> > >>> symbolic substitution dictionaries) for inverting a matrix full of sympy
> > >>> expressions.
>
> > >> The inversion of a 100x100 dense matrix would be quite a big mess,
> > >> wouldn't it?
>
> > > Yes, it would be disaster.
>
> > >> But I see what you mean, and I think it makes sense to cache it, if it
> > >> speeds things up. But see below.
>
> > >>> Does anybody know if this has been done by somebody somewhere, or have
> > >>> any
> > >>> other ideas on how it could be done better than the way I suggested?
>
> > >> I would first try to profile the inversion code to see why it is slow.
> > >> Because for example the adjoint method is just a ratio of two
> > >> determinants, so it may be the determinants calculation that is slow
> > >> (and should be cached). However, there are 100^2 different
> > >> determinants (right?), but I guess caching 1 expressions should
> > >> still be doable. But if this is the case, we should just speed up the
> > >> determinant calculation, imho. Looking at the code, we have 2
> > >> algorithms implemented, bareis and berkowitz.
>
> > > I did some timings of the three matrix inversion (GE, ADJ, LU) using
> > > the timeit module. I also timed the two determinant methods as well
> > > to see they the two perform side by side. It seems that bareis (the
> > > default one) is drastically slower than berkowitz, even when you
> > > call .expand() on the berkowitz determinant to make it end up with
> > > identical symbolic expressions. Maybe this should be the default
> > > method for .det()? I went up to 5x5 matrices, it started to get
> > > really slow after that. Gaussian elimination seems to be the slowest
> > > one, at least for dense matrices like the ones I used. Here is the
> > > code I used to do the timings:
>
> > > import timeit
> > > from numpy import zeros, max
> > > import matplotlib.pyplot as plt
>
> > > # Dimension of matrix to invert
> > > n = range(2, 6)
> > > # Number of times to invert
> > > number = 20
> > > # Store the results
> > > t = zeros((len(n), 5))
>
> > > for i in n:
> > > setup_code = "from sympy import Matrix, Symbol\nM = Matrix(%d,\
> > > %d"%(i,i)+",lambda i,j: Symbol('m%d%d'%(i,j)))"
> > > t[i-2, 0] = timeit.Timer('M.inverse_GE()', se
[sympy] Re: Symbolic matrix inversion
In the formulation I use for PyDy, the equations of motion are generated in first order form. For holonomic systems with n degrees of freedom, there will be 2n first order equations and the first n of these I refer to as the kinematic differential equations. In the simplest case, the form of the first n equations would be like: q_i' = u_i (i = 1,...,n) The remaining n ODE's (dynamic differential equations) would be the ones with the mass matrix in the form like: M(q)*u' = F(q, u, t) Where q, u are n x 1 vectors, M is nxn, and F: n x n x 1 --> n Usually this choice of 'generalized speeds' (a.k.a quasi-coordinates) which result in kinematic differential equations like: "q_i' = u_i" is not the best choice out there. But even with other choices, they can be written as: q' = T(q) * u + R(t) Where in most cases R(t) is zero, except when you have some prescribed motion in your system. The dynamic differential equations would still have the same form as before. So they are linear in the generalized speeds, but the T(q) matrix would in general be comprised of trigonometric entries which depend on the configuration and/or parameters. Such choices of generalized speeds can often lead to reduced complexity of the dynamic equations, with some (usually worthwhile) trade-off in the complexity of the kinematic equations. For nonholonomic systems, depending on how you formulate things you may end up having to integrate fewer dynamic equations. For example, in the rolling disc, there are 5 coordinates that configure the disc, but only 3 degrees of freedom. If you are interested in all of the motion, you would integrate 5 kinematic D.E.'s and 3 dynamic D.E.'s. If you only are interested in the dynamic equations, you would end up only integrating 1 kinematic D.E. (the one associated with the lean of the disc), and 3 dynamic D.E.'s. So in these types of systems, it depends on what your needs are. ~Luke On Sep 29, 4:22 pm, Alan Bromborsky wrote: > Alan Bromborsky wrote: > > Luke wrote: > > >> On Sep 29, 1:09 pm, Ondrej Certik wrote: > > >>> On Tue, Sep 29, 2009 at 12:49 PM, Luke wrote: > > I'm using Sympy from within PyDy to generate the equations of motion for > mechanical systems. At the end of the day, the equations can be most > generally written as: > M(x) * x'' = F(x, x', t) > M(x) is what is known as the mass matrix, and will in general depend on > the > configuration of the system (positions and angles). This matrix needs > to be > inverted in order to solve for x'', which then will allow for numerical > integration or stability analysis. > I am generating M(x) symbolically, and in some case is it sparse, but in > many cases, it isn't. Each entry of the matrix is a Sympy expression, > but > instead of trying to invert a matrix of Sympy expressions, I introduce a > dummy symbol for each non-zero entries and then invert the matrix of > dummy > symbols. Humble systems might have 5 degrees of freedom (so a 5x5 needs > to > be inverted), so this inversion isn't so bad, but beyond this, the matrix > inversions take a really long time, especially when the matrices are full > (no zero entries). > I was thinking that it might be nice to have pre-computed matrix inverses > for n x n matrices. Matrix inversion is O(n^3), so it would be nice to > have > all this precomputed symbolically, and this would greatly speed up > Sympy's > matrix capabilities. Inverses up to say 100x100 could be computed (or > maybe > something smaller), and then when you need an inverse, everything would > be > fast. This could also be used behind the scenes (by introduction of > symbolic substitution dictionaries) for inverting a matrix full of sympy > expressions. > > >>> The inversion of a 100x100 dense matrix would be quite a big mess, > >>> wouldn't it? > > >> Yes, it would be disaster. > > >>> But I see what you mean, and I think it makes sense to cache it, if it > >>> speeds things up. But see below. > > Does anybody know if this has been done by somebody somewhere, or have > any > other ideas on how it could be done better than the way I suggested? > > >>> I would first try to profile the inversion code to see why it is slow. > >>> Because for example the adjoint method is just a ratio of two > >>> determinants, so it may be the determinants calculation that is slow > >>> (and should be cached). However, there are 100^2 different > >>> determinants (right?), but I guess caching 1 expressions should > >>> still be doable. But if this is the case, we should just speed up the > >>> determinant calculation, imho. Looking at the code, we have 2 > >>> algorithms implemented, bareis and berkowitz. > > >> I did some timings of the three matrix inversion (GE, ADJ, LU) using > >> the timeit module. I also timed the two determinant met
[sympy] Re: Symbolic matrix inversion
Yes, this is something I should look into. I am pretty sure that the
Netlib codes have this functionality, but it hasn't been wrapped into
the Python scientific packages that I know of, at least not yet.
scipy.integrate has odeint and ode, but both need everything in first
order form, no mass matrix is allowed.
I figured that if I was going to work symbolically, I may as well go
as far as I can go, so if I can invert the matrices symbolically, I'd
prefer that. This would also be nice because then the equations could
be integrated with any numerical integrator out there, as the ode
formulation would be as generic as possible.
~Luke
On Sep 29, 4:15 pm, Alan Bromborsky wrote:
> Luke wrote:
>
> > On Sep 29, 1:09 pm, Ondrej Certik wrote:
>
> >> On Tue, Sep 29, 2009 at 12:49 PM, Luke wrote:
>
> >>> I'm using Sympy from within PyDy to generate the equations of motion for
> >>> mechanical systems. At the end of the day, the equations can be most
> >>> generally written as:
> >>> M(x) * x'' = F(x, x', t)
> >>> M(x) is what is known as the mass matrix, and will in general depend on
> >>> the
> >>> configuration of the system (positions and angles). This matrix needs to
> >>> be
> >>> inverted in order to solve for x'', which then will allow for numerical
> >>> integration or stability analysis.
> >>> I am generating M(x) symbolically, and in some case is it sparse, but in
> >>> many cases, it isn't. Each entry of the matrix is a Sympy expression, but
> >>> instead of trying to invert a matrix of Sympy expressions, I introduce a
> >>> dummy symbol for each non-zero entries and then invert the matrix of dummy
> >>> symbols. Humble systems might have 5 degrees of freedom (so a 5x5 needs
> >>> to
> >>> be inverted), so this inversion isn't so bad, but beyond this, the matrix
> >>> inversions take a really long time, especially when the matrices are full
> >>> (no zero entries).
> >>> I was thinking that it might be nice to have pre-computed matrix inverses
> >>> for n x n matrices. Matrix inversion is O(n^3), so it would be nice to
> >>> have
> >>> all this precomputed symbolically, and this would greatly speed up Sympy's
> >>> matrix capabilities. Inverses up to say 100x100 could be computed (or
> >>> maybe
> >>> something smaller), and then when you need an inverse, everything would be
> >>> fast. This could also be used behind the scenes (by introduction of
> >>> symbolic substitution dictionaries) for inverting a matrix full of sympy
> >>> expressions.
>
> >> The inversion of a 100x100 dense matrix would be quite a big mess,
> >> wouldn't it?
>
> > Yes, it would be disaster.
>
> >> But I see what you mean, and I think it makes sense to cache it, if it
> >> speeds things up. But see below.
>
> >>> Does anybody know if this has been done by somebody somewhere, or have any
> >>> other ideas on how it could be done better than the way I suggested?
>
> >> I would first try to profile the inversion code to see why it is slow.
> >> Because for example the adjoint method is just a ratio of two
> >> determinants, so it may be the determinants calculation that is slow
> >> (and should be cached). However, there are 100^2 different
> >> determinants (right?), but I guess caching 1 expressions should
> >> still be doable. But if this is the case, we should just speed up the
> >> determinant calculation, imho. Looking at the code, we have 2
> >> algorithms implemented, bareis and berkowitz.
>
> > I did some timings of the three matrix inversion (GE, ADJ, LU) using
> > the timeit module. I also timed the two determinant methods as well
> > to see they the two perform side by side. It seems that bareis (the
> > default one) is drastically slower than berkowitz, even when you
> > call .expand() on the berkowitz determinant to make it end up with
> > identical symbolic expressions. Maybe this should be the default
> > method for .det()? I went up to 5x5 matrices, it started to get
> > really slow after that. Gaussian elimination seems to be the slowest
> > one, at least for dense matrices like the ones I used. Here is the
> > code I used to do the timings:
>
> > import timeit
> > from numpy import zeros, max
> > import matplotlib.pyplot as plt
>
> > # Dimension of matrix to invert
> > n = range(2, 6)
> > # Number of times to invert
> > number = 20
> > # Store the results
> > t = zeros((len(n), 5))
>
> > for i in n:
> > setup_code = "from sympy import Matrix, Symbol\nM = Matrix(%d,\
> > %d"%(i,i)+",lambda i,j: Symbol('m%d%d'%(i,j)))"
> > t[i-2, 0] = timeit.Timer('M.inverse_GE()', setup_code).timeit
> > (number)
> > t[i-2, 1] = timeit.Timer('M.inverse_ADJ()', setup_code).timeit
> > (number)
> > t[i-2, 2] = timeit.Timer('M.inverse_LU()', setup_code).timeit
> > (number)
> > t[i-2, 3] = timeit.Timer('M.det(method="bareis")',
> > setup_code).timeit(number)
> > t[i-2, 4] = timeit.Timer('M.det(method="berkowitz").expand()',
> > setup_code).timeit(number
[sympy] Re: Symbolic matrix inversion
Alan Bromborsky wrote:
> Luke wrote:
>
>> On Sep 29, 1:09 pm, Ondrej Certik wrote:
>>
>>
>>> On Tue, Sep 29, 2009 at 12:49 PM, Luke wrote:
>>>
>>>
I'm using Sympy from within PyDy to generate the equations of motion for
mechanical systems. At the end of the day, the equations can be most
generally written as:
M(x) * x'' = F(x, x', t)
M(x) is what is known as the mass matrix, and will in general depend on the
configuration of the system (positions and angles). This matrix needs to
be
inverted in order to solve for x'', which then will allow for numerical
integration or stability analysis.
I am generating M(x) symbolically, and in some case is it sparse, but in
many cases, it isn't. Each entry of the matrix is a Sympy expression, but
instead of trying to invert a matrix of Sympy expressions, I introduce a
dummy symbol for each non-zero entries and then invert the matrix of dummy
symbols. Humble systems might have 5 degrees of freedom (so a 5x5 needs to
be inverted), so this inversion isn't so bad, but beyond this, the matrix
inversions take a really long time, especially when the matrices are full
(no zero entries).
I was thinking that it might be nice to have pre-computed matrix inverses
for n x n matrices. Matrix inversion is O(n^3), so it would be nice to
have
all this precomputed symbolically, and this would greatly speed up Sympy's
matrix capabilities. Inverses up to say 100x100 could be computed (or
maybe
something smaller), and then when you need an inverse, everything would be
fast. This could also be used behind the scenes (by introduction of
symbolic substitution dictionaries) for inverting a matrix full of sympy
expressions.
>>> The inversion of a 100x100 dense matrix would be quite a big mess, wouldn't
>>> it?
>>>
>>>
>>>
>> Yes, it would be disaster.
>>
>>
>>
>>> But I see what you mean, and I think it makes sense to cache it, if it
>>> speeds things up. But see below.
>>>
>>>
>>>
Does anybody know if this has been done by somebody somewhere, or have any
other ideas on how it could be done better than the way I suggested?
>>> I would first try to profile the inversion code to see why it is slow.
>>> Because for example the adjoint method is just a ratio of two
>>> determinants, so it may be the determinants calculation that is slow
>>> (and should be cached). However, there are 100^2 different
>>> determinants (right?), but I guess caching 1 expressions should
>>> still be doable. But if this is the case, we should just speed up the
>>> determinant calculation, imho. Looking at the code, we have 2
>>> algorithms implemented, bareis and berkowitz.
>>>
>>>
>>>
>> I did some timings of the three matrix inversion (GE, ADJ, LU) using
>> the timeit module. I also timed the two determinant methods as well
>> to see they the two perform side by side. It seems that bareis (the
>> default one) is drastically slower than berkowitz, even when you
>> call .expand() on the berkowitz determinant to make it end up with
>> identical symbolic expressions. Maybe this should be the default
>> method for .det()? I went up to 5x5 matrices, it started to get
>> really slow after that. Gaussian elimination seems to be the slowest
>> one, at least for dense matrices like the ones I used. Here is the
>> code I used to do the timings:
>>
>> import timeit
>> from numpy import zeros, max
>> import matplotlib.pyplot as plt
>>
>> # Dimension of matrix to invert
>> n = range(2, 6)
>> # Number of times to invert
>> number = 20
>> # Store the results
>> t = zeros((len(n), 5))
>>
>> for i in n:
>> setup_code = "from sympy import Matrix, Symbol\nM = Matrix(%d,\
>> %d"%(i,i)+",lambda i,j: Symbol('m%d%d'%(i,j)))"
>> t[i-2, 0] = timeit.Timer('M.inverse_GE()', setup_code).timeit
>> (number)
>> t[i-2, 1] = timeit.Timer('M.inverse_ADJ()', setup_code).timeit
>> (number)
>> t[i-2, 2] = timeit.Timer('M.inverse_LU()', setup_code).timeit
>> (number)
>> t[i-2, 3] = timeit.Timer('M.det(method="bareis")',
>> setup_code).timeit(number)
>> t[i-2, 4] = timeit.Timer('M.det(method="berkowitz").expand()',
>> setup_code).timeit(number)
>>
>>
>> plt.plot(n, t[:,0]/number, label='GE')
>> plt.plot(n, t[:,1]/number, label='ADJ')
>> plt.plot(n, t[:,2]/number, label='LU')
>> plt.plot(n, t[:,3]/number, label='bareis_det')
>> plt.plot(n, t[:,4]/number, label='berkowitz_det.expand()')
>> plt.legend(loc=0)
>> plt.title('Average time to complete 1 matrix inversion/determinant')
>> plt.xlabel('matrix dimension')
>> plt.ylabel('Time [seconds]')
>> plt.xticks(n)
>> plt.axis([2, n[-1], 0, max(t/number)])
>> plt.show()
>>
>>
>> I'd be curious to know if others get similar results as I do. I
>> posted the results of the above sc
[sympy] Re: Symbolic matrix inversion
Luke wrote:
>
> On Sep 29, 1:09 pm, Ondrej Certik wrote:
>
>> On Tue, Sep 29, 2009 at 12:49 PM, Luke wrote:
>>
>>> I'm using Sympy from within PyDy to generate the equations of motion for
>>> mechanical systems. At the end of the day, the equations can be most
>>> generally written as:
>>> M(x) * x'' = F(x, x', t)
>>> M(x) is what is known as the mass matrix, and will in general depend on the
>>> configuration of the system (positions and angles). This matrix needs to be
>>> inverted in order to solve for x'', which then will allow for numerical
>>> integration or stability analysis.
>>> I am generating M(x) symbolically, and in some case is it sparse, but in
>>> many cases, it isn't. Each entry of the matrix is a Sympy expression, but
>>> instead of trying to invert a matrix of Sympy expressions, I introduce a
>>> dummy symbol for each non-zero entries and then invert the matrix of dummy
>>> symbols. Humble systems might have 5 degrees of freedom (so a 5x5 needs to
>>> be inverted), so this inversion isn't so bad, but beyond this, the matrix
>>> inversions take a really long time, especially when the matrices are full
>>> (no zero entries).
>>> I was thinking that it might be nice to have pre-computed matrix inverses
>>> for n x n matrices. Matrix inversion is O(n^3), so it would be nice to have
>>> all this precomputed symbolically, and this would greatly speed up Sympy's
>>> matrix capabilities. Inverses up to say 100x100 could be computed (or maybe
>>> something smaller), and then when you need an inverse, everything would be
>>> fast. This could also be used behind the scenes (by introduction of
>>> symbolic substitution dictionaries) for inverting a matrix full of sympy
>>> expressions.
>>>
>> The inversion of a 100x100 dense matrix would be quite a big mess, wouldn't
>> it?
>>
>>
>
> Yes, it would be disaster.
>
>
>> But I see what you mean, and I think it makes sense to cache it, if it
>> speeds things up. But see below.
>>
>>
>>> Does anybody know if this has been done by somebody somewhere, or have any
>>> other ideas on how it could be done better than the way I suggested?
>>>
>> I would first try to profile the inversion code to see why it is slow.
>> Because for example the adjoint method is just a ratio of two
>> determinants, so it may be the determinants calculation that is slow
>> (and should be cached). However, there are 100^2 different
>> determinants (right?), but I guess caching 1 expressions should
>> still be doable. But if this is the case, we should just speed up the
>> determinant calculation, imho. Looking at the code, we have 2
>> algorithms implemented, bareis and berkowitz.
>>
>>
>
> I did some timings of the three matrix inversion (GE, ADJ, LU) using
> the timeit module. I also timed the two determinant methods as well
> to see they the two perform side by side. It seems that bareis (the
> default one) is drastically slower than berkowitz, even when you
> call .expand() on the berkowitz determinant to make it end up with
> identical symbolic expressions. Maybe this should be the default
> method for .det()? I went up to 5x5 matrices, it started to get
> really slow after that. Gaussian elimination seems to be the slowest
> one, at least for dense matrices like the ones I used. Here is the
> code I used to do the timings:
>
> import timeit
> from numpy import zeros, max
> import matplotlib.pyplot as plt
>
> # Dimension of matrix to invert
> n = range(2, 6)
> # Number of times to invert
> number = 20
> # Store the results
> t = zeros((len(n), 5))
>
> for i in n:
> setup_code = "from sympy import Matrix, Symbol\nM = Matrix(%d,\
> %d"%(i,i)+",lambda i,j: Symbol('m%d%d'%(i,j)))"
> t[i-2, 0] = timeit.Timer('M.inverse_GE()', setup_code).timeit
> (number)
> t[i-2, 1] = timeit.Timer('M.inverse_ADJ()', setup_code).timeit
> (number)
> t[i-2, 2] = timeit.Timer('M.inverse_LU()', setup_code).timeit
> (number)
> t[i-2, 3] = timeit.Timer('M.det(method="bareis")',
> setup_code).timeit(number)
> t[i-2, 4] = timeit.Timer('M.det(method="berkowitz").expand()',
> setup_code).timeit(number)
>
>
> plt.plot(n, t[:,0]/number, label='GE')
> plt.plot(n, t[:,1]/number, label='ADJ')
> plt.plot(n, t[:,2]/number, label='LU')
> plt.plot(n, t[:,3]/number, label='bareis_det')
> plt.plot(n, t[:,4]/number, label='berkowitz_det.expand()')
> plt.legend(loc=0)
> plt.title('Average time to complete 1 matrix inversion/determinant')
> plt.xlabel('matrix dimension')
> plt.ylabel('Time [seconds]')
> plt.xticks(n)
> plt.axis([2, n[-1], 0, max(t/number)])
> plt.show()
>
>
> I'd be curious to know if others get similar results as I do. I
> posted the results of the above script here:
> http://i35.tinypic.com/so09hw.jpg
>
> It looks like Gaussian elimination is suffering from the bottleneck in
> the Bareis determinant since it has an assertion that calls .det() to
> make sure the determinant is non-zero.
>
[sympy] Re: Symbolic matrix inversion
On Sep 29, 1:09 pm, Ondrej Certik wrote:
> On Tue, Sep 29, 2009 at 12:49 PM, Luke wrote:
> > I'm using Sympy from within PyDy to generate the equations of motion for
> > mechanical systems. At the end of the day, the equations can be most
> > generally written as:
> > M(x) * x'' = F(x, x', t)
> > M(x) is what is known as the mass matrix, and will in general depend on the
> > configuration of the system (positions and angles). This matrix needs to be
> > inverted in order to solve for x'', which then will allow for numerical
> > integration or stability analysis.
> > I am generating M(x) symbolically, and in some case is it sparse, but in
> > many cases, it isn't. Each entry of the matrix is a Sympy expression, but
> > instead of trying to invert a matrix of Sympy expressions, I introduce a
> > dummy symbol for each non-zero entries and then invert the matrix of dummy
> > symbols. Humble systems might have 5 degrees of freedom (so a 5x5 needs to
> > be inverted), so this inversion isn't so bad, but beyond this, the matrix
> > inversions take a really long time, especially when the matrices are full
> > (no zero entries).
> > I was thinking that it might be nice to have pre-computed matrix inverses
> > for n x n matrices. Matrix inversion is O(n^3), so it would be nice to have
> > all this precomputed symbolically, and this would greatly speed up Sympy's
> > matrix capabilities. Inverses up to say 100x100 could be computed (or maybe
> > something smaller), and then when you need an inverse, everything would be
> > fast. This could also be used behind the scenes (by introduction of
> > symbolic substitution dictionaries) for inverting a matrix full of sympy
> > expressions.
>
> The inversion of a 100x100 dense matrix would be quite a big mess, wouldn't
> it?
>
Yes, it would be disaster.
> But I see what you mean, and I think it makes sense to cache it, if it
> speeds things up. But see below.
>
> > Does anybody know if this has been done by somebody somewhere, or have any
> > other ideas on how it could be done better than the way I suggested?
>
> I would first try to profile the inversion code to see why it is slow.
> Because for example the adjoint method is just a ratio of two
> determinants, so it may be the determinants calculation that is slow
> (and should be cached). However, there are 100^2 different
> determinants (right?), but I guess caching 1 expressions should
> still be doable. But if this is the case, we should just speed up the
> determinant calculation, imho. Looking at the code, we have 2
> algorithms implemented, bareis and berkowitz.
>
I did some timings of the three matrix inversion (GE, ADJ, LU) using
the timeit module. I also timed the two determinant methods as well
to see they the two perform side by side. It seems that bareis (the
default one) is drastically slower than berkowitz, even when you
call .expand() on the berkowitz determinant to make it end up with
identical symbolic expressions. Maybe this should be the default
method for .det()? I went up to 5x5 matrices, it started to get
really slow after that. Gaussian elimination seems to be the slowest
one, at least for dense matrices like the ones I used. Here is the
code I used to do the timings:
import timeit
from numpy import zeros, max
import matplotlib.pyplot as plt
# Dimension of matrix to invert
n = range(2, 6)
# Number of times to invert
number = 20
# Store the results
t = zeros((len(n), 5))
for i in n:
setup_code = "from sympy import Matrix, Symbol\nM = Matrix(%d,\
%d"%(i,i)+",lambda i,j: Symbol('m%d%d'%(i,j)))"
t[i-2, 0] = timeit.Timer('M.inverse_GE()', setup_code).timeit
(number)
t[i-2, 1] = timeit.Timer('M.inverse_ADJ()', setup_code).timeit
(number)
t[i-2, 2] = timeit.Timer('M.inverse_LU()', setup_code).timeit
(number)
t[i-2, 3] = timeit.Timer('M.det(method="bareis")',
setup_code).timeit(number)
t[i-2, 4] = timeit.Timer('M.det(method="berkowitz").expand()',
setup_code).timeit(number)
plt.plot(n, t[:,0]/number, label='GE')
plt.plot(n, t[:,1]/number, label='ADJ')
plt.plot(n, t[:,2]/number, label='LU')
plt.plot(n, t[:,3]/number, label='bareis_det')
plt.plot(n, t[:,4]/number, label='berkowitz_det.expand()')
plt.legend(loc=0)
plt.title('Average time to complete 1 matrix inversion/determinant')
plt.xlabel('matrix dimension')
plt.ylabel('Time [seconds]')
plt.xticks(n)
plt.axis([2, n[-1], 0, max(t/number)])
plt.show()
I'd be curious to know if others get similar results as I do. I
posted the results of the above script here:
http://i35.tinypic.com/so09hw.jpg
It looks like Gaussian elimination is suffering from the bottleneck in
the Bareis determinant since it has an assertion that calls .det() to
make sure the determinant is non-zero.
> Also how about the creation of so many expressions, it could be slow too.
>
Yeah, this is true. It seems to me though that if you had all the
inverses computed for fully populated matrices, then if so
[sympy] Re: Symbolic matrix inversion
On Tue, Sep 29, 2009 at 12:49 PM, Luke wrote: > I'm using Sympy from within PyDy to generate the equations of motion for > mechanical systems. At the end of the day, the equations can be most > generally written as: > M(x) * x'' = F(x, x', t) > M(x) is what is known as the mass matrix, and will in general depend on the > configuration of the system (positions and angles). This matrix needs to be > inverted in order to solve for x'', which then will allow for numerical > integration or stability analysis. > I am generating M(x) symbolically, and in some case is it sparse, but in > many cases, it isn't. Each entry of the matrix is a Sympy expression, but > instead of trying to invert a matrix of Sympy expressions, I introduce a > dummy symbol for each non-zero entries and then invert the matrix of dummy > symbols. Humble systems might have 5 degrees of freedom (so a 5x5 needs to > be inverted), so this inversion isn't so bad, but beyond this, the matrix > inversions take a really long time, especially when the matrices are full > (no zero entries). > I was thinking that it might be nice to have pre-computed matrix inverses > for n x n matrices. Matrix inversion is O(n^3), so it would be nice to have > all this precomputed symbolically, and this would greatly speed up Sympy's > matrix capabilities. Inverses up to say 100x100 could be computed (or maybe > something smaller), and then when you need an inverse, everything would be > fast. This could also be used behind the scenes (by introduction of > symbolic substitution dictionaries) for inverting a matrix full of sympy > expressions. The inversion of a 100x100 dense matrix would be quite a big mess, wouldn't it? But I see what you mean, and I think it makes sense to cache it, if it speeds things up. But see below. > Does anybody know if this has been done by somebody somewhere, or have any > other ideas on how it could be done better than the way I suggested? I would first try to profile the inversion code to see why it is slow. Because for example the adjoint method is just a ratio of two determinants, so it may be the determinants calculation that is slow (and should be cached). However, there are 100^2 different determinants (right?), but I guess caching 1 expressions should still be doable. But if this is the case, we should just speed up the determinant calculation, imho. Looking at the code, we have 2 algorithms implemented, bareis and berkowitz. Also how about the creation of so many expressions, it could be slow too. So what I want to say is that we should profile it and determine exactly what things are slowing things down, then we should speed them up (or cache them). Ondrej --~--~-~--~~~---~--~~ You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to sympy@googlegroups.com To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sympy?hl=en -~--~~~~--~~--~--~---
