The merits of projecting twice

I started a couple of years ago writing a post called "The trouble with PBD contacts", but never finished. This was based on my practical experience of modeling contact. I hope I will finish that post one day, once I learn all the factors causing all that trouble. But long story short, PBD contacts are pretty crap and usually quite violent. As Stewart puts it: "the resulting discretization behaves as if it had a "random" coefficient of restitution when impacts occur". And the same applies to stabilized velocity contact constraints. In practice, this manifests as objects popping out sometimes at high speed. I call them "violent contacts". Bottom line: the velocity at impact should always be zero and this should be our constraint (unless we want to model non-zero restitution, but we can add that afterwards). If we do that, then all our contacts will be just fine and no explosions will appear. But what about the penetration error? I'm sure you will ask that as nobody likes objects on screen clipping each other. Stay tuned.

I will focus here mostly on what I call position-based methods, which is different from position-based dynamics (PBD), but it includes it. There are generally two types of such methods: explicit and implicit, pertaining of course to the type of integrator used. But, more importantly, is that it also makes the distinction between explicit and implicit constraint directions. Meaning, when is the direction of the constraint forces computed: at the beginning or the end of the time-step? In the first category I put SHAKE and RATTLE which are the darlings of molecular dynamics. They are explicit symplectic integrators of constrained systems. Hence, they can easily blow up or end up without a solution for large time steps, but people live with it as long as energy is conserved. On the other hand, implicit methods tend to be much more stable, as by construction they dissipate energy (and people live with it). This is why they are used so much in graphics and gaming. In this category I put the so-called "projection" methods, where you integrate up to a candidate position (without constraints) and then project the positions back onto the constraint manifold. They are also known as step-and-project (SAP). Here's where I put PBD and "fast projection" too, as they are approximate solvers for SAP. In fact, I will often say PBD instead of SAP, unless I make it clear that I'm referring to the nonlinear Gauss-Seidel projection popularized by Nvidia.

I will not go into the nitty gritty details of equations, but there is one important point I want to make. Neither SHAKE nor PBD preserve the hidden constraint that velocities along constraint directions should be zero (the normal component to the constraint manifold or the first derivative of the constraint equation). And I will prove it with a plot of the velocity constraint drift for the case of a simple pendulum.

RATTLE is the only scheme that calculates the correct velocity. And it does it at the price of an extra projection step at the end for the velocity. Given that the position constraint (explicit) projection at the beginning of RATTLE is nonlinear, and the velocity projection is always linear, it seems like a small price to pay. It's only one linear system to solve, compared to several if using Newton for the positions. So, here's what I did, I took RATTLE and I replaced the first part with PBD. And the results were great in terms of velocity. I called it PROJ2 for lack of a better name and because I project twice: first the positions and then the velocities.

But why should we bother? SHAKE and PBD work perfectly well the way they are. Indeed, if you are dealing with equality constraints. But if you are dealing with contacts (and friction) velocity starts to matter more. That high velocity error in the plot will become the impact reflected velocity. This is not an issue with equality constraints, as no matter what velocity they have along the constraint direction, the bodies will stay pinned to the constrained manifold (up to a solver tolerance). So, this is why I modified RATTLE like I did. It's basically the modification you need to do to PBD in order to make it explosion free for contacts.

I have not not tried PROJ2 on practical contacts yet, but I did something similar in the past. I projected the velocities before running PBD and everything was great. If you think about it, it's the same thing, if you move the last step of RATTLE at the beginning of the next frame. So, I know it works! It made contacts much smoother. The only reluctance I had, like everybody else, was that it was more expensive. And now I know the answer: it doesn't matter! It's worth the price if your contacts are fully inelastic and not exploding. And you also get rid of that stabilization term... 

Think about it: first project velocities and then positions. Does it sound familiar? The latter step is sometimes called "post-projection": where you correct the positions after you solve the problem at velocity level in order to remove any penetration. This technique is usually attributed to a paper by Cline and Pai, but ultimately it's just projection. Usually the problem is linearized, but when dealing with contact, the normals are usually considered constant throughout the solve, so it is linear. Of course, it would be much better if you used PBD instead, but the border between post-projection and PBD is quite gray to be honest anyway.

One more thing, what's wrong with stabilization? Nothing, if you are writing a rigid body simulator where you don't care so much about small penetrations. And then you can use a small stabilization factor and clamp the velocity to avoid explosions. But if you care about having zero penetrations and correct normal velocities, you should use the post-projection step. In fact, this has been a debate for years in the community, and for that reason Box2D uses post-projection, while most other engines uses stabilization or a mixture of the two. Now, I tend to agree with Box2D. By the way, if you care about that little analogy between implicit springs and the stabilization factor and the "constraint force mixing" factor, I would say forget it - I think it's pretty much wrong anyway (or at least not rigorous). 

