 # Sparsity Exploitation of Accelerated Modulus-Based Gauss-Seidel Method for Interactive Rigid Body Simulations

Large-scale linear complementarity problems (LCPs) are repeatedly solved in interactive rigid-body simulations. The projected Gauss-Seidel method is often employed for LCPs, since it has advantages in computation time, numerical robustness, and memory use. Zheng and Yin (2013) proposed modulus-based matrix splitting iteration methods and showed their effectiveness for large problems, but a simple application of their approach to large-scale LCPs in interactive rigid-body simulations is not effective since such a simple application demands large matrix multiplications. In this paper, we propose a novel method derived from accelerated modulus-based matrix splitting iteration methods that fits LCPs arising in interactive rigid-body simulations. To improve the computation time, we exploit sparsity structures related to the generalized velocity vector of rigid bodies. We discuss the convergence of the proposed method for an important case that the coefficient matrix of a given LCP is a positive definite matrix. Numerical experiments show that the proposed method is more efficient than the simple application of the accelerated modulus-based Gauss-Seidel method and that the accuracy in each step of the proposed method is superior to that of the projected Gauss-Seidel method.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In rigid-body simulations, interactions (for example, normal forces) between rigid bodies are often mathematically modeled as constraints. For computing these constraint forces, we usually need to solve certain equations. Two representative constraint formulations for rigid-body simulations are acceleration-based formulations  and velocity-based formulations . In the acceleration-based formulations, the constraints are described with forces and accelerations of rigid bodies; we first compute forces and accelerations, then integrate them to obtain velocity changes. On the other hand, in the velocity-based formulations, the variables in the constraints are impulses and velocities of the rigid bodies. In this paper, we focus on the velocity-based formulations, since the velocity-based formulations are widely used and are known to be superior to the acceleration-based formulations in many aspects (for example, see [6, 8]).

There are two main categories for solving constraints, iterative approaches and direct approaches. We focus on the iterative approaches rather than the direct approaches, since the direct approaches such as pivoting methods often suffer from time complexity and numerical instability as pointed in . In the impulse-based iterative approaches , impulses are applied to the rigid bodies sequentially, until certain convergence conditions are satisfied.

Contact constraints are frequently modeled in a form of complementarity problems ; in particular, linear complementarity problems (LCPs) give mathematical formulations for frictionless contacts. Among various iterative methods for solving LCPs, the projected Gauss-Seidel (PGS) method  has a remarkable flexibility, therefore the PGS method or its extensions are often employed to solve contact constraints [6, 8, 7]. Another iterative approach for solving LCPs is the use of modulus-based methods. Bai  established modulus-based matrix splitting iteration (MMSI) methods, which include modulus-based Jacobi (MJ), modulus-based Gauss-Seidel (MGS), modulus-based successive over relaxation (MSOR), and modulus-based accelerated overrelaxation (MAOR) iteration methods as special cases. Recently, Zheng and Yin  proposed an accelerated modulus-based matrix splitting iteration (AMMSI) method as an improvement of Bai . In a similar way to the MMSI methods, the AMMSI methods include accelerated modulus-based Jacobi (AMJ), accelerated modulus-based Gauss-Seidel (AMGS), accelerated modulus-based SOR (AMSOR), and accelerated modulus-based accelerated overrelaxation (AMAOR) iteration methods. To our best knowledge, no application of the MMSI or AMMSI methods to real-time simulations has been examined before.

Since the AMGS method is not designed for interactive rigid-body simulations, a simple application of the AMGS method causes inefficiency in computation and it is a serious disadvantage for real-time simulations. In many applications of real-time simulations, interactive computer graphics and operations are considered most important. If the computation in each simulation step is considerably slower than real time, the quality of the users’ experience would be seriously degraded.

In this paper, we resolve this difficulty by focusing the update formula in the AMGS method and exploiting the structures related to the generalized velocity vector of rigid bodies.

We also give a theoretical proof on the convergence of the AMGS method. Bai  already discussed the convergence, but the assumption in  is too restrictive to apply the same discussion to rigid-body simulations. We extend the proof of  to cover an important case that the coefficient matrix in the LCP is a positive definite matrix.

Through numerical experiments, we observed that the proposed AMGS method that exploited the sparsity attained shorter computation time than the original AMGS method. Furthermore its convergence rate in each iteration was better than that of the PGS method. These results indicate that the proposed method is useful for practical real-time simulations.

