Physics Simulation Forum

hi,
the short version of the following is:
given the function f(x) = x/|x|, where x is a vector, what is the partial derivative matrix of f wrt the components of x?


the long version:

I've been trying to figure out the "Position Based Dynamics" paper from a few months ago: http://www.continuousphysics.com/Bullet/phpBB2/viewtopic.php?t=577

Dirk Gregorius has kindly helped me via email, and suggested that I post here.

I've only taken a single Calculus course, and that only dealt with scalar-valued functions of a single variable.. so, I'm slightly out of my league. I think I understand how to build a Jacobian/partial derivative matrix of a vector-valued function in theory -- for instance, I understand the Jacobian-based IK method. But in that case, the function which maps joint coordinates to effector coordinates is simple.. I'm having trouble moving beyond that.


The paper requires that I know the gradient of my constraint functions wrt the function inputs. I'm trying things out in 2D, which causes a problem: unless I'm mistaken, many of the gradients given in the paper cease to make sense when moved into 2D naively (i.e by setting the z coords to 0).


For instance: the paper describes a constraint that maintains an angle between two triangles; the angle is measured as the arccos of the dotproduct of two unit direction vectors: the triangle normals.

The normalized cross product (p0 x p1)/|p0 x p1| is used to describe the direction vectors, which doesn't work in 2D, where the normalized cross-product produces a scalar result.. or a unit vector along the z axis, depending on how you look at it.

Instead, I need to describe the directions using normalized vectors: p0/|p0|

Hence, my post. Eran Guendelman's math reference provides a derivation of D(x/|x|), however it disagrees with other sources I've come across:
http://graphics.stanford.edu/~erang/doc/reference.pdf

Does anyone have any experience with this?


Aside: I'm interested in extending the method presented in the paper to deal with rigid bodies; in 2D this seems like it should be relatively straightforward, the constraint functions will just be functions of the position and orientation of the bodies, and the gradient/jacobian of these functions should be very similar to those found in IK literature.

A point-to-point constraint would take the following form: C(p0,a0,p1,a1) = p0 + R_a0(r0) - p1 - R_a1(r1) (where p,a are the position and orientation of a body, and R_a(r) is the function which rotates the vector r through angle a).

Are the derivatives wrt p0 and a0:
Dp0(C) = I, Da0(C) = Perp(r0)? Or am I completely on the wrong track?


Finally, instead of explicitly calculating the partial derivatives, couldn't I simply use finite differences (i.e evaluating the constraint function directly with several inputs) to approximate the jacobian? I've seen some reduced-coordinate articulated-body papers where this approach is used, I'm wondering why it's not more commonly used as an alternative to explicitly figuring out the jacobian.


Please forgive this long post, and any incredibly stupid mistakes I'm making!

thanks,
raigan

First let's build up some identities we can use. Let r be a vector ( x(t), y(t), z(t) )

C = |r| = Sqrt( x*x + y*y + z*z )

dC/dx = 1 / 2 * Sqrt( ... ) * 2 * x
dC/dy = 1 / 2 * Sqrt( ... ) * 2 * y
dC/dz = 1 / 2 * Sqrt( ... ) * 2 * z

-> dC/dr = r/|r|

You can also look in the Baraff course for this one. All those identities can easily be derived looking at the components...

For the building the constraint yo now have to options. You build the time derivative and find the Jacobian by inspection:

dC/dt = dC/dx * dx/dt = J * v

dC/dt = [ v * |r| - r * ( r / |r| * v ) ] / |r|^2
dC/dt = [ v - r / |r| * ( r / |r| * v ) ] / [r]

We substutitue n = r / |r|

dC/dt = [ v - n * ( n * v ) ] / [r]

What you might have missed is how to deal with the triple scalar product which is NOT associative. I will use nT to indicate that I switch to a row vector.

n * ( n * v ) = n * ( nT * v )

This is now a matrix representation which is indeed associative so we can write:

n * ( nT * v ) = ( n * nT )* v

We plug this in our equation:

dC/dt = [ v - ( n * nT )* v ] / [r]

We now multiply the left velocity term by the identity matrix I

dC/dt = [ I * v - ( n * nT )* v ] / [r]

Now we can factor v and find the Jacobian by inspection:

dC/dt = 1 / |r| * [ I * v - ( n * nT ) ] * v

-> dC/dr = J = 1 / |r| * [ I * v - ( n * nT ) ]

You can resubstitute n = r/|r| and get arrive here:

-> dC/dr = J = 1 / |r| * [ I * v - r * rT / rT * r ]

Note that the xT * x is a "kosmetic" formulation for the dot product and could be skipped. It simply looks good.

