Non-linear Gauss-Seidel solver

[Himura78]

Hi,
One more question from me regarding 2D physics Following the evolution of Box2D, I see that different approaches were tried in terms of position stabilization. I understand Baumgarte stabilization (feeding a fraction of the position error back into the velocity solver) and can see some potential issues with that (since the velocities are adjusted, momentum can be artificially affected). I see that the other method Box2D uses is NGS (non-linear Gauss-Seidel), but haven't found a good description of the idea behind this or how it works (I haven't seen any papers or presentations on it). I thought it might be helpful to ask about this before trying to puzzle over the code to work out how it works.
Any insight into this would be appreciated.

[Dirk Gregorius]

Box2D has a separate projection step for the position constraints. As you pointed out already in another thread velocity constraints are linear, but position constraints are not. So you need to solve a *non-linear* system equations. Conceptually you can think of an outer loop for the Newton solver and then an inner loop for the Gauss-Seidel loop. If you merge both loops you get the NGS.

An easier way to think about it is since you are now changing positions you cannot pre-compute effective mass and Jacobians as they depend on the position themselves. So for every constraint you project you need to re-build everything that depends on the involved bodies. This makes this solver pretty expensive.

Conceptually this is identical to position based dynamics, but for rigid bodies. There are a bunch of papers on this. In the end this all goes back to Shake and Rattle solvers from the 70s. E.g. see here
http://www.unige.ch/~hairer/poly-sde-mani.pdf

[Himura78]

Thanks Dirk. I suppose I'm still not clear on it... from Erin Catto's "Modeling and solving constraints" presentation, for the velocity constraints, we use:

Newton's law -> v2 = v1 + M^-1 * p
Virtual work -> p = J^T * lambda
Velocity constraint -> J * v2 = 0

and we form the linear equation

J * M^-1 * J^T * lambda = -J*v

which we solve for lambda. But what's the equivalent system that we're solving for the position constraints?

[Dirk Gregorius]

Short answer:
J * M^-1 * J^T * lambda = -C

This is the inner linear system in a Newton solver. The lambdas are now 'pushes' and get transformed into translations dx and rotations dq similar to impulses. After you applied the pushes the Jacobian and mass matrix (inertia) changes, so you need to rebuild J * M^-1 * J^T.

When you solve on the velocity level the solution will drift away from the constraint manifold. E.g. For a pendulum the mass will not swing on a cicular arc, but slowly drift away. The above approach then performs an orthogonal projection back onto the constraint manifold. For the pendulum this simply means that after you integrated the velocities and found the new position you find the closest point on the circle and put it there.

This is the easiest form of so called projection methods. Check out the paper by Hairer I linked above. Or for a simpler introduction you can look at the Jacobsen cloth paper or Mueller's Position Based Dynamics.

This is a more practical explanation. If you want a more mathematical derivation let me know and I give it a try

HTH,
-Dirk

[RandyGaul]

Wasn't post projection (NGS) sort of hacky? As in J was stolen from the velocity derivation and plugged in to form a psuedo position formula? Description of this by Cline in his nice paper: https://pdfs.semanticscholar.org/476d/f ... cbc781.pdf. In there it seems J is taken from the velocity derivation and used as G for position projection. Seems a bit like math hackery, but since this is all non-physical anyway I suppose it makes quite a bit of practical sense. Thoughts?

[Dirk Gregorius]

No, there is nothing hacky about orthogonal projection. It is simply an alternativ to Baumgart'ish stabilization to keep your solution on the constraint manifold. The linked Hairer paper discusses projection methods in general and in the context of multibody simulation.

The J is not stolen, it comes from here:

C( x + dx ) ~= C( x ) + J * dx

[Himura78]

Hi Dirk, thanks again for your reply.
The Hairer paper is a bit beyond me at this stage but I remember reading the Jacobsen paper a while back and it seemed to make sense, the thrust of it being to resolve constraints on the position level directly by projecting bodies out of penetration.
I understand your description, for the most part, though I'm not clear on the concept of the inner and outer loops, or how you arrive at the equation mentioned. I would definitely be interested in seeing a derivation if you have time