In the past, I feared that post-projection dissipates artificially too much energy. But this phase plot for the pendulum made me feel better.

In fact, it preserves energy much better than PBD. So, it has two merits: it reduces dissipation for PBD and solves penetration problem for velocity level methods. From now on I think this is what I'm going to use as a contact solver (and constraints in general), no matter how expensive it is.

The many faces of projection

What I want to talk about in fact is the Jacobian matrix of the constraint function. All these things that I want to mention are probably clear to people with strong backgrounds in linear algebra and differential geometry. But for me they were revelations of sort many of them.

All right, so maybe the first thing to start with is the fact the Jacobian matrix $\mathbf{J}$ of a function $g$ is the analog of the gradient operator or the normal to a surface. And this is why for many people $\mathbf{J}$ is synonymous to the constraint direction(s). Because the normal to the constraint manifold is transversal to the space of motion, hence along the constrained directions that can't move. Think of a pendulum: the constraint manifold is a circle and the normal to it is along the radius, which is exactly the pendulum's (constrained) direction.

But some math books also call $\mathbf{J}$ the differential (in coordinates). And this takes me to a confusion I always made in the beginning: for a single constraint, should $\mathbf{J}$ be a row or column vector? And it seemed other people weren't sure about it either in their notation. But I stuck with it being a row vector, as that was the majority view. And now it makes sense, because it is a covector (mainly because it has the components of a differential which is a 1-form, same thing). Hence, $\mathbf{J} \equiv dg$ and $\mathbf{J}^T \equiv \nabla g$.

I still don't fully grasp the concepts and that's why I don't want to talk formally about these things. But I do have to stress the fact that the configuration space of our mechanical system (where all the generalized coordinates reside) is not really $\mathbb{R}^n$, but one that has a metric (a curved space) given by the mass matrix $\mathbf{M}$ (i.e., the kinetic energy). It basically tells you that not all bodies are equal but weighed by their mass information. If this weren't the case, then vectors and covectors in $\mathbb{R}^n$ would be the same thing. But they're not! If for velocity (tangent) vectors $\mathbf{v}$ the norm is given by the well known kinetic energy ($\tfrac{1}{2}\mathbf{v}^T \mathbf{Mv}$), then for covectors the matrix is now the inverse of  $\mathbf{M}$. Think about momenta $\mathbf{p}$, the kinetic energy is the same, but now we have to divide by the mass rather than multiply ($\tfrac{1}{2}\mathbf{p}^T \mathbf{M}^{-1}\mathbf{p}$). And you have to multiply velocities by mass to get momenta $\mathbf{p} = \mathbf{Mv}$ and the other way around. Some other covectors are forces (as they work directly on the momenta) and surface normals; and they inhabit a place called the cotangent space.

Now let's fall back from geometry to linear algebra. Locally, $\mathbf{J}$ is a linearization of $g$ and maps the tangent space of the configuration space to the tangent space of the "constraint" space $\mathbb{R}^m$. So, it's a linear transform, thus a rectangular matrix. And it follows all the rules of linear algebra. Ideally, it's full rank and fat (rather than thin, i.e., $m <n$) in order for the system to be solvable. The image of  $\mathbf{J}$ is the whole of $\mathbb{R}^m$, which can be viewed as just a collection of numbers, just as for one constraint you get only one number. The meaning of these numbers comes from the analogy with the normal vector, in fact any vector for that matter. They are the magnitudes or components of the projection of a vector $\mathbf{v}$ on the "normal" $\mathbf{J}$, so $\mathbf{Jv}$ is just that. For one constraint it is a number, in general it is a $m$-dimensional vector because the "normal" $\mathbf{J}$ spans an $m$-dimensional space. You can think of taking each row of $\mathbf{J}$ and dotting it with $\mathbf{v}$ to get the projection length. And in fact, from linear algebra we know that the row space of $\mathbf{J}$ is formed by vectors $\mathbf{J}^T\mathbf{y}$, so we know that the rows of $\mathbf{J}$ span the "normal" or "constraint" space of dimension $m$. This space is orthogonal to the tangent space of the constraint manifold, which is precisely the null space of $\mathbf{J}$, i.e., vectors satisfying $\mathbf{Jv}=0$.

What does this tell us about $\mathbf{J}$? That it takes any vector $\mathbf{v}$ and converts it into its projection components. What about the reverse direction? It will only reconstruct the "normal" part of those vectors. So, you lose some information, but this is fine because you obtain the projected vector $\mathbf{J}^T\mathbf{Jv}$ (almost, stay tuned). If you think about it, this is what you do when you project on a unit vector $\mathbf{n}$: $\mathbf{nn}^T\mathbf{v}$. So you may say that $\mathbf{J}^T\mathbf{Jv}$ is a projection operator. Note that if $\mathbf{v}$ is already in the row space of $\mathbf{J}$ no information is lost and projection leaves it unchanged. Some people call $\mathbf{Jv}$ a "projection" and $\mathbf{J}^T\mathbf{v}$ an "un-projection", but this is not really true, as they do not satisfy the definition of a projection operator; but they do part of the job, so we may call these matrices something like "projectors" (when in fact they define the space we are projecting upon).

