# [ODE] LCP solution methods

Sergiy Migdalskiy migdalskiy at hotmail.com
Wed Jan 5 00:17:41 MST 2005

```Hello:

I'm working on an LCP solver, using JM^-1J^T method in Baraff's terminology
(the LCP w=Mx+q, w,x>=0, wx=0 with matrix M being non singular and strictly
positive definite). I have some pretty good progress, but I'm still far from
a high-quality physics package. I wanted to ask for any suggestions how to
implement iterative O(n) expected running time solver.

I'm looking at running up to 10000 contact problem in real time on the
modern high-end PC with low memory footprint (don't ask me why :)). I'm
pretty successful with this goal so far for frictionless balls :) but I hope
to make it work for everything. The collision detection is not my focus now,
as with this number of constraints the biggest problem is collision response
(that is, when I'm using something more complicated than balls)

The LCP matrix is pretty sparse, it's relatively easy to represent it. But
still, I have about 64000 float entries in the matrix itself, plus control
structures. In essence, it means I'm limited to 50-100 sweeps through the
matrix (matrix-vector multiplications). I tried constructing no-fill
Cholesky factors for the matrix, doesn't seem to work well. Obviously I
can't use any significant fill-in due to the prohibitive size of the
problem.

For some simpler configurations (a pile of 1000 balls, tightly packed), the
matrix is "almost diagonally dominant" in some sense, which means simple
gauss-seidel iteration works just fine and I only need 15-30 matrix-vector
multiplications, and it works beautifully fast.

For more complicated situations however, there are big off-diagonal
elements, which brings gauss-seidel relaxation scheme I'm using to solve the
LCP to its knees, just as the theory suggests. I tried block Gauss-Seidel,
but it has essentially the same problems with the non-diagonally-dominant
matrix.

Can anyone suggest a way to precondition the JM'J'-type matrix so that it
becomes diagonal dominant? I really want to stick to a very simple iterative
method, since it works so fast and well for simple configurations. Does
anyone know any projected iteration methods, besides projected
SOR/Gauss-Seidel, which work for positive-definite LCPs? I'd kill for a
projective analog of Conjugate Gradient for LCP! CG solves the analogous
linear problem (just find x:Mx+q=0) in 30 steps even without
preconditioning - and my budget is 50-100 steps.

Has anyone tried to use multigrid/multilevel to solve contact forces? The
literature on multigrid for LCP is sparse to say the least, any references
would be very much appreciated.

Another way is to use some simplex method. I'm not a specialist in this, but
to my best knowledge simplex methods (Lemke, criss-cross) don't work on
sparse matrices, because the matrix gets more and more filled in during
pivots. And memory is an issue for me. Am I missing something here? Has
anyone tried simplex methods to solve contact forces in big systems (i.e.
where it's impossible to fit the whole dense LCP in memory)?

Dantzig method requires keeping and updating a factor of variable size
(depending on the active basis). I think I heard ODE uses or used this?
Anyway, there's no way for me to keep a dense factor of a submatrix of a
10000x10000 JM'J' matrix with a budget of 50-100 matrix-vector
multiplications. Again, am I missing something? I didn't ever programmed
updating of a factor by adding/deleting a row/column, but I did program a
few sparse (incomplete) factors for this problem and even that was a huge
pain, since the quality of the preconditioner falls dramatically as I try to
shrink memory usage by raising drop tolerances or implementing some moderate
fill-in heuristics. In any case, I suspect Dantzig method requires too many
steps for a big (10000 vars) LCP. Am I wrong?

Another way to solve LCP is to use interior point methods. I experimented
with them using alternative formulations (QP with such objective as energy
minimization). It works pretty well, but again, I need at least 10
outer-loop steps in bad cases, each step consisting of forming a
preconditioner and solving at least one or two linear problems (e.g. for
predictor-corrector type algorithm). I can only use cheap preconditioners,
and it means at least 20-30 steps (very optimistically!) for each linear
problem (I used MinRes, LSQR, QMR, SQMR, CG, CGS to solve the linear
problem, all of them with a few different preconditioners - all just to make
sure they work as expected by theory; I used both the Augmented System and
Normal Equations approaches; AS worked better for me). So, in total I'm
looking at least at 400 matrix-vector computations, which is out of budget.
Besides, in reality, it'll be closer to 1200 due to ill-conditioned systems,
especially closer to smaller complementarity gaps when the interior point
converges to the solution. Plus the time to construct each preconditioner..

Oh yeah, to top it all, the nature of the engine is quite dynamic, and usage
of warm-start (using the solution from the previous frame) is hardly
possible to the extent that can make a big difference...

I'm starting to run out of sane ideas... Any suggestions would be greatly
appreciated.

Sincerely,
Sergiy
```