Implicit velocity projection

During the time I spent looking closer at PBD and what it really is I learned one clear thing: PBD is nothing else than a form of implicit integration (Backward Euler) of the constraint forces. Goldenthal calls it Implicit Constraint Directions (ICD). His optimized version he called Fast Projection (FP). FP solves a linear system at every nonlinear iteration. PBD runs just one Gauss-Seidel (or Jacobi) step at every iteration, resulting in what is known as Nonlinear Gauss-Seidel (NGS). 

The other thing I learned is that PBD is very stable just because it is fully implicit. And this is why you can do cloth with PBD, but not with velocity based methods, which I'll denote hereon as velocity time stepping (VTS). This is because articulated systems can get very unstable with VTS. The short explanation is that VTS is using the constraint direction at the beginning of the frame (i.e. explicit directions) and not the one at the end (implicit directions). Therefore a part of the system is integrated explicitly (the part transversal to constraints) and it is well known that explicit integrators can get easily unstable and even blow up. The good part is that VTS is very similar to SHAKE (pretty much a linearized version of it) which is symplectic. This means energy gets conserved better in the transverse direction - I will detail these things in another post.

Unfortunately conserving energy is not always what we want. Stability comes a lot of time for free from implicit integrators because they dissipate energy artificially. This may be bad for science projects, but not for games, graphics and real-time simulation. Therefore I've been looking for a way to make VTS stable without going really the PBD way. The differences between the two comes a lot from what we are trying to solve for: VTS solves for velocity constraints, PBD for position constraints. So how can we keep the VTS way without altering it too much, and also get the stability of PBD? The answer: integrate using implicit constraint directions but keep the velocity constraints.

Let's throw in some equations for this case:

\begin{align}
Mv^{l+1} &= Mv^l +hf_{ext} + hJ(x^{l+1})^T \lambda \\
x^{l+1} &= x^l + hv^{l+1} \\
0 &= \tfrac{\beta}{h} c(x^l) + J(x^l)v^{l+1}
\end{align}

We can express this as a system in the unknowns $v$ and $\lambda$:
\begin{align}
Mv &= Mv^l +hf_{ext} + hJ(x^l + hv)^T \lambda \\
0 &= \tfrac{\beta}{h} c(x^l) + J(x^l)v = G(v)
\end{align}

While the second equation is linear, the first one is not and we need to apply Newton's method to it. The Jacobian for the whole system then becomes:

\begin{equation}
\begin{bmatrix}
M-h^2 H(x)\lambda & -h J(x)^T \\
J(x^l) & 0
\end{bmatrix},
\end{equation}
where $x = x^l + hv$.We can proceed just like in the case of FP and ignore the geometric stiffness term $K=H\lambda$, where $H=\nabla J = \nabla^2 c$, and also set the residue of the first of the first equation to be always 0. We then get the following linear system:
\begin{equation}
\begin{bmatrix}
M & -h J(x_0)^T \\
J(x^l) & 0
\end{bmatrix}
\begin{pmatrix}
\Delta v_1 \\ \Delta \lambda_1
\end{pmatrix} =
\begin{pmatrix}
0 \\ -G(v_0)
\end{pmatrix}
\end{equation}

We can now take directly a Schur complement and obtain the usual constraint solver equation:

\begin{equation}
hJM^{-1}J^T\Delta\lambda_1 + G(v_0) = 0,
\end{equation}

where we can take by approximation $J$ to be the same on both side, i.e. $J(x_0)$. Then the change  in velocity is $\Delta v_1 = hM^{-1}J^T\Delta\lambda_1$ and the update step is given by:
\begin{align}
v_1 &= v_0 + \Delta v_1, \\
x_1 &= x^l + h v_1.
\end{align}
If you iterate this scheme you get implicit VTS and it's way more stable than VTS. I only tested this on a 2D thread/cable made of distance constraints and I haven't got to make measurements, but it looks just like PBD. But the algorithm looks a lot like VTS - there are two main differences though:

1. at every iteration, after you compute the new velocities, you compute new positions too, and
2. at every iteration you recompute the Jacobians of the constraints based on the new positions.

Otherwise, everything stays the same. You could use an exact solver for the Lagrange multiplier or better just use one Gauss-Seidel step, which results basically in NGS. The only trick is that you don't set $\beta$ too close to 1 or it explodes (why, I don't know yet).

Here's how it looks like for some extreme initial velocity conditions:

The same scenario, using VTS.

And for completeness, the PBD case.