The outline of this paper is as follows. In Section 2, we briefly introduce a formulation of velocity-based constraints as an LCP. We also discuss the projected Gauss-Seidel method and the AMGS method to solve LCPs. The application of the AMGS method to rigid-body simulations is developed in Section 3, and we prove convergence theorems of the AMGS method in Section 4. In Section 5, we will show numerical results to verify the efficiency of the AMGS method. Finally, we will give a conclusion in Section 6.

## 2 Preliminaries

### 2.1 Linear complementarity problem with velocity-based constraints

For the latter discussions, we briefly introduce an LCP that arises from velocity-based constraints. For more details, the readers can refer to [4, 9, 10].

During a rigid-body simulation, we keep tracking movements of the rigid bodies in multiple time periods, therefore, an entire simulation is divided into a sequence of simulation steps, and each simulation step corresponds to a small time step. Since the constraints on the rigid bodies should be satisfied at each time, we solve the following LCP in each simulation step:

 (1)

In this paper, we use the superscript to denote the transpose of a vector or a matrix.

The decision variable in this LCP is which is the impulse vector applied to the rigid bodies in the constraint space. The first constraint in (1) requires to be nonnegative to ensure that the constraint impulse must be repulsive.

The second constraint in (1) corresponds to the velocity constraints in the rigid-body simulation. The vectors and

are the bias vector in a constraint space and the generalized velocity vector of rigid bodies, respectively. More precisely, when we have

rigid bodies, is a vector that consists of linear and angular velocities, i.e.,

 v=(vT1ωT1vT2ωT2⋯vTNωTN)T

where and are the linear velocity and the angular velocity of the th rigid body, respectively, for ; thus the length of is . The matrix is the Jacobian matrix corresponding to the velocity constraints. The generalized mass matrix of the rigid bodies

consists of masses and inertia tensor matrices in the diagonal positions:

 M=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝m1E3I1m2E3I2⋱mNE3IN⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (2)

where and are the mass and the inertia tensor matrix of the th rigid body, respectively. We also use

to denote the identity matrix of order

. The inertia tensor matricies are symmetric, so is .

The third constraint in (1) is a complementarity condition. We can understand this complementarity condition as follows. If holds for some , then the rigid bodies are moving away from each other in the direction of the th constraint, therefore, the th constraint should be “inactive”. However, is the impulse vector, thus implies that the th constraint must be “active”. Hence, and should not hold simultaneously, and this requirement is implemented in the complementarity condition.

By denoting and and introducing an auxiliary variable , the LCP (1) can be expressed in a general LCP as follows:

 LCP(q,A)⎧⎪⎨⎪⎩Aλ+q=wwTλ=0λ,w≥0 (3)

It is known that if the constraints are non-degenerate, the Jacobian matrix is full row rank and the matrix is positive definite (see , for example). Throughout this paper, we assume that is positive definite.

At the end of this subsection, we should note that the input data and vary in accordance with simulation steps. If we express the time dependence explicitly, they should be and where is the simulation step. However, in this paper, we mainly focus on solving (1) in each simulation step, therefore, we usually drop the simulation step from and .

### 2.2 Projected Gauss-Seidel method

In the LCP (3) from the rigid-body simulation, the matrix has a structure such that . The projected Gauss-Seidel (PGS) method  is designed to solve more general LCPs (3) in the sense that the assumption for is only positive definiteness. The PGS method is an iterative method and generates a sequence .

A key property in the PGS method is to decompose into such that is a diagonal matrix, a strictly lower triangular matrix, and a strictly upper triangular matrix. Since we assume is a positive definite matrix, is invertible. Due to this decomposition, is equivalent to

Taking this formula and the complementarity condition into consideration, the PGS method computes the next iteration by the following update formula:

 λk+1=max{0,D−1(Lλk+1+Uλk−q)}. (4)

Throughout this paper, we use () to denote the element-wise maximum (minimum, respectively) of two vectors and . The PGS method continues the update by (4) until the sequence converges enough, or the number of the iterations reaches a certain limit.

In the rigid-body simulation, the initial vector is usually set as a zero vector or the impulse vector obtained in the previous simulation step. Since the initial vector often affects the performance of iterative approaches, the use of the solution from the previous simulation step makes the convergence faster . Such a technique is called warm start.

### 2.3 Accelerated modulus-based Gauss-Seidel method

For solving the general LCP (3), Bai  devised the following implicit fixed-point equation that is essential for the modulus-based matrix splitting iteration (MMSI) methods:

 (M0Γ+Ω1)x=(N0Γ−Ω2)x+(Ω−AΓ)|x|−q (5)

