A Unified Analysis of Penalty-Based Collision Energies

We analyze a wide class of penalty energies used for contact response through the lens of a reduced frame. Applying our analysis to both spring-based and barrier-based energies, we show that we can obtain closed-form, analytic eigensystems that can be used to guarantee positive semidefiniteness in implicit solvers. Our approach is both faster than direct numerical methods, and more robust than approximate methods such as Gauss-Newton. Over the course of our analysis, we investigate physical interpretations for two separate notions of length. Finally, we showcase the stability of our analysis on challenging strand, cloth, and volume scenarios with large timesteps on the order of 1/40 s.


INTRODUCTION
Robust collision response is an integral part of any deformable body simulation, such as those used in feature film production, games, and biology. Perhaps the simplest approach is a penalty-based response, which applies forces to primitives that have violated some distance threshold.
In order to obtain stable motion under large timesteps, implicit solvers require the gradients of these forces [Baraff and Witkin 1998]. However, these gradients can produce indefinite systems, so they are commonly filtered to be positive-semi-definite (PSD) by removing negative eigenvalues at each quadrature point using a numerical factorization [Teran et al. 2005], or by applying a Gauss-Newton approximation [Choi and Ko 2002;Kim and Eberle 2022;Zehnder et al. 2021]. A better approach is to use a fast, analytic eigendecomposition, but while these have been constructed for a wide variety of deformation energies [Kim 2020;Lin et al. 2022;Smith et al. 2019], penalty-based collisions have not undergone the same treatment. This is unfortunate, because we have found that they can provide significant gains in speed and robustness.
We perform this analysis by recasting collisions in terms of a novel reduced frame, C. From this perspective, a unified treatment emerges that encapsulates both vertex-face and edge-edge collisions. From there, we can perform a generic eigenanalysis that yields closed-form expressions for the eigendecomposition of penalty-based force gradients. Over the course of this analysis, we present a notion of signed and unsigned length that is applicable to a variety of collision energies. We then recast several commonly-used energies in terms of these length measures.
Our final analysis is performed on general penalty energies involving unsigned and signed length. We show a variety of scenarios where applying our approach is significantly faster than computing a numerical eigendecomposition, and more robust than using approximate force gradients.

Penalty-Based Collision Energies
2.1.1 Spring Energies. Penalty-based collision forces are a well-established method for handling contact in deformable body simulation [Baraff and Witkin 1998;Jimenez and Luciani 1993]. Object inter-penetrations are suppressed by inserting a spring between vertex-face or edge-edge pairs that cross a distance threshold. They have similarly been used for cloth at least since Baraff and Witkin [1998], though additional challenges arise when inside/outside queries become unavailable. For example, Teran et al. [2005] use quadratic energies to generate filtered force gradients for self-collisions that occur during muscle and flesh simulation.
Strand collisions can be viewed as a collection of edge-edge resolutions [Spillmann and Teschner 2007], and have additionally been used to simulate cohesion [Ward et al. 2004] in addition to contact repulsion. Many works [Chang et al. 2002;Choe et al. 2005;Selle et al. 2008] use penalty forces during hair simulation, and stiff springs to adhere joints in woven cloth [Michels et al. 2015].
Specific simulation goals have given rise to a variety of spring penalty formulations. Gast et al. [2015] introduce a linear falloff in the spring constant, yielding a cubic energy, and McAdams et al.
[2011] use a control variable to interpolate between isotropic and anisotropic springs. Spring-like forces are also used to simulate clumping in wet hair, where Lin et al. [2011] use "wetness" to determine the activation distance for adhesive springs at strand joints, while Fei et al. [2017] use cross-sectional surface areas to govern the strength of inter-strand cohesion.
2.1.2 Barrier Energies. Barrier-type energies contain terms that grow to infinity under special conditions, effectively prohibiting specific configurations while maintaining a penalty-like approach [Schüller et al. 2013]. Recently, Incremental Potential Contact (IPC) has made extensive use such energies for collisions, and utilize similar per-element PSD projections [Li et al. 2020]. This underlying framework is extensible to codimensional settings featuring volumes and cloth [Li et al. 2021], rigid bodies , and generalized solvers [Chen et al. 2022]. Our eigenanalysis extends equally well to the vertex-face and edge-edge energies used in these IPC-style solvers.

