A new PGS solver

[ngbinh]

[This is a repost from ODE list]

Hi.

I stumble across this paper: "An Algorithm for the Approximate and Fast
Solution of Linear Complementarity Problems"

http://www-unix.mcs.anl.gov/~leyffer/li ... orales.pdf

It has a description of a new PGS that should work directly with ODE and (should be) more stable.
If any of you have time to implement it, I may help a bit in theory part.

[Antonio Martini]

[This is a repost from ODE list]

Hi.

I stumble across this paper: "An Algorithm for the Approximate and Fast
Solution of Linear Complementarity Problems"

http://www-unix.mcs.anl.gov/~leyffer/li ... orales.pdf

It has a description of a new PGS that should work directly with ODE and (should be) more stable.
If any of you have time to implement it, I may help a bit in theory part.

the link above doesn't seem to work, here's a link that works for me, however the title is slightly different...

http://www.optimization-online.org/DB_F ... 6/1675.pdf

cheers,
Antonio

[ngbinh]

Should be the same algorithm!

My link worked before.

[Erin Catto]

Hmm, this seems to require some direct solvers. It doesn't look practical for games.

[Erwin Coumans]

For small subsystems (say 6x6) I think direct methods are practical, Doom 3 SDK has very efficient SIMD optimized implementations available.

I think the paper has two ideas to improve the regular PGS convergence/performance:

• keep two sets: an active set and a non-active set (of constraints that are already satisfied)
• only solve the active set using a direct method (isn't that a effictively using a block-sparse matrix?)

I would guess that a direct method for solving a 6x6 problem converges faster and is more efficient then using a 6x6 iterative method. Also I can see how introducing block-wise solving makes the solution converge faster. Solving larger direct subsystems could be broken up into such smaller blocks, isn't it?

At what size of n the direct method breaks down?

What do you think?
Erwin

[Erin Catto]

You can't apply a simple 6x6 block solver for the direct part because there may be kinematic loops. They are using some sort of sparse Cholesky solver on the full matrix.

In other words they need to solve the big matrix problem exactly (non-iteratively).

[Erwin Coumans]

I meant using a direct (non-iterative) MLCP solver for a 6x6 blocks, just solving the sub-problems with 6 rows at a time, rather then 1 row at a time.

You can't apply a simple 6x6 block solver for the direct part because there may be kinematic loops.

Not sure if I understand, can you explain this further in detail, not too math heavy
Thanks,
Erwin

[Dirk Gregorius]

I think Erwin only wants to solve pairwise (binary) constraints as a block, e.g. 5x5 for a hinge or 3x3 for ball-socket joint. This is indeed possible as Erin shows in Box2D for example in the ball-socket joint. We are also discussing this in the LCP thread, so you might want to create a new thread from this two

Solving chains and trees of binary constraints is possible in linear time as e.g. shown by Baraff. Erin derived such a solver in the context of SI and I implemented it for a test case. This works great and gives very stiff joints. This is another advantage of SI since it doubt that this could be easily added to the ODE PGS Quickstep solver. At least I don't see an efficient implementation of this ad hoc.

Of course you can also just solve the (not necessarily sparse) linear system J*W*JT * lambda = -J * v. For kinematic loops such a system could become singular though.

[Erin Catto]

ODE has a direct solver based on Dantzig. Inside the Dantzig solver it has to do some LU factorizations of a big matrix. The matrix has dimensions equal to the number of active constraints. It would be nice if we could factor that big matrix by just solving small blocks, but this is not possible in general.

Here's a basic problem in linear algebra (not math heavy):

Code: Select all

A * x = b

We want to solve for x and A is n-by-n. In general we cannot get an exact solution for x by solving blocks. In our case A consists of a bunch of blocks corresponding to individual constraints, but there may be no pattern for solving all those blocks simultaneously without factoring A as a whole.

If we could solve any mechanical system exactly by just using blocks and solving in order-n time then I would dump SI immediately. It is just not possible.

[Erwin Coumans]

Well, usually we solver just a single row a time, so solving say 6 rows at a time (6x6 blocks) would be an improvement. Blockwise solving would be easily implemented in both SI and PGS. For Russel's quickstep, just add a type and a switch statement. Not very elegant, but trivial.

I have been working with one of the authors of the paper, so I'll ask them to make the source code available.
Hope this helps,
Erwin

[Antonio Martini]

Well, usually we solver just a single row a time, so solving say 6 rows at a time (6x6 blocks) would be an improvement. Blockwise solving would be easily implemented in both SI and PGS.

i have not read the paper, however hasn't block solving already been used for a very long time now? block solvers are intuitive and easy to implement. For example isn't Bullet solving for 3 contact normal forces in the manifold using a small direct 3x3 LCP solver?

cheers,
Antonio

[Erin Catto]

Are we discussing this paper: "An Algorithm for the Fast Solution of Linear Complementarity Problems"?

The paper does not talk about using a simple block solver.

The same linear algebra solver, namely PARDISO, is employed to solve the linear systems arising in Algorithm 3.1 and the interior-point method.

They want an exact big matrix solver to overcome the mass ratio problem. They offer no improvement to the PGS algorithm itself. They are combining PGS with a big matrix solver.

PARDISO is a complex sparse matrix solver which I doubt would fit on an SPU. I again assert that this new algorithm is not useful for games.

[Dirk Gregorius]

Well, usually we solver just a single row a time, so solving say 6 rows at a time (6x6 blocks) would be an improvement

Well, you emphasize the local block, but this doesn't necessarily improve the convergence of the whole problem. Theoretically you can get an improvement up to sqrt( 2 ) in the convergence rate, but you only get speed improvements if you efficiently solve the sub-blocks. This is e.g. the case for for a 3x3 system. Solving a 5x5 is a different story. Then there is the idea of solving big blocks like chains or trees. This is a totally different story, but the general idea still applies here as well.

For example isn't Bullet solving for 3 contact normal forces in the manifold using a small direct 3x3 LCP solver?

No, but there is some code that tries to solve 2 constraints simultaneously. I am not sure whether this is used in praxis though.

Are we discussing this paper: "An Algorithm for the Fast Solution of Linear Complementarity Problems"?

I think not. Looks more like we are discussing stiff subsolvers

[Erwin Coumans]

Let's brainstorm a bit, and keep this discussion a bit more general.

They are combining PGS with a sub solver (big matrix solver). Mixing (sub) solvers can be very useful, several games use a combination of Featherstone, small direct solvers and PGS.

I agree with Erin that introducing a slow solver such as PARADISO (or Dantzig/Lemke) doesn't sound appealing, but mixing (sub) solvers definately does.

The Dynamo physics engine had an interesting mechanism that allowed dynamic (automatic at run-time) switching between constraint solving methods. This was based on convergence if I remember correctly, but the switch could be based on arbitrary other rules (a fallback such as running out of memory on SPU etc).

As far as Bullet goes, it should be a good playground to experiment with the various ideas in 3D, just like Box2D (either using SI or quickstep). Hopefully someone gets some time to experiment with this PGS-stiff subsolvers combo, but there are still other improvements from Box2D that I'd like to see ported over to Bullet first (NGS, split impulse etc).

[KenB]

Yes, hybrid solvers is the way to go, and PGS/Direct blocking is very efficient
for many problems. I suppose this paper could be of interest.
"A Parallel Block Iterative Method for Interactive Contacting Rigid Multibody
Simulations on Multicore PCs" by Claude Lacoursière, Umeå university.

http://www.hpc2n.umu.se/para06/papers/paper_270.pdf

The full paper is in the PARA06 proceedings, but is also covered in Claude's thesis.