Equivalence sequential impulses & Projected Gauss Seidel

 Posts: 324
 Joined: Fri Jul 01, 2005 5:29 am
 Location: Irvine
 Contact:
Equivalence sequential impulses & Projected Gauss Seidel
So ODE got the PGS optimizations from Croteam. I wonder where Croteam got it from.
If you study the mathematics of PGS, you can find that they are equivalent to the sequential impulse algorithm of Mirtich. In the sequential impulse approach, the sparsity and auxilliary variables are implicit. The impulse method has low/high limits for friction so that's covered too.
If you study the mathematics of PGS, you can find that they are equivalent to the sequential impulse algorithm of Mirtich. In the sequential impulse approach, the sparsity and auxilliary variables are implicit. The impulse method has low/high limits for friction so that's covered too.

 Posts: 135
 Joined: Wed Jul 27, 2005 10:28 am
 Location: SCEE London
please correct me if im wrong, in PGS the contact constraints are solved at the same istant of time, in Mirtich after an impulse has been applied, bodies are numerically integrated, collision detection is performed again etc. So the system is stepped in time while solved. Are you may be referring to the propagation method mentioned in Baraff's and in the more recent Kenny Erleben's book or in code the Bullet solver?Erin Catto wrote: If you study the mathematics of PGS, you can find that they are equivalent to the sequential impulse algorithm of Mirtich. In the sequential impulse approach, the sparsity and auxilliary variables are implicit. The impulse method has low/high limits for friction so that's covered too.
(see also section 8 )
http://www.cs.cmu.edu/~baraff/papers/sig89.pdf
if they are the same it would be interesting to see whats the equivalent of general preconditioning, warm starting and overrelaxation.
Last edited by Antonio Martini on Thu Jan 05, 2006 4:44 pm, edited 1 time in total.

 Posts: 324
 Joined: Fri Jul 01, 2005 5:29 am
 Location: Irvine
 Contact:
Right, I forgot that Mirtich moves time forward after each impulse.
But Nonconvex Rigid Bodies with Stacking http://graphics.stanford.edu/papers/rigid_bodiessig03/ uses sequential impulses which is equivalent to PGS, and that paper is from 2003.
But Nonconvex Rigid Bodies with Stacking http://graphics.stanford.edu/papers/rigid_bodiessig03/ uses sequential impulses which is equivalent to PGS, and that paper is from 2003.
There is a subtle difference in the approaches... It dependens on how the righthand side vector b in the complementarity formulation is set up.
I believe Projected Gauss Seidel as done by Moreau and Jeans Sweeping Process (blocked nonsmooth gauss seidel) is exactly the same as the ``propagation method'', this is because the righthand vector in the subproblem is reevaluated each time... this is the same as would happen in a propagation method.
However, in the nonblocked approach used by nearly everybody, one computes the bvector at a single instant in time. Therefore one can observe a difference in simulation results, if you don't believe me take a look at these comparisons that I did sometime ago
http://www.diku.dk/~kenny/visionday05.pdf
Of course the Guendelmann et al. approach is not allowed to run until convergence, but neither is the PGS solver:)
/Kenny
I believe Projected Gauss Seidel as done by Moreau and Jeans Sweeping Process (blocked nonsmooth gauss seidel) is exactly the same as the ``propagation method'', this is because the righthand vector in the subproblem is reevaluated each time... this is the same as would happen in a propagation method.
However, in the nonblocked approach used by nearly everybody, one computes the bvector at a single instant in time. Therefore one can observe a difference in simulation results, if you don't believe me take a look at these comparisons that I did sometime ago
http://www.diku.dk/~kenny/visionday05.pdf
Of course the Guendelmann et al. approach is not allowed to run until convergence, but neither is the PGS solver:)
/Kenny

 Posts: 874
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Guendelman  Mirtich:
The major difference  at least in my understanding  is that Guendelman doesn't need any velocity threshold since he uses the new integration scheme:
 Candidates
 Collision
 Velocity Update
 Contact
 Position Update