Analytic Eigenanalysis
Several techniques have been developed for analytically computing the eigensystems of various force gradients. However, these have all been for specific cases of hyperelastic materials, not the collision energies that are the focus of this paper. Smith et al. [2019] performed an eigenanalysis of hyperelastic deformation energies through the introduction of novel stretch-based invariants, which was later extended to membrane stretching energies [Kim 2020;Panetta 2020]. Lin et al. [2022] used a polynomial-based method on the Cauchy-Green invariants to perform an alternate analysis of ARAP-style [Sorkine and Alexa 2007] materials. All of these analyses offer broad techniques for constructing fast, analytic eigendecompositions, but to our knowledge, ours is the first to attempt a similar analysis on collision energies.

Complementary Approaches
Complementary collision response strategies have been developed to address more extreme collision scenarios, and deployed in addition to penalty-based energies. Such treatments include linear complementarity formulations [Bridson et al. 2003;Gascón et al. 2010], with careful considerations toward preventing artifacts that arise from direct vertex corrections. Eulerian grids have also been investigated for collision resolution, leveraging incompressibility conditions [McAdams et al. 2009;Teng et al. 2016] or velocity constraints [Levin et al. 2011] to maintain a penetration-free state. The success of these methods are then linked to grid resolution. Müller et al. [2015] construct volumetric "air" elements between objects that strictly prohibit negative volume, but element construction for complex scenes introduces additional difficulties.
Perhaps the closest approach to our reduced frame is in Jiang et al. [2017], where the R matrix from the QR factorization is used to construct a deformation-specific frame. However, their factorization is used to construct a new elastic potential, not analyze existing ones. We investigate directions, while they generate forces.

A REDUCED COLLISION FRAME
In the following, we will denote matrices capital bold (X), vectors with lower-case bold (x), and scalars with unbolded lower-case ( ). Normalized vectors are denoted with a hat (x). We also use the tensor unfolding operator Vec(·) [Golub and Van Loan 2013; Kim and Eberle 2022].
In collision processing, length is the most basic geometric computation, because it asks: are two primitives close to each other? We characterize this with a generic signed length function (e , e , e ) = e · (e × e ) The challenge is then to define the e , , edges appropriately. In order to arrive at a unified analysis of both vertex-face and edge-edge collision energies, we will define the edges in terms of a projected collision frame that contains all the important geometric information for collisions, while discarding everything else. This frame will then form the basis for our analysis.

Vertex-Face and Edge-Edge Lengths
For vertex-face collisions, we denote e 0 and e 1 as the two edges of the triangle face, and e 2 as the vector from any vertex on the triangle to the collision vertex (Fig. 2, left). With these edges defined, (e 0 , e 1 , e 2 ) gives us the desired point-plane distance. We observe that {e 0 ,e 1 , e 2 } do not generally form an orthogonal frame. This is most visible with e 2 in Fig. 2 (left), which contains components that are tangential to the triangle face. This extraneous information is not required to compute point-to-plane distance.
For two edges in collision, we denote e 0 and e 1 as the two nearby edges and e ab as the vector of least distance between them (Fig. 2, right). This is usually found with a geometric test [Rhodes 2001] and the signed length is computed as (e 0 , e 1 , e ab ). Again, since {e 0 ,e 1 , e ab } is not generally orthogonal, we see that the edges contain extraneous information.

