Symbolic matrix inversion

[Luke]

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.  

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?

Thanks,

~Luke

[Ondrej Certik]

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 10000 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

[Luke]

On Sep 29, 1:09 pm, Ondrej Certik <ond...@certik.cz> wrote:

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 10000 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 some of the
entries were zero these things would simplify some (or a lot), and
then those expressions would be used rather than the full ones.

Maybe it could be an optional module that gets imported in when you
need to work with large matrices and want to have inverses precomputed
so that you don't have to sit and wait for them for hours? That way
it wouldn't affect the normal behavior of Sympy, but could be imported
if desired.

> 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).

Should I use cProfile for this? Or do you recommend another profiler?

~Luke

>
> Ondrej

[Alan Bromborsky]

Are there differential equation solvers where you don't have to invert
the matrix?

[Alan Bromborsky]

Also, what does the system look like if written as a first order system
of 2n equations?

[Luke]

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

[Luke]

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

[Luke]

I also check the GSL (GNU Scientific Library). They have a nice
numerical integrator, but it doesn't allow for a mass matrix.
~Luke

[Alan Bromborsky]

Would you give an example set of equations you wish to solve?

[Tim Lahey]

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

[Luke]

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

[Luke]

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.

[Luke]

On Sep 30, 8:39 am, Luke <hazelnu...@gmail.com> 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

[Ondrej Certik]

I use either cprofile, or line_profile:

http://packages.python.org/line_profiler/

the line_profiler is nice, that it tells you which lines take so long.

So you need 3x3 matrices, those are fine as it is in sympy, isn't it?
And then you need 6x6 sparse matrices? Are they block diagonal?

Ondrej

[Vinzent Steinberg]

On Sep 29, 9:49 pm, Luke <hazelnu...@gmail.com> 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