How to handle large mass ratios ?

[Erwin Coumans]

Iterative solvers are very popular nowadays. However most methods don't copy well with large mass ratios. How do people deal with this ?

Perhaps an LCP version of conjugate gradient may be the answer to this. Someone on the list mentioned PCG which may be what I'm looking for. Any good references on that?

Also, diagonalization was mentioned, has any tried this ?

[kenny]

Iterative solvers are very popular nowadays. However most methods don't copy well with large mass ratios. How do people deal with this ?

Well, there is a bunch of iterative methods, I guess you refer to matrix splitting methods? I will consider only these in the following discussion.
If you rephrase the iterative method as a matrix equation, for GS this is something like

x = T x + c

with

T = (D+L)^{-1} U
c = (D+L)^{-1} b

In order for GS to converge at all the spectral radius of T must be less than 1. General speaking the smaller spectral radius the better it is (or should I say the faster the convergence).

Large mass ratios seem to affect the magnitude of the eigenvalues of the system matrix. Thus the spectral radius of the A-matrix is affected.
So the trick is to get a small spectral radius! However, this is not easy.

General speaking preconditioning and diagonal scaling (diagonal scaling can be seen as a left-right preconditioning) can be used to ``change'' the spectrum of the eigenvalues. General speaking, the new preconditioned system looks like

A' = (C^{-1} A C^{-T}) C^{ T}
b' = C^{-1} b

I have tried using C = sqrt( 1.0 / diag(A) ) (this is diagonal scaling). If one compares a plot of the residual error of the new system with the old system then one see a small displacement of the curves. Diagonal scaling seems to lower the convergence rate by a constant but it does not change the asymptotic behavior.

To really get this preconditioning to work one needs to find a C-matrix which is a good approximation to A^{-1}) If this can be achieved one would have A' \approx I. Using diagonal scaling C is of course not a really good approximation to the inverse of A. Diagonal scaling seems to work best for diagonal dominant systems. In my experience multibody dynamics problems does not appear to be that. In genral matrix splitting methods seem to work best on diagonal dominant systems. Just consider the contrived case of using Jacobi on a diagonal system. Here T will be the identity matrix, and you will hit the solution in one iteration.

Shock-propagation can also be understood as a two stage blocked iterative scheme with a kind of preconditioner. The outer stage consist of deviding your problem into stack-layers. Each stack layer corresponds to a single block. Then one iterates over the blocks in a bottom-up fastion (this is a strange mixture of Jacobi/Gauss-Seidel like iteration. The diagonal blocking is similar to a Jacobi scheme, however layers must be processed in a certain order, which makes the processing of blocks sequential in nature, i.e. similar to Gauss-Seidel). For each block the mass ratios are changed by fixating some of the bodies, before solving the block.

Relaxation or damping also seem to work rather well. That is use

A' = A + epsilon I

For non-negative epsilon. Now as n-> infty we have epsilon -> zero. This is termed constaint force mixing in the ODE community. To really see what it means consider the equivalent QP. That is minimize

f(lambda) = lambda^T b + lambda^T A lambda + epsilon lambda^T lambda

s.t

lambda >= 0

Using KKT-conditions leads to an LCP with the A' matrix above. Looking closely at f one see that the magnitude of lambda is minimized by the last term. Also note lambda always will be dampen by this scheme. This leads to weaker constraint forces..... This can be used to model soft constraints but if you really want hard-constraints then you have lost!!!

However in my experiece relaxation seem to really change the slope of the convergenece rate. That is to say for the convergence rate p given by

lim_{n \rightarrow \infty} \frac{|e^n|}{|e^{n-1}|^p} = c

where c is called the convergence constant. Larger epsilon usually in my experience means lower c. Try looking at log-plots of the residual error and the different slopes becomes apparent. For GS type of solvers p=1 (this is really bad:-()

[Antonio Martini]

has anyone tried the method proposed in the following paper?:

http://www-evasion.imag.fr/Publications ... 3/jnrr.pdf

it seems to be a sort of projected biconjugate gradient solver. So it should have fewer requirements than PGS and PCG in order to converge.

Antonio

[Erin Catto]

It is odd that they are using biconjugate gradient when conjugate gradient would do. They appear to be restarting the BCG whenever the active set changes. This might be inefficient. A good PCG algorithm should work on the whole system, not just the active set. Then the PCG does not restart and needs less bookkeeping.

This paper presents a PCG method but does not go into much detail:

http://www.inrialpes.fr/bipop/publicati ... bd2005.pdf

[Erwin Coumans]

I already asked Stephane Redon, who is a collegue of 3 of the authors, for a reply from one of them. He mailed them, so lets cross the fingers

[Antonio Martini]

It is odd that they are using biconjugate gradient when conjugate gradient would do. They appear to be restarting the BCG whenever the active set changes. This might be inefficient.

they mention a Baraff 94's paper in regard to the LCP formulation and
in the initial part of section 6(dynamic friction) of the following mentioned paper:

http://citeseer.csail.mit.edu/baraff94fast.html

it can be seen that the friction model employed leads to an unsymmetric LCP matrix and a PCG solver would not work in that case. However i personally prefer the less accurate friction models that lead to a symmetric system.

talking about performances they say that their method is slightly faster than a standard BCG applied to equality constraints. A standard BCG should be twice as slower than a CG:

from:
http://mathworld.wolfram.com/Biconjugat ... ethod.html

"For symmetric positive definite systems, the method delivers the same results as the conjugate gradient method, but at twice the cost per iteration"

however it would be very interesting to see the paper's authors direct comments.

[kenny]

they mention a Baraff 94's paper in regard to the LCP formulation and
in the initial part of section 6(dynamic friction) of the following mentioned paper:

http://citeseer.csail.mit.edu/baraff94fast.html

it can be seen that the friction model employed leads to an unsymmetric LCP matrix and a PCG solver would not work in that case. However i personally prefer the less accurate friction models that lead to a symmetric system.

In this paper Baraff use a splitting method for handling friction in 3D. Strictly speaking this is in fact a NCP not a LCP.

[butcher]

I feel ashamed for admitting this, but we sidestepped the whole problem of mass ratios in Halo 2. Havok allows you to mess with the collision response directly, so we simply pin the relative masses of objects so that there is never more than a 16:1 ratio between them.