Our Reduced Frame
In order to distill the problem down to its fundamental geometry, we propose a reduced collision frame, which we denote C = e 0⊥ e 1⊥ e 2⊥ . Here, e 0⊥ , e 1⊥ , and e 2⊥ are orthogonalized versions of the edges from the above collision cases. This frame is distinct from a QR factorization, because we orthogonalize the vectors but do not normalize them, as that would discard important geometric information. This frame closely mirrors the deformation gradient commonly used in the finite element analysis [Bonet and Wood 2008], where irrelevant translation modes are discarded from the analysis by taking the deformation derivative.
The reduced vertex-face and edge-edge collision frames can then be formulated as: where , , and are modified Gram-Schmidt orthogonalization constants. These edges are then sent to (Eqn. 1).
One immediate objection to this approach is that , , and depend on the edges, so the complexity of the derivatives could potentially explode when computing forces. This may be true for general functions, but not for the case of collisions. We will see in §4.3 that, aside from an easy-to-filter term, all potential sources of complexity when calculating gradients and Hessians of energy functions go to zero.

AN ANALYSIS OF COLLISION ENERGIES
We can now use our reduced collision frame to analyze signed length, and from there, generic collision energies.

Existing Analysis Methods
We begin by observing that existing analysis methods are insufficient. We tried to adapt existing deformation gradient-based methods ], but it only yielded an analysis for unsigned lengths = ∥e 2⊥ ∥. This is insufficient because many collision energies use the sign to determine whether a vertex is approaching (positive) or already in collision (negative) with a face.
We show the results of this analysis in Appendix A. In general, they use matrix invariants along the lines of 3 (C) = det C and 5 (C, v) = v ⊤ C ⊤ Cv. For unsigned length specifically, the cross-terms that arise from differentiating these invariants are numerous and nontrivial, and result in rank-one updates that necessitate numerical root-finding. Overall, we did not find this to be a promising path towards an analytic eigensystem.

Signed Length Analysis
Instead, we will analyze a generic collision energy that depends on the length between two primitives Ψ( ), which we abbreviate Ψ . Our goal is to gain a detailed understanding of the derivatives of this function with respect to our reduced collision frame C. In §4.3.1, we will show how to transform these generic results back to the vertex-face and edge-edge cases.
In particular, we seek the eigensystem of the Hessian of Ψ : where c = Vec (C), H = 2 c 2 , and g = c . We will see later on that g resides in the nullspace of H , so the main challenge is to find the eigensystem of H .
In contrast to the matrix invariant approach, we have found that direct inspection of the Hessian reveals its underlying structure. We directly examine H : where the blocks are composed of second derivatives of : A detailed derivation of these blocks is given in the supplemental materials. The matrix is symmetric, so we have fully described its entries. We observe that all of the individual blocks are either outer product or projection terms, which allows us to state the eigensystem of H directly: where 0 3 and 0 3×3 respectively denote a vector-3 and 3 × 3 matrix of zeros, and = 0,1,2,3 0,1,2,3 − ∥e 2⊥ ∥ 2 Thus, while H ∈ R 9×9 , we can now see that it is rank-4, and contains a rank-5 nullspace. c 2 . The q 0 and q 1 vectors respectively represent buckling and rotation in the direction of e 1⊥ , while q 2 and q 3 are analogous modes along e 0⊥ , and q 4 is a purely compressive mode.
These simple expressions are possible because of the orthogonality of our reduced frame. Without it, Eqns. 5-7 balloon in complexity, making an analytic eigensystem intractable. For more details on how we derive H , see the supplement.
With the eigensystem (Eqns. 8 and 9) in hand, we can return to the generic problem of 2 Ψ c 2 (Eqn. 3). In that equation, the Ψ term simply scales the eigenvalues of H . Additionally, is orthogonal to H , because its sole entry e 2⊥ is known to be orthogonal to the corresponding e 0⊥ or e 1⊥ terms in q 0,1,2,3 . Thus, g must form a fifth eigenpair whose eigenvalue is 2 Ψ 2 . The complete rank-5 eigensystem for generic collision energies is then: 4.2.1 Discussion. Physically, q 0 , q 2 , and q 4 correspond to compression modes in the reduced frame, while q 1 and q 3 are rotation modes. Eigenvalue filtering can then be interpreted as disambiguating between two particular rotation and compression modes, along with filtering compression in the normal direction (Fig. 3).

