As an exercise I wrote a solver for a non-peneration constraint of three 1D spheres (intervals). Basically the solution is found by projecting to the nearest surfaces or edges on the configuration manifold. Adding mass to the picture was tricky, but it's possible by intersecting the desired solution lines with the 'mass plane' passing through the initial configuration point. This is a plane with a normal in the direction given by creating a vector with the mass of each sphere for each component. When all the spheres are equal mass, the normal direction is (1.0, 1.0, 1.0) as in the original mass-less case. I think the approach may extend to more than three bodies in 1D but I haven't tried it. Anyway there are $n^2-n$ possible solutions of an n-body system where all bodies are affected (i.e. the solution results in all 'n' bodies in contact, and there are $n^2-n$ possible orders), that might be a not-so-great asymptotic compared to sequential offsets.

This was possible to visualize since configuration space for three 1D positions is a 3D space (attached example). Visualizing the configuration space of two or more spheres in 2D or 3D is more difficult since it's either 4D or 6D minimum. One difference with position constraints in 1D is the configuration space is disconnected since spheres must cross through each other to exchange ordering (there's no other dimension to go "around" an obstruction in 1D, you can see this as the six wedge-shaped regions outside the manifold in the attached image). Does anyone know if a closed form solution (using just a finite number of projections, intersections, etc.) is even possible in the 2D (or 3D), for three or more spheres?

Well I think the answer might be the general case can be formulated as some kind of Linear Programming problem solvable by simplex method (I think this is what's called "Dantzig's algorithm"). Coutinho 2013 discusses using Dantzig's algorithm for computing frictionless contact forces, which is to say its solving a 2nd or 3rd-order problem when here I'm just considering a 1st-order constraint, so maybe it doesn't even need to be that complex. Under the assumptions of the problem (positive definite matrix, existence of a solution), it's claimed that this algorithm will terminate in a "finite number of steps", so I take that to mean it's not an "iterative method" that approximates, but an exact solution. It doesn't seem too bad but there is some bookkeeping going on with tracking what contacts are "active" or not and updating their state accordingly and the algorithm proceeds to add more constraints one by one. Then there's the problem of adding friction (Baraff's algorithm) which makes things even more ugly so I guess that's why iterative impulse-based methods are preferred to the force-acceleration model.

Edit: One final thing I have been searching for an answer on relating to iterative methods, what is "Nonlinear" referring to in "Nonlinear Gauss-Seidel" (NGS)? Is it just "normal" Gauss-Seidel applied to a nonlinear system of equations? If so when does the system of equations become non-linear, is it due to introducing rotation specifically, or does it arise when the constraints are linearly independent? (Linearly independent constraints arises in the 2D or 3D case, but in the 1D case the constraints are collinear if I am thinking about this correctly.)