But this not the whole truth. In 3D you also have to divide by the norm squared $\mathbf{n}^T\mathbf{n}$ if we're dealing with a non-unit vector. And you have to do the same in $n$-dimensional space. Linear algebra and least square problems give you the formula $\mathbf{A}^{-1}\mathbf{Jv}$ for computing the actual projection lengths, where $\mathbf{A}=\mathbf{JJ}^T$. Here $\mathbf{A}$ plays the role of the norm squared of  $\mathbf{J}$ that you divide by. Finally, the projection that takes these lengths back to configuration space is $\mathbf{P}=\mathbf{JA}^{-1}\mathbf{J}^T$ and this is a genuine projection operator. Of course, if you're lucky and $\mathbf{A}$ is the identity then the previous paragraph still works, but that's rarely the case.

The projection $\mathbf{P}$ or rather its complement $\mathbf{I}-\mathbf{P}$ is usually called the "canonical" projection on the constraint manifold. And you may recognize the formula from the analytical solution of linearized velocity-level constrained dynamics. In practice you usually do not invert $\mathbf{A}$, but rather solve a linear system where the unknowns are the projection lengths. But $\mathbf{A}$ is not really the same, is it? 

So, what is the meaning behind $\mathbf{A}$ anyway? Some people call it the reduced or effective mass, others call it the Delassus operator, others compliance and so on. I think I hinted already at it: $\mathbf{A}$ determines the metric of the "constraint" space. In fact, it's the inverse $\mathbf{A}^{-1}$ that is the true effective mass, the one seen by the degrees of freedom along the constraint directions, and helps define the kinetic energy there. (I picture it as the transformation of the mass matrix $\mathbf{M}$ into "constraint" space). And $\mathbf{A}$ is actually the metric of the cotangent "constraint" space. This makes sense as $\mathbf{J}$ spans the cotangent space of the constraint manifold (in configuration space) and there the metric is $\mathbf{M}^{-1}$, so our "norm squared" becomes $\mathbf{A}=\mathbf{JM}^{-1}\mathbf{J}^T$. This is exactly the matrix that appears in constrained dynamics.

I want to end by noting the fun fact that by deriving a projection operator (pretty similar to what you would do in 3D), we actually solved for the Lagrange multipliers in constrained dynamics. In fact, projection is an optimization problem: find the shortest distance from a point to a surface. And the closest point happens to be on the perpendicular to the surface, one like our "normal" $\mathbf{J}$. Pretty much what we talked about. So, deriving the projector $\mathbf{P}$ amounts to solving the projection (or least squares) problem. Or it does at least for velocity vectors. The general case for projecting a point in configuration space on the constraint manifold is a bit trickier as it's a nonlinear problem, but it can be broken down into such linear steps. But that's another story...

To summarize, in order to solve the (linearized) constrained dynamics problem one simply needs to project velocities into the tangent space which is the null space of $\mathbf{J}$. This is done by projecting them onto the space spanned by the columns of $\mathbf{J}^T$ and then taking the orthogonal complement of that projection.

Mixed FEM and constrained dynamics

 I recently got a paper accepted, which in itself was an achievement and culmination of a long period of work, hardship and even failures. But the idea I want to bring forward in this post is that the paper only scratches the surface of the implications it can have and could be developed further. Such is the way of academic publishing: not always the results I or you deem important are considered results worth publishing by the academic community. And I won't go that far into saying that my results were revolutionary. I'm merely saying that they are interesting and they may settle a decades-long debate in the engineering, graphics and simulation circles.

But let me tell you more before I get to that. The mixed formulation of the finite element method (FEM) is a technique that has been around at least from the 70's. Its main purpose was to solve incompressible problems in fluid flow and elasticity. This was due to a problem called locking that is still affecting our simulations. As a simple example, locking happens when you simulate a cantilever bar bending under some load; and although you know it should deflect by quite a bit, the numerical results tell you the deflection is zero or very close to zero. But let's not worry about locking, as even if avoiding it may seem like the main result, I argue that it isn't.