Filtering in State Space
Our generic analysis can be applied to vertex-face collisions by introducing an appropriate changeof-basis. If we denote the stacked vector of vertex positions v as x ∈ R 12 , the force and stiffness densities can respectively be written using the chain rule: where : is a tensor double contraction [Kim and Eberle 2022]. Our eigenanalysis of 2 Ψ c 2 now seems to be insufficient, because the double contraction term appears to add a difficult new term that must also be analyzed. However, the orthogonality properties of C allow us to express everything in terms of 2 Ψ c 2 . By treating the orthogonalization scalars , , and as though they were constant, we decompose c x into the sum of a simple term c x and a remainder term c x Δ . The gradient and Hessian can then be written in the following forms: This transformation is shown in detail in Appendix B.1. These expressions avoid the need to compute the unwieldy 3rd order tensor 2 c x 2 ∈ R 9×12×12 , and an analytic filtering strategy immediately becomes apparent. Equipped with the analytic eigensystem for 2 Ψ c 2 , we can pin negative eigenmodes to zero in the first term in Eqn. 17, and pin the positive eigenmodes to zero for the second term, resulting in a PSD matrix.
Replacing c with c vf or c ee in Eqns. 16 and 17 yield the full gradients and Hessians.

Signed Length Energies
Once a penalty energy has been rewritten in terms of signed length , Eqns. 11-13 can be used to obtain the analytic eigensystem of the Hessian for both the edge-edge and vertex-face case. Letting be the activation distance and be a force constant, several common penalty energies can be rewritten into this form.
4.4.1 Quadratic Penalty. Inserting a stiff penalty spring at a small distance threshold prevents interpenetration, and can reduce the system sizes needed for constraint-based collision resolution [Bridson et al. 2002;Harmon et al. 2008]. The quadratic energy in terms of signed length is: Fully working out the analytic eigensystem, we have where and ( , ) are as defined in Eqn. 10. Here, we see that the fourth eigenmode corresponds precisely with simple harmonic motion in the direction of e 2⊥ and is always positive. Analyzing signs, we additionally note that ( , ) is always positive and greater than 1, meaning that one of 0,1 and 2,3 is always negative due to the multiplication by (1 ± ( , )). As a result, one eigenmode of each pair always gets filtered. Furthermore, since < is guaranteed, only the sign of determines which vectors should be filtered between them.
We also make use of an analogous quadratic penalty energy using unsigned length Ψ uq : where is a regularization constant to avoid singularities in the square root under differentiation. By applying the general equations, the eigensystem is as follows: Looking closely at the signs, filtering between q 0,1 and q 2,3 once again comes down to the sign of , while 4 is always positive.

Barrier Energies.
IPC solvers use barrier energies to simulate a wide variety of collisions [Li et al. 2021]. These can be rewritten as: The analytic eigensystem is then Interestingly, 0 , 2 , and 4 are always positive while the other eigenvalues are always negative, meaning that filtering the barrier energy does not require any sign checks. However, since barrier energies are undefined when signed length is negative, they become unsuitable for solvers that are expected to deal with unavoidable interpenetration, as with the scripted intersection of kinematic objects [Baraff et al. 2003;Daviet 2020].

Hybrid Energies
McAdams et al.
[2011] characterize stiff springs using an parameter to control the anistropy of the response, observing that higher anisotropy enables more tangential sliding. Reformulating in terms of signed and unsigned length, we see that this tangential sliding corresponds to the difference in eigensystems between signed and unsigned length Hessians. For signed length, any tangential motion comes with a corresponding motion in one of the face vertices (Fig.  3). For unsigned length, pure tangential motion make up two of the three eigenvectors (Fig. 8).
If we characterize hybrid length as ℎ = · + (1 − ) · , with 0 ≤ ≤ 1 as the anisotropy parameter, the mixed energy becomes and the Hessian becomes: where H = 2 c 2 . Our analysis gives us the full eigensystem for the sum of the first two terms with Eqns. 11-13, with the slight modification of 4 = 2 Ψ 2 + 2 Ψ 2 + 2 2 Ψ . Invariant analysis in Appendix A further reveals the eigenpairs of the final term (Eqn. 52). These terms can thus be filtered separately and then summed, yielding an analytic filtering strategy.