You can reformulate this as you like to the other formulas in the Guendelman paper.

The other method is to directly build the gradient. We have a vector valued function V devided by a scalar function f. As mentioned in the e-mail classification is the key to success here. So the formula you need is

D( V/f ) = ( f * Trans( DV ) - V * Trans( Df ) ) / f^2

V = r -> DV = dr/dr = I (Identity matrix - good practise for you to proove this)
f = |r| -> Df = r/|r| (We derived this at the beginning of the post)

Now you simply plug this in the formula:

dC/dr = [ |r| * Trans( I ) - Trans( r/|r| ) * r ] / |r|^2

-> dC/dr = 1 / |r| * [ I * v - r * rT / rT * r ]

Let me know if this helps. It is quite late here. If you want I derive you the gradient of the normalized cross product tomorrow...

HTH,
-Dirk

I've just worked through this several times and I'm pretty sure I understand! At least, the answers were the same each time

I was tripping up on the transpose; for some reason whenever i see a vector transposed I assume that the author is just trying to make the dotproduct look fancy/nice; once I realized that xTx is dotproduct, xxT is outer product, things made more sense. The other mistake I was making was completely ignoring the transposes altogether (i.e using the quotient rule for scalar functions.. whoops).

I find directly building the gradient to be much easier, that might be because I have little experience with the time derivatives of vector functions -- for instance it took me a while to realize that the reason you went from (d|r|/dt) to (r / |r| * v) was the d|r|/dt = d|r|/dr * dr/dt.


I'd be interested to know your thoughts on simply using finite differences to approximate the Jacobian. It may be less accurate, but you're guaranteed to not have any bugs in your Jacobian


Also, I've come across Dean Macri's GDC2000 paper, which has a vector calculus intro, as well as formulas for bending constraints: http://www.gamasutra.com/features/gdcarchive/2000/macri.doc

Unfortunately, I've also found David Pritchard's "Implementing Baraff & Witkin's Cloth Simulation", where he claims Macri got several things wrong -- including using arccos to formulate the bending constraint: http://davidpritchard.org/freecloth/docs/report.pdf

The reason this is a problem is that the more recent "Position Based Dynamics" paper also uses arccos.. tomorrow I'm going to try to work through both versions and see which makes more sense.

thanks again,
raigan