Here, the matrices and are a splitting pair of such that . and are two diagonal matrices whose diagonal entries are positive. and are nonnegative diagonal matrices such that . We should emphasize that the variable in (5) is , and we use to denote the element-wise absolute values of . The relation between and the pair of and in (3) will be discussed in Theorem 2.1.

By setting , , and with a parameter in (5), we can obtain a simplified implicit fixed-point equation:

 (M0+γΩ)x=N0x+(γΩ−A)|x|−γq. (6)

Based on (6), the iteration of the MMSI method can be derived as follows:

 (M0+γΩ)xk+1=N0xk+(γΩ−A)|xk|−γq. (7)

We decompose into in the same way as the PGS method. Set , and let and be two parameters. Then, we can derive the update formula of four methods from (7); the modulus-based Jacobi (MJ) by setting in (7), the modulus-based Gauss-Seidel (MGS) by , the modulus-based successive over relaxation (MSOR) by , and the modulus-based accelerated overrelaxation (MAOR) iteration method by , respectively.

Zheng and Yin  utilized two splitting pairs of the matrix such that , and devised a new equation based on (5):

 (M1Γ+Ω1)x=(N1Γ−Ω2)x+(Ω−M2Γ)|x|+N2Γ|x|−q. (8)

Zheng and Yin  established the following theorem to show an equivalence between (8) and in (3). Since a detailed proof is not given in , we give the proof here.

###### Theorem 2.1.

 The following statements hold between (8) and :

1. if is a solution of , then satisfies (8).

2. if satisfies (8), then the pair of and is a solution of .

###### Proof.

We first prove (i). Since is a solution of , satisfies the four constraints, , , and . The first constraint is equivalent to

 (Ω+AΓ)(Γ−1λ−Ω−1w) =(Ω−AΓ)(Γ−1λ+Ω−1w)−2q.

From the rest three constraints and the fact that and are diagonal matrices whose diagonal entries are positive, if , it holds that . Therefore, satisfies

 (Ω+AΓ)x =(Ω−AΓ)|x|−q (9)

and this is equivalent to (8).

To prove (ii), from (9), it holds that . By the relations and , we obtain . Since and are positive diagonal matrices, it is easy to check that and are nonnegative vectors. Finally, it is also easy to show the element-wise complementarity between and . ∎∎

We may use Theorem 2.1 to establish some iterative methods for solving , but we need to set appropriate matrices for the implicit fixed-point equation (8) in actual computations. In particular, the splitting pair of is not unique. By fixing , and , we derive a simplified update equation of (8) as follows:

 (M1+γΩ)x=N1x+(γΩ−M2)|x|+N2|x|−γq. (10)

Based on this equation, Zheng and Yin  provided an update formula of the AMMSI methods:

 (M1+γΩ)xk+1=N1xk+(γΩ−M2)|xk|+N2|xk+1|−γq (11)

When the sequence converges enough, the AMMSI methods output the impulse vector by using the relation .

By changing the splitting pairs of , the update formula (11) above yields variant methods; MMSIM ( and ), the accelerated modulus-based Jacobi (AMJ) iteration method (, , and ), the accelerated modulus-based SOR (AMSOR) iteration method (, , and ), and the accelerated modulus-based accelerated overrelaxation (AMAOR) iteration method (, , and ).

In particular, the update formula of the accelerated modulus-based Gauss-Seidel (AMGS) method in  is derived with , , and as follows:

 (D+γΩ−L)xk+1=Uxk+(γΩ−D+U)|xk|+L|xk+1|−γq. (12)

Let be the difference between and . Then, (12) is equivalent to

 (D+γΩ)Δxk=Lxk+1−(γΩ+D−U)xk+(γΩ−D+U)|xk|+L|xk+1|−γq. (13)

By Theorem 2.1 and , the sequence for the LCP (3) can be associated with the sequence generated by (12) by the relation , thus is a multiple of the positive part of . This motivates us to split into the positive and negative parts such that , where and . From the relations and , (13) is equivalent to

 (D+γΩ)Δxk=γLλk+1−(γΩ+D−U)(γ2λk−xk−)+(γΩ−D+U)(γ2λk+xk−)−γq=γLλk+1−γ(D−U)λk+2γΩxk−−γq.

Therefore, for computing , we only need the th component of , which will be denoted as , since and are diagonal matrices. This simplifies the computation of . Recalling the decomposition of , we compute for each by

 Δxki =γDii+γΩii(i−1∑j=1Lijλk+1j−Diiλki+m∑j=i+1Uijλkj+2Ωii(xk−)i−qi) =−γAii+γΩii(i−1∑j=1Aijλk+1j+m∑j=iAijλkj+qi−2Ωii(xk−)i). (14)