Matlab Verification
To ensure that our analyses are correct, we have provided Matlab implementations of key formulas, along with numerical tests. Verify_General_Length.m verifies the eigensystems for H (Eqns. 8 and 9) and H (Eqns. 51 and 52) by generating random reduced frames and comparing the analytic outer-product decomposition with finite differences. We also provide Hessian eigensystems for Ψ lq and Ψ uq with the same numerical test, along with numerical verification for Eqns. 16 and 17 for each case of signed/unsigned length and vertexface/edge-edge collision (see Verify_[lq/uq]_Collision_[vf/ee].m).

Cable Stretching.
Beginning with an implicit strand solver [Bergou et al. 2010], we pull two collections of cables against each other by moving their endpoints in opposite directions and resolve collisions with various penalty energies. We used an implicit solver [Baraff and Witkin 1998] with a large time step of Δ = 1 /30. To better view the unadorned stability behavior of the collision energies, no line search was used.
We initialized the cables with a Young's modulus = 10 4 Pa and a Poisson's ratio = 0.36. All quadratic penalties had the same spring constant = 10 N/m and distance envelope = 0.09 m. We set a gravitational force to induce an acceleration of 9.81 m/s 2 . To prevent excessive oscillations, we also incorporated Rayleigh damping with a factor = 0.006.
We see that Gauss-Newton approximations of the Hessian discards important information. When equipped with Ψ uq , the jittering in Gauss-Newton causes cables under tension to unrealistically pop through (Fig. 1). Analytically clamping the same scene correctly resolves the contact.
Stretching the cables again with Ψ lq , Gauss-Newton persistently contains jittering artifacts. The jittering diminishes when using our approach, but still remains present. However, these artifacts arise from the collision energy itself, not the filtering strategy. When we use a hybrid energy along the lines of Cheng et al. [2022], which falls back to Ψ uq when the closest point on the edge corresponds to an endpoint, we obtain less jittery results. Simply passing in the unfiltered Hessian with the same scene parameters results in an explosion immediately after the first few contacts.
The energies exhibited qualitative differences that correspond to their formulations. Unsigned penalty forces that do not explicitly cancel off tangential sliding were "stickier" than their signed counterparts (Fig. 4).

Fabric Stretching.
We ran a test similar to the cable example with an implicit cloth solver [Baraff and Witkin 1998], Δ = 1 /40, and no line search, pulling two pieces of cloth against each other (Fig. 5). We initialized the fabrics with an ARAP membrane stretching energy ] with stretch = 10 and a dihedral bending energy with bend = 0.01. All quadratic penalties  were set with a spring constant = 4000 N/m and a distance envelope = 0.02 m. Gravity was set to be 1 m/s 2 . To prevent excessive oscillations, we used Rayleigh damping with = 1.
Picking between Ψ lq or Ψ uq for vertex-face and edge-edge collisions results, we simulated two cloth pulling tests. As seen in the video, Gauss-Newton approximation results in jittery cloth across the board, while our analytic clamping smoothly settles.

Hammock Drop.
We evaluated the qualitative behavior of a variety of energy combinations in a more complex collision scenario (Fig. 6) where several patches of cloth are dropped onto a hammock and then slide off. The differences in "stickiness" between energies is visible in the way that the cloth bunches and the rates at which different patches fall away. The simulation parameters were the same as §5.2.2, except that Rayleigh damping was set to zero.