The main achievement of mixed FEM is its beauty and unifying power. So what is it? I will make an analogy with springs. The opposite of mixed FEM is irreducible FEM, or I call it standard FEM, because it's used all over the place. Spring forces (like stresses in a finite element) are given by a constitutive law like Hooke's law: the elongation times a spring constant. These forces are then plugged in some equilibrium conditions or Newton's second law to obtain a simulation of the system. The main takeaway here is that you always need a constitutive law of some sorts besides the laws of statics and dynamics to completely determine the behavior of a mechanical system.

What if we don't substitute the constitutive law into the equations of mechanics as we usually do? We get the mixed formulation of springs, and FEM and everything else. We just get (roughly) double the number of equations that we need to solve, but the problem is essentially the same. In fact, it can be argued this formulation is "more" fundamental, and the irreducible form is one that has been reduced from it (and cannot be reduced any further). But why wouldn't you substitute and reduce the equations? Apparently, no reason at all.

This takes us a long lasting debate between choosing reduced coordinates over maximal coordinates. That is, if a system is constrained in any way (like a pendulum attached to a point) then it makes complete sense to consider that some degrees of freedom were removed by the constraints and use only the remaining ones. Or not... The view of constrained dynamics is that you should still use all the initial degrees of freedom of the pendulum and add extra equations and variables for each constraint. Understandingly, a lot of people don't like this, especially in their computer implementations. But it just so happens that this is the best way to simulate constrained dynamics when you don't know from the start what the constraints will be (i.e. model an arbitrary system).

On the other hands, constraints are only ideal and do not really exist in nature. They are just really strong forces in the end. So, there is a middle way: use all the degrees of freedom and model constraints using very very strong forces (like the ones in springs). This keeps the number of variables and equations in the middle. It is also the standard or irreducible form.

But something quite bad happens when the spring constants grow very large. In theory they should tend to infinity so that the constraints are exactly fulfilled. And this is not possible in practice. Numerics become ill-conditioned and, on top of it, extra problems appear like locking. And these won't go away by just putting more power into the number crunching. You need to roll back to the reducible or mixed form of the equations. Or the one that had both constraint equations and equations of motion. You need to make the problem larger in order to make it solvable. And here comes the big revelation: the constitutive laws in the mixed formulation are just soft constraints in constrained dynamics.

And of course people knew about it. At least many do or have a strong intuition about it. But the rest act as if they don't. There are actually two camps in solving contacts and other topics which divide the world into penalty vs. constraint formulations. When they are in fact two sides of the same coin. In their defense, the constraint camp does assume that their constraints are infinitely hard which is not possible when using finite spring constants or material parameters (and despite being ill-conditioned their problems are still solvable, even if by a bit of cheating). But the mixed formulation is just a place where these two different views of the world meet. In technical terms it is called a saddle point problem and keeps popping up whenever you do simulation (see Strang's textbooks for example). And optimization people have known it for a long a time: they call it the KKT optimality conditions. For them the degrees of freedom of the system are the primal variables and the tensions inside the springs are dual variables or Lagrange multipliers.

My personal breakthrough was to show that any simulation problem, no matter how complex (including PDEs), can be formulated as such a mixed formulation involving both the equations of motion and the soft constraints (aka constitutive laws). If any of the constraints needs to be hard, we just put a zero value next to it and everything still works. So, we can build a single unified solver for all mechanical phenomena, including say FEM and contacts, by using the same framework. Of course, that doesn't mean it will be optimal; but sometimes elegance and generalization is better than some optimizations here and there.

And, ok, you would argue, this was still known. Yes, yet nobody used that knowledge. For example, you can build a whole multibody simulator by just treating every constitutive law used as a soft constraint and then use a single constraint solver like XPBD (or a more accurate one, I don't care). Or you can use very stiff springs (the penalty approach) as long as you find a nice stable way to solve the problem. And in the end it may not even matter as optimization theory is already using techniques that rely on both facets of the problem, like the augmented Lagrangian method. You add a bit of penalty, a bit of constraints and you get to solve your problem efficiently.

So what's left to do? I would say a lot. We need to build these unified multibody simulators using the mixed framework and see how they perform. What is the best approach for solving them? Sometimes, you have to pay a bigger price or be a bit slower than the state of the art to reach this degree of generality. Also, pay attention to traps like locking, which become visible in the mixed formulation, but require deep changes at the discretization level in order to be fixed. Overall, I would say it gives us a much higher vantage point on things. What else? Maybe people using implicit integration of stiff springs will stop claiming that they are doing a very different thing from people using constraints. Or claiming you can not simulate some fancy nonlinear dynamics using constraints (which you can by using the three field mixed formulation). Or constraint people realizing that regularization is good from time to time, that relaxation solvers like Gauss-Seidel are in fact already doing it and that not every problem is a nail. And that at some point you can write a solver that takes the best of both worlds. For example, one that tackles high stiffness through constraint solvers and large mass ratios through primal solvers (at the same time, in a single pass, inside one monolithic hybrid solver).