We can summarize a framework of the AMGS method as follows.

###### Method 2.2.

(the AMGS method for (3))

Choose a nonnegative vector as an initial vector. Generate the iteration sequence by the following procedure:

repeat
for  do

.
until  satisfies a certain convergence threshold

## 3 Accelerated modulus-based matrix splitting iteration methods for interactive rigid-body simulation

In this section, we propose a numerical method to solve LCPs that arises from interactive rigid-body simulations using the AMGS method. Since the AMGS method is a generic method for LCPs, it is possible to simply apply the AMGS method to (3). However, explicit evaluation of is inefficient even though the matrices and are sparse. Thus, such a simple application of the AMGS method is not practical. As we will discuss in later, this is mainly because the number of the non-zero elements in becomes very large. To overcome this inefficiency, we modify the AMGS method so that it does not require the explicit evaluation of the matrix . In a similar way to the AMGS method, we can also modify the AMSOR method for solving LCPs in the rigid-body simulations.

As already pointed out at above, the direct computation of is unfavorable for real-time simulations. In the viewpoint of the computation cost, we should avoid the computations of and , which involve all the off-diagonal elements of .

To improve the computation efficiency, we introduce an intermediate variables and . In particular, stores information of applied impulse [6, 10]. For , let

 λk+1,i =(λk+11…λk+1i−1λki…λkm)T (15) vk+1,i =v+ˆMλk+1,i (16)

where By the definitions of and , it holds

 i−1∑j=1Aijλk+1j+m∑j=iAijλkj+qi=(Aλk+1,i)i+qi = (JM−1JTλk+1,i)i+(Jv)i+bi=(Jvk+1,i)i+bi=m∑j=1Jijvk+1,ij+bi,

 xk+1i=xki+Δxki=xki−γAii+γΩii(m∑j=1Jijvk+1,ij+bi−2Ωii(xk−)i).

Due to the relation , we can compute by updating only the th position of ,

Thus, from (16), we obtain

 vk+1,i+1=vk+1,i+ˆM(λk+1,i+1−λk+1,i)=vk+1,i+ˆM∗i(λk+1,i+1i−λk+1i),

where is the th column of .

The following method summarizes the proposed AMGS method for rigid-body simulations.

###### Method 3.1.

(the proposed AMGS method for rigid-body simulations)

Choose a nonnegative vector as an initial vector. Generate the iteration sequence by the following procedure:

.
repeat
for  do

if  then

else

until  satisfies a certain convergence threshold

We compare the computation costs of Method 2.2 and Method 3.1. In Method 2.2, we update , and for each , we compute and in (14). Therefore, if is fully dense, the computation cost to obtain is .

In Method 3.1, we can exploit the structures of and . We use to denote the number of nonzero elements in . Since the matrix is composed with and at the diagonal positions as shown in (2), we know that . Actually, each row of is a multiple of one column of or a linear combinations of three columns of . To obtain from , the computation of is required for , thus the computation cost in the th outer iteration amounts to . Similarly, to obtain from , we need , thus the computation cost of this part for each is , and this is same as . Consequently, the computation cost to obtain from in Method 3.1 is . In the rigid-body simulation, is and is smaller than , thus we can expect Method 3.1 is much faster than the direct use of in Method 2.2. We will verify this efficiency in the numerical experiments of Section 5.

We should mention that the matrix is not always fully-dense, and we can actually find some zero elements in of the test instances that will be used in Section 5. However, the positions of nonzero elements cannot be determined before the multiplication . In addition, even if we skip the zero elements of in Method 2.2, the effect is less significant than Method 3.1 and the numerical efficiency of Method 2.2 is still insufficient for real-time simulations.

### 3.1 Linear Complementarity Problems with Lower and Upper Bounds

In this section, we discuss a method for solving LCPs with lower and upper bounds on , which often occur in the formulations of contact constraints with friction [10, 8]. We call such LCPs with lower and upper bounds “Boxed LCPs (BLCPs)”. Consider the following BLCP:

 BLCP(q,l,u,A)⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩Aλ+q=wl≤λ≤ufor each i=1,…,m,⎧⎨⎩wi≥0if \ λi=liwi≤0if \ λi=uiwi=0if \ li<λi

Here, and are the lower and the upper bounds, respectively. Without loss of generality, we assume for each .

We define a projection function of to the interval and by

 pBLCP(λ)=min{max{l,λ},u}.

Since is a nonnegative vector, so is . With this projection, we can give a AMGS method for Boxed LCPs as follows.

###### Method 3.2.

(an AMGS method for Boxed LCPS in rigid-body simulations)