Volumetric Meshes.
We showcase the necessity of negative eigenvalue filtering in a volumetric example, where we lower a set of volumetric rods onto a collection of thicker rods with a time step Δ = 1/40 s and no line search.
We used a Stable Neo-Hookean distortion energy [Smith et al. 2018] with = 10 4 Pa and = 0.45. We set gravity to 1 m/s 2 as well. For collision response, we used Ψ lq and Ψ uq for the vertex-face and edge-edge collision penalties, respectively. Both were initialized with a spring constant = 1200 N/m and a distance envelope = 0.02 m. Rayleigh damping was at = 0.25.
Disabling eigenvalue filtering in Ψ lq and Ψ uq results in indefinite systems that lead to obvious artifacts around frame 170 and a subsequent explosion (Fig. 7). All the other unfiltered scenes exhibited similar behavior. The specific vertex-face (VF) and edge-edge (EE) energies used are listed along the left. From left to right, the "stickiness" that the collision energy gives the red cloth determines whether the blue cloth falls away at frame 240. Similarly, the stickiness determines whether the red cloth remains by frame 410 (right). We analytically filtered all of the collision energies, and ran the simulations at Δ = 1 /40 s These observations are in line with previous work, which stress the necessity of eigenvalue filtering for PSD-projection not only for collisions [Li et al. 2020], but also several other deformation energies Lin et al. 2022;Smith et al. 2019;Teran et al. 2005].

Performance
We initialized 10 7 random vertex arrangements and compared the performance of our analytic approach with the numerical eigensolver solver in Eigen [Guennebaud et al. 2010] and Gauss-Newton approximation. Our technique is 4.5× -7.2× faster than the numerical eigensolver, and requires only 8% − 70% more computation than the single outer product for Gauss-Newton. As shown in Fig. 1, this modest overhead significantly reduces jittering and pop-through.

CONCLUSIONS AND FUTURE WORK
We have presented an eigenanalysis of energies used for vertex-face and edge-edge collisions in implicit solvers. By using analytic methods to filter negative eigenvalues, we obtain faster matrix assembly times and more robust simulations. The analytic eigenpairs represent distinct collision modes, and clarify how filtering certain eigenvalues correspond to discounting specific responses. We could extend our analysis to vertex-vertex and vertex-edge, mollified energies in near-parallel, and other degenerate cases, allowing for analytic eigensystems in more hybridized penalty energies. Potential-based frictional contact [Li et al. 2020] and damping [Brown et al. 2018]   are other interesting directions, since the energies we analyzed only penalize interpenetration, not relative sliding motion.

ACKNOWLEDGMENTS
We thank the reviewers for their comments. This work was supported by Adobe, the Bungie Foundation, Teng and Han Family Fund, and National Science Foundation (IIS-2132280). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

A INVARIANT-BASED ANALYSIS
The unsigned length and signed length of a point translated e 2⊥ from the plane whose normal is in the direction of e 0⊥ × e 1⊥ are: = ∥e 2⊥ ∥ (37)

A.2 Unsigned Length Energies Are Special Cases
The eigenanalysis from  can be applied to any energy of the form Ψ ( 5 (v)) where v is kept constant. Energies that are functions of only fit these requirements, since Ψ ( ) = Ψ √︁ 5 (ẑ) . The eigenpairs of 2 Ψ c 2 are then:   8. Eigenvectors of Ψ : q 0 represents compression, while q 1,2 represent sliding in the parallel plane.
1,2 simplify similarly. As for v T 1 ,T 2 , the reduced frame consists of orthogonal vectors, so we take V as the identity, with U and as Plugging these into v 1 and v 2 yields simpler expressions for the eigenpairs of unsigned length energies: Physically, v 0 represents spring compression along the normal, while v 1,2 represent buckling in the orthogonal directions (Fig. 8).
In practice, we only apply Eqns. 51-52 for stiff springs using unsigned length in the vertexface collision case, while we use Eqns. 41-42 for the edge-edge case. This is because a simpler non-orthogonal frame unique to the edge-edge case exists: where e ab = e 2 − e 0 + e 1 is the vector of least distance between e 0 and e 1 , and and are clamped and taken as constant, matching Kim and Eberle [2022]'s implementation. Using this frame, 2 c x 2 is zero, circumventing the need for any extra term considerations.