What I also like is, that with this method you don't need the exact TOI.
Impulse Based Simulation vs. Velocity Constraints (LCP  here PGS)
There is a big difference between them and I don't agree here. In the impulse based formulation you seek a sequence of impulses such that the relative normal velocity at each contact point becomes zero. Here you solve one contact point after the other reevolving the velocity after each contact point as opposed to the PGS where you solve simultan. The iterative nature of the PGS makes them similar, but the idea is totally different.
Bullet  Guendelman
A very good idea of Guendelman (and probably a difference to Bullet) is that instead of applying the full impulse at each contact point he uses negative coefficients of restitution which physically is equivalent to slowing down the object instead of suddenly stop it. It also gives you a better sampling. What also makes a difference is using the current position for the impulse, but still using candidates for finding contacts. I also apply the impulses at the contact with highest penetrating velocity first, what is the same to Guendelman who applies the impulse at the contact with the deepest penetration. Since I don't refresh my contact information I use this heuristic. See:
K. Erleben book: Chapter 5 "Impulse Based Dynamics"
Chatterjee: http://ruina.tam.cornell.edu/research/t ... d_body.pdf
Guendelman: PowerPoint Contact processing ( Improving accuracy or so )
On average impulses:
In order to avoid the awful tilting which happens through sequential applying impulses (s. Kenny) I use an averaged impulse in the collision phase. Instead of using the arithmetic mean of the contact points I use a weighted sum scaled by the relative normal velocity  only considering points with penetrating velocity. John Schultz at GameDev.net uses something similar. Instead of using the relative normal velocity he scales by the normal impulse (length) at each contact point.
Both schemes give good result for well defined stacking, but have problems for twisted and offset stacks. Also note that you totally miss the friction when the upper box only rotates around the contact normal. The averaging will cancel out the "Torsion Momentum". This is work in progress here...
@Kenny:
(1) Novodex uses Block LCP (at least for the joints). I have a paper of the ETHZ where Adam is so nice to explain quite some details. IIRC they have removed the paper, so if you can't find it I can send it Erwin and he maybe makes it accessable to a bigger audience...
(2) When relaxing the contacts in the Guendelman approach instead of correcting the whole relative velocity at each contact point directly Guendelman starts with smaller impulses (s. above). Do you think that this could be carried over to the LCP formulation.
It can be shown that a SOR paramater must lie between 0 < w_sor < 2 for convergence:
 How do I coose this parameter ? The ODE uses 1.3  why?
 Wouldn't make underrelaxation make sense here?
 What about making the relaxation parameter be a function of the iteration count, e.g. w1 = 1.0, w2 = 1.0, w3 = 0.9, ... Or maybe ascending, e.g. w1 = 0.1, w2 = 0.2, w3 = 0.3 . Is there something like "adaptive" SUR/SOR around? I think of two opposing constraints here, where once cancels out the other. When you use smaller w parameter later in the solver, you maybe could smear the error a bit this way...?
Regards,
Dirk
HTH,
Dirk
The major difference  at least in my understanding  is that Guendelman doesn't need any velocity threshold since he uses the new integration scheme:
 Candidates
 Collision
 Velocity Update
 Contact
 Position Update