Choose a nonnegative vector as an initial vector. Generate the iteration sequence by the following procedure:

)
repeat
for  do

if  then

else

until  satisfies a certain convergence threshold

Note that with and being a very large vector represents the same problem as , so Method 3.2 can be considered as a generalization of Method 2.2.

## 4 Convergence theorem

In this section, we focus on the convergence of the accelerated modulus-based Gauss-Seidel (AMGS) method. Zheng and and Yin  discussed only the two cases: (i) is a positive definite matrix and is a multiple of the identity matrix ( for some ), and (ii) is an -matrix.

Since the matrix is always positive definite in the rigid-body simulation, we extend the proof in  so that we can handle a more general form of . For example, if we can take , we may be able to improve the convergence, as mentioned later in Remarks 4.4 and 4.5. Here, is a positive constant, and is the diagonal matrix whose diagonal elements are those of , thus is not covered by .

In order to establish the convergence theorem of the AMGS method with which is not always a multiple of the identity matrix, we need Lemma 4.1 below. We use to denote the Euclidean norm of . We also use to denote the spectral norm of a matrix , thus we can employ . In addition let denote the diagonal matrix whose diagonal elements correspond to those of a matrix , i.e.,

 diag(X)=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝X11X22⋱Xmm⎞⎟ ⎟ ⎟ ⎟ ⎟⎠.
###### Lemma 4.1.

Let be lower triangular matrices, and be a non-singular lower triangular matrix. Then, the following equations hold:

 ∥L1∥ =∥diag(L1)∥ (17) diag(L1L2) =diag(L1)diag(L2) (18) diag(ˆL−1) =(diag(ˆL))−1 (19)
###### Proof.

We first prove (17). As

is a lower triangular matrix, the eigenvalues of

are , and these values are also those of the diagonal matrix . Since the spectral norm of a matrix is the maximum absolute value of its eigenvalues, we have .

Next, we consider (18). Focusing the th diagonal element of , we have

 (L1L2)ii =m∑j=1(L1)ij(L2)ji =∑j∈{1,2,…,i}∩{i,i+1,…,m}(L1)ij(L2)ji =(L1)ii(L2)ii,

therefore, it holds that .

Finally, we prove (19). By (18) and the fact that the inverse matrix of a lower triangular matrix is also a lower triangular matrix, we obtain

 diag(ˆL)diag(ˆL−1)=diag(ˆLˆL−1)=Em.

Therefore, is the inverse matrix of . ∎

We are now ready to establish the convergence theorem. The theorem covers the case that , which will be actually used in numerical experiments of Section 5. For the subsequent discussion, we define .

###### Theorem 4.2.

Let

 δ=2∥(D+γΩ)−1D∥τ+∥(D+γΩ)−1(γΩ−D)∥.

If , the iteration sequence generated by Method 3.1 with an arbitrary nonnegative initial vector converges to the unique solution of .

###### Proof.

In the AMMSI methods  from which the AMGS was derived, we chose and for the simplified fixed-point equation (10). Let be a solution of , and let . From Theorem 2.1, satisfies (10), and the convergence of to can be guaranteed by that of to . Therefore, we focus the convergence of the sequence . Subtracting (10) with from the update formula (11), we obtain

 (M1+γΩ)(xk+1−x∗)= N1(xk−x∗)+(γΩ−M2)(|xk|−|x∗|) +N2(|xk+1|−|x∗|).

In the setting of the AMGS method ( and ), we can further evaluate this inequality as follows:

 (xk+1−x∗)= (D−L+γΩ)−1U(xk−x∗) (20) +(D−L+γΩ)−1(γΩ−D+U)(|xk|−|x∗|) +(D−L+γΩ)−1L(|xk+1|−|x∗|).

Using the triangular inequality for any two vectors and of the same dimension, we have

 ∥xk+1−x∗∥≤ ∥(D−L+γΩ)−1U∥∥xk−x∗∥ (21) +∥(D−L+γΩ)−1(γΩ−D+U)∥∥xk−x∗∥ +∥(D−L+γΩ)−1L∥∥xk+1−x∗∥ ≤ ∥(D−L+γΩ)−1U∥∥xk−x∗∥ +∥(D−L+γΩ)−1(γΩ−D+L−L+U)∥∥xk−x∗∥ +∥(D−L+γΩ)−1L∥∥xk+1−x∗∥ ≤ 2∥(D−L+γΩ)−1U∥∥xk−x∗∥ +∥(D−L+γΩ)−1(γΩ−D+L)∥∥xk−x∗∥ +∥(D−L+γΩ)