p.s - my results for a 2D angle constraint (assuming it's valid to use arccos) that maintains the angle theta between the vectors (b-a) and (c-a):

C(a,b,c) = arccos( (b-a)/|b-a| . (c-a)/|c-a| ) - theta

-translate b,c to be relative to a (the relative origin), to get:

C(a,b,c) = arccos( b/|b| . c/|c| ) - theta

-let Qx = partial derivative of C wrt vector x
-let d = b/|b| . c/|c|

Qb is the column:
[bx*by*cx]
[bx*by*cy]

scaled by: 1 / ( sqrt(1 - d*d) * |b|^3 * |c|)

Qc is just Qb with c,b swapped, Qa = -Qb - Qc

There is a paper that deals with automatic differentiation. I believe it was from Grinspun. Maybe you look for a paper called "Thin Shells". It might be mentioned there. I think there was also example code for it...


HTH,
-Dirk

A lot of stuff related to automatic differentiation can be found on autodiff.org
(if I have the same understanding of automatic differentiation you have)

Noehrgel

I meant this here (s. bottom of page):

http://www.multires.caltech.edu/software/

The paper this is related to is also quite interessting...

So, has anyone implemented a bend constraint using the "Position Based Dynamics" paper?

The fact that they use Dot(normalA,normalB) to measure the bend angle seems wrong -- dotproduct's 180deg symmetry means that their constraint won't be able to distinguish between +90 and -90 for instance. Is this just something that never comes up with a small enough step?

Thanks for the automatic differentiation links; I actually came across autodiff.com yesterday, but didn't see the discrete shells paper.

Yes, I implemented the bending constraint, but don't use if because of performance reasons.

I remember that there were problems with the acos. Can you shortly fresh up my mind what might happen?


Cheers,
-Dirk

hi,
Pritchard's partical derivative formulas for the bend constraint are different than Macri's; in section 5.1 Pritchard talks about this, but I haven't worked through them yet.

My problem is that dotproduct isn't really a measurement of the angle between two vectors. The bending constraint proposed in posBasedDyn.pdf is: C = acos( Dot(n1,n2) ) - theta

but, consider: a=(1,0), b=(0,1), c=(0,-1)
a.b == a.c, but the angle between a and b is -1* the angle between a and c!

It makes more sense to me to describe the constraint as:
C = (atan2(n1) - atan2(n2)) - theta

Otherwise, there will be two configurations which satisfy the constraint, and the particles will be pushed towards the closest of the two.

I just finished writing a small test; using finite differences to approximate the Jacobian behaves just as well as using an explicit formula. It does seem to be slightly slower, and things might not behave as well with other constraint functions. I'll know soon, I'm going to try implementing bend constraints..

<new: so far, my angle constraint isn't behaving correctly.. looks like i need to work on this a bit more>

I had the problem with the dot product as well in rigid body dynamics. The constraint can indeed flip, since there are two answers. There are several formulations for the same constraint, but with different numerical properties. This is something I am heavily thinking about currently. E.g. a simple distance constraint

C(p1, p2) = |p1 - p2| - L0
C(p1, p2) = 0.5 * [ (p1 - p2)^2 - L0^2 ]

Which one is better? What are guidelines when coosing a constraint function? Gleicher tries to give some "intuitive" recommendations in his PhD, but I am not sure how correct this is. For the bending constraint we could also use:

C(n1, n2) = n1 * n2 - cos(alpha0)

or reformulate to use the sin() like Fedkin/Bridson. I don't find this easy to answer correctly.

When you use a dot product as position constraint I don't think in terms of angle. The constraint is driven by the projection of a onto b and you want that this projection should vanish...


Cheers,
-Dirk

I'm glad I'm not the only one wondering about this stuff!

I fixed my angle constraint: as usual, I forgot that when dealing with 2D angles you always have to wrap deltas so that they take the shorter arc.

I'm amazed by how well the finite difference approach works (Pritchard calls it "numeric derivatives"); without much tweaking (i.e finding the optimal delta to apply for each difference, or using backward/central differences instead of simple forward difference) it's only a factor of 2-4 times slower than the hand-derived version for distance constraints.

It's amazing because it's very easy to try arbitrary constraint functions and see how they behave.. all with basically no math! Tomorrow I'll try several angle formulations; my current angle constraint function is:
C = (atan2(n1) - atan2(n2)) - theta

Of course, this is probably not the easiest formulation to derive, but since I didn't have to do any derivation (yet), that's no problem.

I'm sorry for being so slow earlier, today I finally realized that something like: C = atan2(n), n = x/|x|
isn't simply a function, it's a function (C) of a function (atan2) of a function (n) of x.. which of course necessitates a more involved derivative. Did I mention I barely made it through first year calc?


Here's a fast demo, just a simple chain of sticks and 90deg angles (click near the top-left to set focus): http://www.harveycartel.org/raigan/pbdTest05.swf
turn on capslock to run it, or click with caps off to tick; hold z to grab the first particle w/ the mouse, spacebar to reset.

This is 5 iterations of constraint projection; one iteration works perfectly fine, except things get a bit "bouncy"/excitable. They never blow up though; some short tests suggest that choosing a better delta for the finite differences might help this a lot.

I'll be happy to share the code, but frankly it's exactly what you'd expect: each constraint is just a function that returns a scalar given a set of particle positions. That's all that's needed to build approximate gradients via differences, which are then used in the general constraint formula given by eq.8 in the "PBD" paper.


I've never gotten this far with angle constraints; previously I've had a simple impulse-based system with Baumgarte joint limits, but not only was that much less well-behaved, it was also a lot harder for me to understand.

Simply pushing the particles down the gradient of the constraint function not only makes perfect sense, it's braindead to implement with finite differences! This seems too good to be true.

I'll let you know if I come up with anything interesting, thanks again for all of the help, I really appreciate it.

About choosing between different constraint functions: a useful tool would be to graph the "contour lines" (or whatever that's called) of the functions. Do you know of any tools for that? I guess Maple would probably do.

Since we're just moving down the gradient, the more circular/less irregular the contours, the better the convergence should be.

Aside from that I have no idea. I've seen that 0.5*(curlen^2 -restlen^2) distance formula before, and wondered what the point of it was. Maybe just avoiding a sqrt?

You are right. Maybe one should look onto the isolines of the constraint function. Maple, Mathematica and the like should do this...

This first one is better because:
1. Other position constraints are likely to have units of length, so your position stabilization will be balanced.
2. The associated Lagrange multiplier will have units of force. This is important if you are making a distance motor or need to specify a breaking force.

I'd never thought about the units thing before, that's a good point.

I've been reading through some "geometric constraints" papers (it seems like they're mostly dealing with constraining segments and vertices in CAD programs), several have claimed that using error squared results in faster convergence. However, I haven't found any that actually explain why this might be or prove it in any way.

I guess squaring the error would make the constraint function "steeper", which might make a single constraint converge faster, but does that still hold when there are multiple competing constraints?