This post just needed to be done, to mark again something that was known for some time and that uses once again the concept of projection. I am calling "projective constraints" what Baraff and Witkin called conjugate gradient (CG) filters in their break of the century seminal paper on cloth. I took the name from the SOFA framework documentation where it appears in a slightly different (but simpler) form, sign that this concept was previously known (which is almost always the case). In fact, a very similar form of CG appears in the book by Nocedal and Wright on numerical optimization.
I will try to be brief and acknowledge my main source as being the excellent tutorial on deformable body simulation from Pixar written by David Eberly and Theodore Kim. They are the ones to extract an important idea from a 2015 paper on multigrid cloth simulation by Tamstorf et al. And that idea is that you don't necessarily need CG to filter (or project) constraints, but you can do it with any linear system solver of choice as long you transform your system into the right space. This was an idea that was explored before by Ascher, Boxermann and Ye, but they all focused on the idea of a reduced system (of a smaller and even varying size). The pre-filtered preconditioned CG (PPCG) method, as it was later called, got one thing right: the correct form of the transformed system that was well-conditioned and works great in practice. In retrospective this looks moderately easy to derive, but somehow no one managed to derive it correctly before.
Hopefully, I will not mess up equations too much and add to the confusion, but here are the main lines of reasoning. You start with a linear system $\mathbf{Ax}=\mathbf{b}$. Then, you have a projection matrix $\mathbf{S}$ that takes you from the full configuration space to a smaller one where motion is restricted by some constraints. The letter $\mathbf{S}$ is chosen historically and also to differentiate from its orthogonal complement $\mathbf{P} = \mathbf{I} - \mathbf{S}$ which takes you onto the constraint space as was described on this blog before. Keep in mind also that orthogonal projection matrices are symmetric.
By the rules in which matrices and vectors transform under a change of basis, we can move our linear system into the "free" or "unconstrained" space: $\mathbf{SASx}=\mathbf{Sb}$. Or at least it would seem so. The more correct form is the one from the PPCG paper:
$(\mathbf{SAS} + \mathbf{I} - \mathbf{S}) \mathbf{y}=\mathbf{S}(\mathbf{b} - \mathbf{Az})$.
We assumed $\mathbf{z}$ is a prescribed solution for the constrained dimensions of the system. This means the final solution is a mixture between these and the free degrees of freedom: $\mathbf{x} = \mathbf{y} + \mathbf{z}$. Hence, $\mathbf{y} = \mathbf{x} - \mathbf{z}$. You need to be a bit extra careful to project $\mathbf{z}$ by $\mathbf{I} - \mathbf{S}$ whenever you are adding it to $\mathbf{x}$ so that you don't add contributions from the wrong dimensions, but if we consider it from the start as being in the right space we shouldn't worry about it. Finally, we are solving for $\mathbf{y}$ which is in the same space as $\mathbf{x}$ but it simply doesn't move in the constrained directions. We took care of it not happening by modifying the system matrix. Moreover, the system is perfectly solvable as the $\mathbf{I} - \mathbf{S}$ term guarantees no rank deficiency (and one can use a direct solver too).
I admit I don't fully grasp every little reasoning and subtlety in the derivation. And it's always easy to get yourself confused between one projection matrix and its complement. I swapped them by mistake many times and it takes a bit of practice and study of the original papers to get it right, but it's doable.
The only thing that I want to draw attention to is this remark from the PPCG paper: "However, the above method easily generalizes to any kind of equality constraint, including the face-vertex and edge-edge constraints used in other collision response algorithms." Indeed, this is a step beyond the original CG filters, but it begs the question how to do it, as it's not directly obvious. And this is where the constraints projection matrix comes in handy and the relation $\mathbf{S} = \mathbf{I} - \mathbf{P}$. Now, in our previous discussions and the constrained dynamics literature, this matrix involves the inversion of the Delassus operator (or Schur complement matrix). But, here, we can consider the metric of the configuration space as being the identity, as all the authors above have considered, and we get a simpler form: $\mathbf{P} = \mathbf{J}^T(\mathbf{JJ}^T)^{-1}\mathbf{J}$, where $\mathbf{J}$ is the Jacobian of the constraint functions.
And immediately we notice that we may run into problems if $\mathbf{JJ}^T$ is rank deficient (i.e., the constraints are redundant). But this can be solved by performing a QR or SVD decomposition of $\mathbf{J}$ in order to compute the required inverse. Moreover, in my brief experiments I found that taking only some of the first singular values is enough to make it work. There is still much work to be done here, but it shows in principle that this method is feasible in practice. And if the (approximate) factorization of $\mathbf{J}$ is done fast enough it can outcompete the Gauss-Seidel solver for the Schur complement system usually employed for constraints. (I think there is a paper by Jernej Barbic that did it some time ago too.) But I'm not sure if in the end it's really worth it (in terms of speed or solving inequality constraints).
All I care about is the beauty of the connections. And I'm just glad that I know now how to better understand CG filters and how to (potentially) extend them to arbitrary constraints. (The only paper that I know of that goes on a similar path is "An Implicit Frictional Contact Solver for Adaptive Cloth Simulation" by Li et al.) Who knows? Maybe in the future we will find out that solving constraints this way is more robust, and cures all the problems of dual relaxation solvers like slow convergence to zero constraint error or large mass ratios.
But should we use the proper mass matrix in the projection matrix? And how do we find the prescribed $\mathbf{z}$ values? Isn't that a full projection step that we do to solve constraints? And how are the free degrees of freedom interacting with the constraints? (It seems we are freezing the constraint solutions when in fact there should be an equilibrium between them and the rest of the system.) I think these are question that have to be answered in the wider scope of coupling elasto-dynamics and constraints (which may still offer some surprise new solvers over the next years). And projective constraints may play an important role in solving this problem, or at least an alternative.