[Dirk Gregorius]

Ok, I give it a try.

We integrate the forces and velocities to find a new position p using some unstabilized method of choice. Given the new p we now want to find a correction dp such that C(p + dp) = 0. We can approximate this as shown above by

C(p + dp) ~= C(p) + J * dp = 0

We also want to restrict the correction to be orthogonal to the constraint manifold and scale the corrections relative to the masses such that

dp = M^-1 * J^T * lambda

Combining the two equations yields

J * M^-1 * J^T * lambda = -C(p)
dp = M^-1 * J^T * lambda

This is a single Newton-Raphson step for an iterative solution to a non-linear system of constraints. So we essentially get

for ( n Newton iterations )
{
J * M^-1 * J^T * lambda = -C(p)
dp = M^-1 * J^T * lambda
p += dp
}

This builds the outer iterative loop. We could now also solve the inner linear system J * M^-1 * J^T lambda = -C(p) using an iterative Gauss-Seidel solver. This is what I referred to as the inner loop earlier. NGS simply merges the two loops solving a single constraint at a time constantly updating the mass, Jacobians and constraint error to account for the non-linearity of the problem. The key observation is simply that both the inertia and Jacobian are functions of the position and need to be updated when iterating on the position level.

Of course you can make the argument now that for small errors you can keep M and J constant and only iterate on the constraint error. You might even want to reuse the Jacobian and mass from the velocity phase. This is essentially what split impulse does. My experience is identical to Erin's in Box2d that this works relatively good for contacts, but not so much for joints.

I am not a mathematician and I have more of a practical engineering view on the problem, but hopefully this still makes sense!

Cheers,
-Dirk

[Himura78]

Thanks Dirk, I think I mostly understand it now.
Just to clarify, are you saying that this

Code: Select all

For n iterations
  Update J, M and C using p
  Solve J * M^-1 * J^T * lambda = -C
  Update dp
  Update p

(Assuming p, dp, C and lambda are vectors and J and M are matrices, ie we're solving all constraints simultaneously each iteration)

Is equivalent to this

Code: Select all

For n iterations
  For m constraints
    Update J, M and C using p
    Solve lambda = -C / (J * M * J^ T)
    Update dp 
    Update p

(Assuming J is a vector and M and C are scalar, ie we're solving a single constraint per inner loop iteration)? I guess this is similar to the similarity between the matrix algorithm used in the "iterative dynamics with temporal coherence" paper vs box2d's sequential impulse?

Also from my limited understanding of Newton raphson, it doesn't always converge to a solution, how are we sure here that this won't be an issue?

[Dirk Gregorius]

I would add one step for clarity:

Code: Select all

For n Newton iterations
  Update J, M and C using p
  For m Gauss-Seidel iterations
    For k constraints
      Solve lambda = -C / (J * M * J^ T)
      Update dp 
      Update p

And then you merge it like you show. I would not say it is equivalent, but it can converge against the same solution if you iterate enough.

Yes, in the final implementation we usually solve a single scalar constraint. You can also solve several constraints together. This is called blocked solving. You have no guarantee for convergence, but things work usually fine in practice.

[Himura78]

Thanks Dirk. Really appreciate all your replies.
After about 9 years of studying this stuff in my spare time I almost feel ready to code a simulator

[Dirk Gregorius]

Sure thing! I am glad I could help.

One thing I found over the years is that the theory can be quite involved, but the final result is relatively simple. On the downside there are quite some problems in the detail and implementation. But you will figure those out along the way and there is often no general solution. E.g. a game developer can solve many problems in content, while a middleware provider has to solve most problems in code.

I think you should give it a try. The way you asked questions and how you presented them I think you will be fine. I always recommend to use Box2D Lite as an initial guideline and then iterate from there. The magic ingredient for physics engines is to try to always keep a workable solution and solve one problem at a time. Physics and collisions get really hard if you try to solve too much at a time and need to figure out what is broken. Baby steps are a good way to learn to walk.

HTH,
-Dirk