Baraff/Witkin's integration matrix not positive: what to do?
Baraff/Witkin's integration matrix not positive: what to do?
I implemented the algorithm described by Baraff and Witkin in "Large Steps in Cloth Simulation" (SIGGRAPH 98) and encountered severe instabilities with large stiffness and damping values. I found out the reason, namely that the system matrix used by the implicit integration is NOT positive definite for such values, and thus the conjugate gradient method is not applicable.
For example, consider the case of a springlike energy E = k/2 (x  d)^2. The second derivative is
d2E/dx2 = k (1  d/x) + k (d/x^3) x x^T
Now if the current extension x is less than the rest length d, take a direction y perpendicular to x. Then
y^T d2E/dx2 y = k (1  d/x) < 0,
so d2E/dx2 is not positive definite (the function E is not convex).
For timestep dt, the system matrix is M + dt^2 * d2E/dx2. For large dt, this matrix loses its positive definiteness due to the d2E/dx2contribution. This is in contradiction to the claim of Baraff/Witkin that the matrix is positive definite (which they didn't proof).
Now my question: Does anyone know of a paper where this issue (handling a nonpositive matrix) is addressed? I could reduce the time step, but that would foil the advantages of implicit integration. Or I could use normal equations (A^T A instead of A as matrix), but then I need a lot more iterations in the CG solver.
For example, consider the case of a springlike energy E = k/2 (x  d)^2. The second derivative is
d2E/dx2 = k (1  d/x) + k (d/x^3) x x^T
Now if the current extension x is less than the rest length d, take a direction y perpendicular to x. Then
y^T d2E/dx2 y = k (1  d/x) < 0,
so d2E/dx2 is not positive definite (the function E is not convex).
For timestep dt, the system matrix is M + dt^2 * d2E/dx2. For large dt, this matrix loses its positive definiteness due to the d2E/dx2contribution. This is in contradiction to the claim of Baraff/Witkin that the matrix is positive definite (which they didn't proof).
Now my question: Does anyone know of a paper where this issue (handling a nonpositive matrix) is addressed? I could reduce the time step, but that would foil the advantages of implicit integration. Or I could use normal equations (A^T A instead of A as matrix), but then I need a lot more iterations in the CG solver.

 Posts: 875
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Re: Baraff/Witkin's integration matrix not positive: what to do?
Are you sure you use the correct gradients? I lately checked the paper of DiMacri and from a quick scan I think they are not correct. Basically I took another root. Instead of building the gradients directly I built the time derivatives and found the Jacobian (gradient) and Hessian (second derivativess) by inspection. Basically:
dC/dt = del_C/del_x * dx/dt + del_C / del_t
The constraints don't depent explicitely on time so we can simplify:
dC/dt = J * v
Now we build the second time derivative:
d2C/dt^2 = dJ/dt * v + J * dv/dt
The problem is the time derivative of the Jacobian which is a rank three tensor (actually an arry of Hessian matrices). From there you can find the Hessian by inspection.
Cheers,
Dirk
PS: del_ should indicate partial derivatives
dC/dt = del_C/del_x * dx/dt + del_C / del_t
The constraints don't depent explicitely on time so we can simplify:
dC/dt = J * v
Now we build the second time derivative:
d2C/dt^2 = dJ/dt * v + J * dv/dt
The problem is the time derivative of the Jacobian which is a rank three tensor (actually an arry of Hessian matrices). From there you can find the Hessian by inspection.
Cheers,
Dirk
PS: del_ should indicate partial derivatives

 Posts: 875
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Re: Baraff/Witkin's integration matrix not positive: what to do?
I quickly checked you derivative for the spring. For my derivations I used the following identities:
C(x) = x => del_C/del_x = x / x
C(x) = x / x => del_C/del_x = 1/x * ( E3  x*x^T / (x^T*x) )
Here E3 is a 3x3 identity matrix.
Now for the energy function:
E = k/2 * ( x  L )^2
del_E/del_x = k * ( x  L ) * x / x
Now I use: del_( phi * F ) / del_x = phi * del_F / del_x + F * ( del_phi / del_x )^T where phi = ( x  L ) and F = x / x
del_2E/del_x2 = k * ( x  L ) * 1/x * ( E3  x*x^T / (x^T*x) ) + ( x / x ) * ( x / x )^T
Note that ( x / x ) * ( x / x )^T = x*x^T / (x^T*x) using x^2 = x^2 = x^T * x
del_2E/del_x2 = k * ( x  L ) * 1/x * ( E3  x*x^T / (x^T*x) ) + x*x^T / (x^T*x)
del_2E/del_x2 = k * ( x  L ) * E3 + ( 1  k * ( x  L ) ) * x*x^T / (x^T*x)
HTH,
Dirk
C(x) = x => del_C/del_x = x / x
C(x) = x / x => del_C/del_x = 1/x * ( E3  x*x^T / (x^T*x) )
Here E3 is a 3x3 identity matrix.
Now for the energy function:
E = k/2 * ( x  L )^2
del_E/del_x = k * ( x  L ) * x / x
Now I use: del_( phi * F ) / del_x = phi * del_F / del_x + F * ( del_phi / del_x )^T where phi = ( x  L ) and F = x / x
del_2E/del_x2 = k * ( x  L ) * 1/x * ( E3  x*x^T / (x^T*x) ) + ( x / x ) * ( x / x )^T
Note that ( x / x ) * ( x / x )^T = x*x^T / (x^T*x) using x^2 = x^2 = x^T * x
del_2E/del_x2 = k * ( x  L ) * 1/x * ( E3  x*x^T / (x^T*x) ) + x*x^T / (x^T*x)
del_2E/del_x2 = k * ( x  L ) * E3 + ( 1  k * ( x  L ) ) * x*x^T / (x^T*x)
HTH,
Dirk
Re: Baraff/Witkin's integration matrix not positive: what to do?
Thanks for the reply, Dirk. In your second post, the derivation is wrong (check the physical units). I also get
del_E/del_x = k * ( x  L ) * x / x
Then write this as:
del_E/del_x = k * ( 1  L/x ) * x
With del_(1/x)/del_x = 1/x^3 * x I get
del_2E/del_x2 = k * ( 1  L/x ) * E3 + (k * L / x^3) x*x^T
Now the point is that Baraff/Witkin construct the system matrix M + h^2 del_2E/del_x2 with time step h and mass matrix M, and they claim that this is positive definite. Positive definiteness is also required for the CG solver. But if you look at del_2E/del_x2, you see that the factor of E3 becomes negative if x is less than the rest length L, and then the matrix del_2E/del_x2 is no longer positive definite: You can choose a direction y perpendicular to the current difference vector x, and then y^T del_2E/del_x2 y = k * ( 1  L/x ) * y^T y < 0! This means that for large time steps h, the whole system matrix is not positive definite.
del_2E/del_x2 is positive definite for simpler energies like k/2 (xc)^2 with a "rest vector" c. Such an energy is typically used in 1dimensional models, but not for 3 dimensions, where we want to use a rest length and thus have to compute the modulus.
I am currently trying to split the system matrix into a positive and a negative part, and then to estimate a maximum time step where the matrix is still positive definite.
del_E/del_x = k * ( x  L ) * x / x
Then write this as:
del_E/del_x = k * ( 1  L/x ) * x
With del_(1/x)/del_x = 1/x^3 * x I get
del_2E/del_x2 = k * ( 1  L/x ) * E3 + (k * L / x^3) x*x^T
Now the point is that Baraff/Witkin construct the system matrix M + h^2 del_2E/del_x2 with time step h and mass matrix M, and they claim that this is positive definite. Positive definiteness is also required for the CG solver. But if you look at del_2E/del_x2, you see that the factor of E3 becomes negative if x is less than the rest length L, and then the matrix del_2E/del_x2 is no longer positive definite: You can choose a direction y perpendicular to the current difference vector x, and then y^T del_2E/del_x2 y = k * ( 1  L/x ) * y^T y < 0! This means that for large time steps h, the whole system matrix is not positive definite.
del_2E/del_x2 is positive definite for simpler energies like k/2 (xc)^2 with a "rest vector" c. Such an energy is typically used in 1dimensional models, but not for 3 dimensions, where we want to use a rest length and thus have to compute the modulus.
I am currently trying to split the system matrix into a positive and a negative part, and then to estimate a maximum time step where the matrix is still positive definite.

 Posts: 875
 Joined: Sun Jul 03, 2005 4:06 pm
 Location: Kirkland, WA
Re: Baraff/Witkin's integration matrix not positive: what to do?
You are right. I just assumed that the mistake was in the gradient because if you read the cloth papers people always state that the matrix is PD.
Viel Erfolg!
Dirk
Viel Erfolg!
Dirk
Re: Baraff/Witkin's integration matrix not positive: what to do?
Yes, this is also my impression. These papers seem to be a bit optimistic.Dirk Gregorius wrote:if you read the cloth papers people always state that the matrix is PD.
The case of nonPD is not academic, it occurs already for typical time steps together with stiff cloth. For positive definiteness, you get the condition m + h^2 k (1  L/x) > 0, so for x < L:
h^2 < m / (k * (L/x  1))
This is the typical relationship between stiffness and time step which is already known from explicit euler integration, but with an additional division by (L/x  1). Now one could argue that for stiff springs, x is always close to L, so the division is by a small number and large time steps are possible. But in a test with damped stiff structural and bending springs (bending springs being based on angles), I had to use time steps of about 1/1000 second.
Re: Baraff/Witkin's integration matrix not positive: what to do?
Hi:
I roughly read the paper, but it seems baraff only mentioned that the K matrix(d2E/dx2 )
is a symmetric Matrix not a PD matrix (cha 4.1).
and in 6.1 he mentioned that CG methods, however, require symmetric matrices and there is also a label "In fact, they work best on positive definite symmetric matrices", I don't knwon CG Method so I suppose that CG should work not only with a PD matrix but symmetric is enough?
I roughly read the paper, but it seems baraff only mentioned that the K matrix(d2E/dx2 )
is a symmetric Matrix not a PD matrix (cha 4.1).
and in 6.1 he mentioned that CG methods, however, require symmetric matrices and there is also a label "In fact, they work best on positive definite symmetric matrices", I don't knwon CG Method so I suppose that CG should work not only with a PD matrix but symmetric is enough?
Re: Baraff/Witkin's integration matrix not positive: what to do?
They state: "The matrices we ultimately hand to our CG method are positive definite." and "Let us define the symmetric positive definite matrix A by A = M  h df/dv  h^2 df/dx."
If you read mathematical papers on the CG method, there is the requirement of positivity. Maybe that you are lucky and CG converges also for some special situations with indefinite matrices, but that's not guaranteed. Theoretically, if A is indefinite for large h, you can also find a h where A is not invertible (the h where A switches from positive definite to indefinite). But with CG we want to solve Ax=b, so we need an invertible A.
If you read mathematical papers on the CG method, there is the requirement of positivity. Maybe that you are lucky and CG converges also for some special situations with indefinite matrices, but that's not guaranteed. Theoretically, if A is indefinite for large h, you can also find a h where A is not invertible (the h where A switches from positive definite to indefinite). But with CG we want to solve Ax=b, so we need an invertible A.
Re: Baraff/Witkin's integration matrix not positive: what to do?
Is it possible that they don't use the springlike energy . I remeber they only use stretch ,bend and shear energy. May be these energy will cause a PD matrix but not all energy you derive will be a PD matrix?
Re: Baraff/Witkin's integration matrix not positive: what to do?
Yes, they only use stretch/bend/shearterms, but the stretchterm is a spring. They have a more general setting with triangles and uvcoordinates so that the behaviour is not so dependent on the tesselation, but for a triangle with one u and one vparallel edge, the formula for the stretchterm is exactly that of a spring (or actually two springs, one for the uedge and one for the vedge, but that doesn't matter).
Re: Baraff/Witkin's integration matrix not positive: what to do?
Hi mathis :
Have you ever tested to derive the derivation from Baraff's stretch enegry directly but not use your spring form? Because I think there is still a little difference beacasue baraff use the particle's position as varaible rather than the extension of the spring , could this induce the K matrix to be different of yours,just my guess .
Have you ever tested to derive the derivation from Baraff's stretch enegry directly but not use your spring form? Because I think there is still a little difference beacasue baraff use the particle's position as varaible rather than the extension of the spring , could this induce the K matrix to be different of yours,just my guess .
Re: Baraff/Witkin's integration matrix not positive: what to do?
The formulas are the same for a uvaligned triangle. Let (x0, x1, x2) be such a triangle with the uvcoordinates (0, 0), (L1, 0), (0, L2). L1 and L2 are the rest lengths of the edges x0x1 and x0x2, respectively. This fulfills the requirement of the paper that dw/du = dw/dv = 1 if w is the rest state w(u,v) = (u,v,0), where w in general maps (u,v) into the current 3d vertex position.
Now for the current state w, the paper approximates
(dw/du, dw/dv) = (x1x0, x2x0) * ((u1u0, u2u0), (v1v0, v2v0))^1
In our case, the matrix ((u1u0, u2u0), (v1v0, v2v0)) is a diagonal with the entries L1, L2, so the formula simplifies to
dw/du = (x1  x0) / L1
dw/dv = (x2  x0) / L2
Now the paper uses C(x0,x1,x2) = a * (dw/du  bu, dw/dv  bv) as constraint and k/2 C^T C as energy. This leads to the energy
E = a^2 k / 2 * ( ((x1  x0)/L1  bu)^2 + ((x2  x0)/L2  bv)^2 )
which is the sum of the energy of two springs with rest lengths bu*L1, bv*L2 and spring constants a^2 k / L1^2, a^2 k / L2^2. And this energy leads to a nonPD system matrix for large time steps and extensions less than the rest lengths.
Now for the current state w, the paper approximates
(dw/du, dw/dv) = (x1x0, x2x0) * ((u1u0, u2u0), (v1v0, v2v0))^1
In our case, the matrix ((u1u0, u2u0), (v1v0, v2v0)) is a diagonal with the entries L1, L2, so the formula simplifies to
dw/du = (x1  x0) / L1
dw/dv = (x2  x0) / L2
Now the paper uses C(x0,x1,x2) = a * (dw/du  bu, dw/dv  bv) as constraint and k/2 C^T C as energy. This leads to the energy
E = a^2 k / 2 * ( ((x1  x0)/L1  bu)^2 + ((x2  x0)/L2  bv)^2 )
which is the sum of the energy of two springs with rest lengths bu*L1, bv*L2 and spring constants a^2 k / L1^2, a^2 k / L2^2. And this energy leads to a nonPD system matrix for large time steps and extensions less than the rest lengths.
Re: Baraff/Witkin's integration matrix not positive: what to do?
Oh, I now realize they are the same, so this mean the matrix in baraff 's original paper of stretch enegy is not a PD matrix ..strange
Re: Baraff/Witkin's integration matrix not positive: what to do?
Well, it is not PD only for large time steps and large stiffness/damping constants, but as their paper is entitled "large steps in cloth simulation", they should have realized this
I realized it already with my second test of my implementation where the simulation just exploded. After some hours of further testing, I discovered the nonPDness of the matrix.
I realized it already with my second test of my implementation where the simulation just exploded. After some hours of further testing, I discovered the nonPDness of the matrix.
 Erwin Coumans
 Site Admin
 Posts: 4235
 Joined: Sun Jun 26, 2005 6:43 pm
 Location: California, USA
 Contact:
Re: Baraff/Witkin's integration matrix not positive: what to do?
Finally, it took 11 years but it was worth waiting: a new paper by Theodore Kim addresses the problem, and one of the references points to this forum thread. See http://www.tkim.graphics/FEMBW or https://twitter.com/_TheodoreKim/status ... 0997393409