Real-Time Physics Simulation Forum

Re: Non-linear Gauss-Seidel solver

2017-05-09T02:52:11+00:00

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

Statistics: Posted by Dirk Gregorius — Tue May 09, 2017 2:52 am

Re: Non-linear Gauss-Seidel solver

2017-05-09T01:38:42+00:00

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

Statistics: Posted by Himura78 — Tue May 09, 2017 1:38 am

Re: Non-linear Gauss-Seidel solver

2017-05-08T16:45:13+00:00

I would add one step for clarity:

Code:

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.

Statistics: Posted by Dirk Gregorius — Mon May 08, 2017 4:45 pm

Re: Non-linear Gauss-Seidel solver

2017-05-08T04:29:19+00:00

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

Code:

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:

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?

Statistics: Posted by Himura78 — Mon May 08, 2017 4:29 am

Re: Non-linear Gauss-Seidel solver

2017-05-04T17:02:39+00:00

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

Statistics: Posted by Dirk Gregorius — Thu May 04, 2017 5:02 pm

Re: Non-linear Gauss-Seidel solver

2017-05-04T11:47:53+00:00

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

Statistics: Posted by Himura78 — Thu May 04, 2017 11:47 am

Re: Non-linear Gauss-Seidel solver

2017-05-04T03:52:00+00:00

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

Statistics: Posted by Dirk Gregorius — Thu May 04, 2017 3:52 am

Re: Non-linear Gauss-Seidel solver

2017-05-04T02:04:12+00:00

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?

Statistics: Posted by RandyGaul — Thu May 04, 2017 2:04 am

Re: Non-linear Gauss-Seidel solver

2017-05-03T15:00:17+00:00

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

Statistics: Posted by Dirk Gregorius — Wed May 03, 2017 3:00 pm

Re: Non-linear Gauss-Seidel solver

2017-05-03T12:23:40+00:00

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?

Statistics: Posted by Himura78 — Wed May 03, 2017 12:23 pm