What I also like is, that with this method you don't need the exact TOI.
Impulse Based Simulation vs. Velocity Constraints (LCP  here PGS)
There is a big difference between them and I don't agree here. In the impulse based formulation you seek a sequence of impulses such that the relative normal velocity at each contact point becomes zero. Here you solve one contact point after the other reevolving the velocity after each contact point as opposed to the PGS where you solve simultan. The iterative nature of the PGS makes them similar, but the idea is totally different.
Bullet  Guendelman
A very good idea of Guendelman (and probably a difference to Bullet) is that instead of applying the full impulse at each contact point he uses negative coefficients of restitution which physically is equivalent to slowing down the object instead of suddenly stop it. It also gives you a better sampling. What also makes a difference is using the current position for the impulse, but still using candidates for finding contacts. I also apply the impulses at the contact with highest penetrating velocity first, what is the same to Guendelman who applies the impulse at the contact with the deepest penetration. Since I don't refresh my contact information I use this heuristic. See:
K. Erleben book: Chapter 5 "Impulse Based Dynamics"
Chatterjee: http://ruina.tam.cornell.edu/research/t ... d_body.pdf
Guendelman: PowerPoint Contact processing ( Improving accuracy or so )
On average impulses:
In order to avoid the awful tilting which happens through sequential applying impulses (s. Kenny) I use an averaged impulse in the collision phase. Instead of using the arithmetic mean of the contact points I use a weighted sum scaled by the relative normal velocity  only considering points with penetrating velocity. John Schultz at GameDev.net uses something similar. Instead of using the relative normal velocity he scales by the normal impulse (length) at each contact point.
Both schemes give good result for well defined stacking, but have problems for twisted and offset stacks. Also note that you totally miss the friction when the upper box only rotates around the contact normal. The averaging will cancel out the "Torsion Momentum". This is work in progress here...
@Kenny:
(1) Novodex uses Block LCP (at least for the joints). I have a paper of the ETHZ where Adam is so nice to explain quite some details. IIRC they have removed the paper, so if you can't find it I can send it Erwin and he maybe makes it accessable to a bigger audience...
(2) When relaxing the contacts in the Guendelman approach instead of correcting the whole relative velocity at each contact point directly Guendelman starts with smaller impulses (s. above). Do you think that this could be carried over to the LCP formulation.
It can be shown that a SOR paramater must lie between 0 < w_sor < 2 for convergence:
 How do I coose this parameter ? The ODE uses 1.3  why?
 Wouldn't make underrelaxation make sense here?
 What about making the relaxation parameter be a function of the iteration count, e.g. w1 = 1.0, w2 = 1.0, w3 = 0.9, ... Or maybe ascending, e.g. w1 = 0.1, w2 = 0.2, w3 = 0.3 . Is there something like "adaptive" SUR/SOR around? I think of two opposing constraints here, where once cancels out the other. When you use smaller w parameter later in the solver, you maybe could smear the error a bit this way...?
Regards,
Dirk
HTH,
Dirk

 Posts: 324
 Joined: Fri Jul 01, 2005 5:29 am
 Location: Irvine
 Contact:
If you study PGS, you can discover that what it does is solve one constraint at a time and update the net force on the two bodies (assuming only twobody constraints). Since PGS works at the velocity level, this is equivalent to applying an impulse and updating the current velocity of the bodies. This is the mathematical similarity.Impulse Based Simulation vs. Velocity Constraints (LCP  here PGS)
There is a big difference between them and I don't agree here. In the impulse based formulation you seek a sequence of impulses such that the relative normal velocity at each contact point becomes zero. Here you solve one contact point after the other reevolving the velocity after each contact point as opposed to the PGS where you solve simultan. The iterative nature of the PGS makes them similar, but the idea is totally different.
Now PGS and sequential impulses may use different impulses based on different friction models, position error resolution, and so forth. But the basic idea of solving one constraint at a time and propagating the result is identical.
If you consider a frictionless system with inelastic contacts, the two algorithms give identical results assuming collision and position updates are handled identically. Of course, the order of collision processing and position updates are independent of how velocities are updated.

 Posts: 874
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
I see  you mean since you update the lambda (impulse) the velocity change is implicit. Thanks for pointing this out!
@Erin and Kenny:
Did you try different relaxation methods, like SOR or SUR.
And maybe another question. Would it be possible to solve the contact situation the "classical" way, while solving the joints in a block wise fashion? Here is a quote from a Novodex presentation at the ETHZ:
Then two questions on setting up the Jacobian:
1) Would I gain an advantage when setting up the first friction direction parallel to the tangential direction of the tangential velocity at the contact point instead of choosing two arbitrary perpendicular vectors in the contact plane?
2) For long lever arms the rotation term gets predominant. Are there more clever ways to setup the joint Jacobians such that they are more PGS friendly? I mean are there some consideration on choosing the basis for the joint constraint? When one body is fixed I see that I should choose one basis vector in the direction of the lever arm, but how should I choose the joint basis for two dynamic bodies?
Dirk
@Erin and Kenny:
Did you try different relaxation methods, like SOR or SUR.
And maybe another question. Would it be possible to solve the contact situation the "classical" way, while solving the joints in a block wise fashion? Here is a quote from a Novodex presentation at the ETHZ:
How could this be achieved?? Because of strong interdependency of
Jacobian rows, pure GS has problems.
? Use blocked GS for joints.
Then two questions on setting up the Jacobian:
1) Would I gain an advantage when setting up the first friction direction parallel to the tangential direction of the tangential velocity at the contact point instead of choosing two arbitrary perpendicular vectors in the contact plane?
2) For long lever arms the rotation term gets predominant. Are there more clever ways to setup the joint Jacobians such that they are more PGS friendly? I mean are there some consideration on choosing the basis for the joint constraint? When one body is fixed I see that I should choose one basis vector in the direction of the lever arm, but how should I choose the joint basis for two dynamic bodies?
? Choice of Basis for constraint
Bad Good Good
? Rotation term predominant:
dv = n*j/m; dw = j*inv(I)*r x n, r being lever arm.
? If both objects are movable then we can?t do #2
Dirk
I am not sure if I entirely understand your question/problem?DonDickieD wrote:Did you try different relaxation methods, like SOR or SUR.
Yes, but for my purpose it was not worth the trouble.
DonDickieD wrote: And maybe another question. Would it be possible to solve the contact situation the "classical" way, while solving the joints in a block wise fashion?
/Kenny

 Posts: 324
 Joined: Fri Jul 01, 2005 5:29 am
 Location: Irvine
 Contact:
Here's a paper that makes a similar claim about the GaussSeidel nature of sequential impulses:
http://www.ep.liu.se/ecp/010/004/ecp01004.pdf
I find it interesting that with PGS it is viewed as a clever optimization to use an accumulator, but with sequential impulses the optimization is builtin.
http://www.ep.liu.se/ecp/010/004/ecp01004.pdf
I find it interesting that with PGS it is viewed as a clever optimization to use an accumulator, but with sequential impulses the optimization is builtin.

 Posts: 874
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Thanks for the link Erin. I found the paper very interessting.
I recently looked at the Doom3 SDK source and they seem to have an impulse based formulation as well as a LCP solver. I have not enough knowledge or experience with LCP solvers to recognize the kind of the solver. Do they use a GS or similar method or is this already a more advanced solver that is superior to the GS in the context of computer games?
Dirk
I recently looked at the Doom3 SDK source and they seem to have an impulse based formulation as well as a LCP solver. I have not enough knowledge or experience with LCP solvers to recognize the kind of the solver. Do they use a GS or similar method or is this already a more advanced solver that is superior to the GS in the context of computer games?
Dirk
 Erwin Coumans
 Site Admin
 Posts: 4187
 Joined: Sun Jun 26, 2005 6:43 pm
 Location: California, USA
 Contact:
A few years ago, Jan Paul van Waveren (idSoftware's Doom3 physics programmer) told me he used impulsebased method for contacts and an LCP solver for ragdolls (noncontact joints). Probably he uses a direct /pivoting LCP method.I recently looked at the Doom3 SDK source and they seem to have an impulse based formulation as well as a LCP solver. I have not enough knowledge or experience with LCP solvers to recognize the kind of the solver. Do they use a GS or similar method or is this already a more advanced solver that is superior to the GS in the context of computer games?
I'll ask if he can explain here on the forum.
It is typical to mix a direct LCP solver with iterative method. Especially high speed Ragdolls don't behave well on impact, their limbs detach. Typical scenario is a high speed car against a pedestrian, or a character standing next to a barrel (explosion). I've seen this in a game I worked on in the past, and we used Featherstone for the ragdolls, and impulse based for contacts.
Erwin

 Posts: 874
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Good idea  this would be very interessting.A few years ago, Jan Paul van Waveren (idSoftware's Doom3 physics programmer) told me he used impulsebased method for contacts and an LCP solver for ragdolls (noncontact joints). Probably he uses a direct /pivoting LCP method.
I'll ask if he can explain here on the forum.
Oh yes, we have exactly this problem. When a dragon takes a knight into his mouth and throws him around the ragdoll is always torn apartIt is typical to mix a direct LCP solver with iterative method. Especially high speed Ragdolls don't behave well on impact, their limbs detach. Typical scenario is a high speed car against a pedestrian, or a character standing next to a barrel (explosion). I've seen this in a game I worked on in the past, and we used Featherstone for the ragdolls, and impulse based for contacts.
What exactly do you mean by "mixing"? Are the nonarticulated bodies handles in an iterative solver, while the articulated bodies are solved by some kind of LCP solver (PGS or whatever).
Regarding Featherstone:
How is the collison impulse propagated through the body? While maybe not accurate, the impulse could maybe only be applied to colliding body of the hierarchy. While physical not correct, this could maybe look good enough in a game.
Here is a paper I once spotted, but never read actually  silly me . Maybe interesting in this context:
http://www.gdconf.com/conference/archiv ... ngelis.pdf
Dirk

 Posts: 324
 Joined: Fri Jul 01, 2005 5:29 am
 Location: Irvine
 Contact:
Conceptually the idea of combining iteration on contacts with joint coordinates for bilateral constraints is appealling. But I wonder how well such a system would perform.
There is nothing that says the ragdoll bone has to get full control of the animation bone. You can just apply the rotation if you like.
BTW, I have worked on such a system and the external impulse just propagates through the joint coordinates. This can be done in O(n) time using the Featherstone mass matrix factorization.
There is nothing that says the ragdoll bone has to get full control of the animation bone. You can just apply the rotation if you like.
BTW, I have worked on such a system and the external impulse just propagates through the joint coordinates. This can be done in O(n) time using the Featherstone mass matrix factorization.

 Posts: 135
 Joined: Wed Jul 27, 2005 10:28 am
 Location: SCEE London
regarding the optimized algorithm of section 4.3 of the above paper, it should be noticed that if we use the propagation method for solving the contacts the algorithm simplifies even further. We would need to know only the effect of a contact on itself informally speaking and so the last loop(but also the first) will not be from "1 to n" but rather "1 to l" if i have understood the paper's notation. We would need to traverse only the subchain starting from the body where the impulse is applied to the root node.DonDickieD wrote: http://www.gdconf.com/conference/archiv ... ngelis.pdf
Dirk
If PGS and the propagation methods are "equivalent", here the latter is a clear winner when it come to down to integrating the Featherstone method, given that PGS need to know the effect that the impulse applied at one contact point has on all the other contacts rather than just one.
However if by 'equivalent' we don't mean 'similar' i assume that the results for one method(preconditioning,overrelaxation,warm starting,contacts handling) transfer automatically to the other one.
Last edited by Antonio Martini on Mon Jan 09, 2006 7:42 pm, edited 6 times in total.
For articulated figures (like ragdolls) DOOM III uses a direct O(n) solver for treelike structures with bilateral constraints. Furthermore a direct Box Constrained Mixed LCP solver based on an incremental LDLt decomposition is used for all other constraints (including contacts). The advantages of the separate solver for treelike structures with bilateral constraints are better performance but most of all the bilateral constraints are accurately enforced without being pulled apart due to typical numerical issues of other constraints (like too many or conflicting contacts). All the matrix math is SIMD optimized and numerous tricks are used to save work in the direct solvers. The complete source code is available in the DOOM III and the Quake 4 SDK (http://www.iddevnet.com/).
For simple rigid bodies (like boxes and barrels) DOOM III uses a separate impulse based solver. This solver is rather inaccurate but also very cheap and only solves one contact per object every frame (at 60 fps). Objects sometimes feel a little sticky because so few impulses are applied over time. However, in combination with the continuous collision detection system you can still properly stack hundreds of objects without real problems.
JP,
For simple rigid bodies (like boxes and barrels) DOOM III uses a separate impulse based solver. This solver is rather inaccurate but also very cheap and only solves one contact per object every frame (at 60 fps). Objects sometimes feel a little sticky because so few impulses are applied over time. However, in combination with the continuous collision detection system you can still properly stack hundreds of objects without real problems.
JP,