A mathematical theory for mass lumping and its generalization with applications to isogeometric analysis

Explicit time integration schemes coupled with Galerkin discretizations of time-dependent partial differential equations require solving a linear system with the mass matrix at each time step. For applications in structural dynamics, the solution of the linear system is frequently approximated through so-called mass lumping, which consists in replacing the mass matrix by some diagonal approximation. Mass lumping has been widely used in engineering practice for decades already and has a sound mathematical theory supporting it for finite element methods using the classical Lagrange basis. However, the theory for more general basis functions is still missing. Our paper partly addresses this shortcoming. Some special and practically relevant properties of lumped mass matrices are proved and we discuss how these properties naturally extend to banded and Kronecker product matrices whose structure allows to solve linear systems very efficiently. Our theoretical results are applied to isogeometric discretizations but are not restricted to them.


Introduction and background
Mass lumping dates from the early days of the finite element method, in the 1970s, when practitioners were facing the burden of inverting non-diagonal mass matrices in explicit time integration schemes for solving wave propagation problems. They quickly discovered that reasonably good solutions could still be recovered after replacing the mass matrix by some easily invertible (typically diagonal) ad hoc approximation, a procedure which became known as mass lumping. The mass lumping approach relied on physical intuition of mass conservation [1,2]. Numerous strategies were proposed but all later attention in modern literature has focused on only three of them: the rowsum technique, the nodal quadrature method and a so-called "special mass lumping" scheme [2], more commonly referred to as the diagonal scaling procedure [3].
The row-sum method was among the first techniques proposed. As the name suggests, the row-sum lumped mass matrix is a diagonal matrix where the ith diagonal element is simply the sum of all entries on the ith row. Apart from being straightforward to implement, it is widely used and appreciated for the improved stability in connection with explicit time integration schemes, a property we will formally prove in this article. However, this lumping technique is prone to failure by generating singular or indefinite lumped mass matrices, leading to erroneous dynamical behaviors. Nonnegative partition of unity methods are exempt from such breakdowns as they always deliver positive definite row-sum lumped mass matrices. Therefore, the emergence of new discretization methods based on a nonnegative partition of unity, such as isogeometric analysis [4,5], has sparked renewed interest in the row-sum technique [6,7]. Unfortunately, the row-sum technique in isogeometric analysis leads to a loss of convergence rate for the smallest generalized eigenvalues. As a matter of fact, it was shown in [6] for 1D problems and for quadratic and cubic B-splines that the convergence rate is downgraded to quadratic order when using the row-sum lumped mass matrix. To this date, high order mass lumping in isogeometric analysis is an open research problem and has only been achieved (to our knowledge) in [8] using dual basis functions.
Another mass lumping technique specific to classical finite element methods is the nodal quadrature method. In this method, the interpolation property of Lagrange basis function is cleverly used to construct diagonal matrices by choosing quadrature nodes that coincide with the finite element nodes. Although this technique might preserve the higher order convergence of the consistent mass [9,10], it also breaks down in case nonpositive quadrature weights are used. Positive weights can be guaranteed by either choosing specific quadrature rules and sacrificing some accuracy [3,11] or enlarging the finite element space by introducing additional degrees of freedom [12]. Both solutions come with their own drawbacks: the first one does not preserve the rate of convergence while the second one requires solving a larger system of equations. In either case, research is hampered by tedious element specific analysis. Other authors, instead, have suggested modifying the time integration scheme to accommodate nonpositive weights. In [13], a mixed implicit-explicit time integration scheme is designed to cope with positive semidefinite (potentially singular) lumped mass matrices.
Finally, the diagonal scaling method [14] was designed to circumvent breakdowns in the row-sum and nodal quadrature methods. Positive definite lumped mass matrices are indeed guaranteed by construction. The method scales the diagonal of the consistent mass matrix, which is necessarily positive, by a positive factor that depends on the mass of the element.
Despite appealing numerical results, the mathematical theory supporting these methods is rather limited and usually restricted to specific discretization methods, orders and element types. Although a completely unified analysis is indeed probably intractable, a framework with the least restrictive assumptions is sought. For this reason, our analysis will be restricted to the row-sum technique because it allows us to draw very general conclusions. It is based on purely algebraic considerations and is therefore independent of the discretization method used. Moreover, although the three aforementioned mass lumping strategies generally yield different lumped mass matrices, they might be identical in some cases. In [3], for instance, equivalence conditions for quadrature rules have been established for the Lagrange basis. Thus, although our conclusions will generally not apply to other mass lumping strategies, they might still hold under certain circumstances.
The rest of the paper is structured as follows: in Section 2, we briefly recall the discretization of a model wave equation and the problem arising from explicit time integration methods for the semi-discrete problem. The generalized eigenvalues of the matrix pencil formed by the stiffness and mass matrices are critical for the dynamics and stability of explicit methods. Section 3 is designed to bring the reader up to date by recalling well-known theoretical results for generalized eigenproblems and paves the way for the upcoming sections. Section 4 is devoted to the theoretical analysis of perturbed generalized eigenproblems. Immediate conclusions can be drawn regarding the stability of explicit time integration methods in connection with mass lumping. Some of our results for high dimensional discretizations assume the mass matrix is expressed as a Kronecker product, a structure inherited by tensor product basis functions. This is generally not the case in practice for a variety of reasons. Therefore, Section 5 explains how the mass matrix may nevertheless be approximated by a Kronecker product. Established results can then be applied to the approximation itself. Numerical experiments are presented along the way and closely follow the theory. They serve multiple purposes by illustrating the findings, highlighting their limitations and drawing additional insight. In order to foster applications to other problems, theoretical results are stated for general matrices satisfying certain properties. Finally, Section 6 concludes our article and summarizes our findings.

Mass lumping in explicit dynamics
We consider as model problem the classical wave equation. Let Ω ⊂ R d be an open connected domain with Lipschitz boundary and let I = [0, T ] be the time domain with T > 0 denoting the final time. We look for u : Ω × [0, T ] → R such that where u 0 and v 0 are some initial displacement and velocity, respectively, and we prescribe homogeneous Dirichlet boundary conditions for simplicity. A Galerkin discretization of the spatial variables (see for instance [2,15]) leads to solving the semi-discrete problem where M, K ∈ R n×n are the mass and stiffness matrices, respectively, and x(t) ∈ R n is a vector containing the coefficients of the approximate solution u h (x, t) in the basis of some finite dimensional subspace. The time-dependent right hand side f (t) ∈ R n accounts for the forcing function f and potential non-homogeneous Neumann and Dirichlet boundary conditions. Naturally, the properties of M and K depend on the discretization method. Nevertheless, they are always symmetric with M being positive definite for the standard Galerkin method. Provided Dirichlet boundary conditions are prescribed on some portion of the boundary, which is necessary for the well-posedness of the problem, K is also positive definite. Therefore, we will assume throughout the paper that both M and K are symmetric positive definite.
There exists a host of numerical methods for approximating the solution of the semi-discrete problem (2.2). The Newmark method [16], proposed in 1959 and summarized in Algorithm 2.1, has been historically used and still remains very popular today for structural problems with fast dynamic loadings such as blasts, impacts or earthquakes. The time domain [0, T ] is first discretized using a uniform time step ∆t = T N such that t s = s∆t, s = 0, 1, . . . , N with N ∈ N * . The method then consists in approximating the solution x(t) as well as its first and second time derivative (interpreted as the displacement, velocity and acceleration, respectively) at the discrete times t s . These quantities are denoted x s , v s and a s in Algorithm 2.1.

7:
Compute x s+1 =x s+1 + β∆t 2 a s+1 8: end for Algorithm 2.1: Newmark method Both explicit and implicit methods are recovered for specific combinations of two scalar parameters β and γ in [0, 1]. Explicit methods are recovered for β = 0 while implicit unconditionally stable methods arise for β ≥ γ 2 ≥ 1 4 . Whereas explicit methods are necessarily conditionally stable, implicit methods may or may not be unconditionally stable. The central difference method for instance is an explicit second order method and its stability analysis reveals that the critical time step is given by ∆t c = 2 √ λn(K,M ) , where λ n (K, M ) is the largest eigenvalue of the generalized eigenvalue problem Ku = λM u. Unfortunately, rapidly changing boundary conditions require a sufficiently small time step for us to resolve the correct physical behavior. Thus, the unconditional stability of some implicit methods cannot be fully exploited while the computational burden prevails in Line 5. Indeed, for very large problems, solving the linear system in Line 5 at each time step becomes a major computational bottleneck. However, for explicit methods (β = 0), this linear system only involves the mass matrix. Therefore, replacing the consistent mass matrix in Lines 1 and 5 with some kind of easily invertible approximationM becomes very appealing. Mass lumping stems from this observation. In fact, in engineering practice all kinds of ad hoc replacements may be made involving also the stiffness matrix or the right-hand side vector but we will not consider those here.
While the computational benefits of mass lumping are obvious, the substitution of M byM both effects the stability of the time integration method and the accuracy of the approximate solution. Strong numerical evidence suggests that mass lumping improves the stability by increasing the critical time step [7,8]. Verifying this claim amounts to guaranteeing that λ n (K,M ) ≤ λ n (K, M ). On the other hand, regarding accuracy, the smallest generalized eigenvalues of (K, M ) usually play a prominent role although the largest ones might also sometimes be involved depending on the initial conditions and loading [17]. Thus, the smallest eigenvalues of (K,M ) and (K, M ) must be very close. Consequently, both questions of stability and accuracy can be tackled by comparing the eigenvalues and eigenspaces of some matrix pencils. The upcoming analysis requires some theoretical understanding of generalized eigenvalue problems. All necessary background knowledge will be covered in the next section.

Generalized eigenvalue problems
Given the important role generalized eigenproblems play in this work, this section provides a brief reminder of some important properties. Whereas standard eigenproblems are clearly a particular case of generalized eigenproblems, results for standard eigenproblems generally cannot be straightforwardly extended to the generalized case. Some of the differences are summarized below. Since all matrices we will be considering for our applications are symmetric, we will restrict the discussion to symmetric generalized eigenvalue problems. Let A, B ∈ R n×n be symmetric matrices. We look for pairs (u, λ) ∈ C n × C such that Au = λBu. In analogy to standard eigenproblems, the generalized eigenvalues are the roots of the characteristic polynomial p(λ) = det(A − λB). Clearly, matrix pencils formed by matrices A and B whose null space intersects nontrivially are necessarily singular. It can be achieved very easily, for example by taking Moreover, contrary to standard eigenvalue problems, symmetry of the coefficient matrices is not sufficient to guarantee real eigenvalues. This point must be emphasized to avoid common misconceptions. A sufficient condition for real finite eigenvalues, which will cover all our applications, is given in the following lemma [18]. The proof is included for completeness.  We conclude this section by introducing a few notations. Λ(A, B) will denote the spectrum (i.e. the set of all eigenvalues) of the matrix pencil (A, B) whereas λ(A, B) will refer to individual eigenvalues. For standard eigenvalue problems (when B = I), we will simply write λ(A) instead of λ(A, I). The matrices will be explicitly specified when referring to individual eigenvalues or to the entire spectrum if there is potential confusion. The eigenvalues will always be numbered in ascending algebraic order such that   The next theorem will play a central role in our analysis. Its proof is a direct consequence of the celebrated Courant-Fischer theorem. Theorem 4.3. Let A, B, C ∈ R n×n be symmetric positive definite matrices and let all eigenvalues be numbered in ascending algebraic order. Then Proof. Let u i , v i and w i denote the eigenvectors associated to the eigenvalues λ i (A, B), λ i (A, C) and λ i (C, B), respectively. We define the subspaces Since all matrices are assumed symmetric positive definite, the generalized eigenvectors of all matrix pencils are linearly independent (see Lemma 3.2).
We already know that the minimum in the second expression is attained for S = U k and the maximum in the third expression is attained for S = U k . Moreover, we note that for any given subsets X k ⊆ C n and X k ⊆ C n with dim(X k ) = k and dim(X k ) = n − k + 1. The strategy now consists in cleverly choosing these subsets. In particular, for the upper bound, we obtain and for the lower bound The last theorem is useful in understanding the relationship between the eigenvalues of (K, M ) and those of (K,M ) following a worst case analysis. ProvidedM is symmetric positive definite, a straightforward application of Theorem 4.3, inequalities (4.1a), with A = K, B =M and C = M leads to This result can also be alternatively deduced starting from Ostrowski's theorem (see for instance [20,Theorem 4.5.9]). An immediate conclusion can be drawn from these bounds: ifM is very good preconditioner for M , then the entire spectrum is well approximated since both the lower and upper bounds are close to 1. They reveal the close connection between the quality of a preconditioner and the approximation of the spectrum. However, not all preconditioners are appropriate for our problem. Indeed, one of our main objectives when designing a preconditioner M is to ensure that λ k (K,M ) ≤ λ k (K, M ) for all k = 1, . . . , n. This requirement is sometimes easily fulfilled by simply understanding the relation between M andM in the Loewner partial order. The next lemma gathers some equivalent sufficient conditions guaranteeing that λ k (K,M ) ≤ λ k (K, M ) for all k = 1, . . . , n. Proof. We first prove the equivalence between the statements: We now show that λ k (A,B) ≤ λ k (A, B) for all k = 1, . . . , n. Using for instance the third characterization combined with Theorem 4.3, inequalities (4.1a), we obtain As we will show, these conditions are actually satisfied for an appropriate definition of mass lumping.
Definition 4.6 (Lumping operator). Let A ∈ R n×n . The lumping operator L : R n×n → R n×n is defined as where d i = n j=1 |a ij | for i = 1, . . . , n. Remark 4.7. Definition 4.6 defines more specifically the row-sum lumping operator. Some readers might be concerned about the absolute values included in our definition. They are absolutely necessary for our results to hold as they guarantee positive definiteness of the lumped mass matrix provided A does not contain a row of zeros. For nonnegative matrices (i.e. matrices with nonnegative entries), our definition obviously coincides with the usual row-sum mass lumping. Some typical examples include isogeometric analysis [4], classical low order finite element methods and nonnegative partition of unity methods applicable to generalized or extended finite element methods [21]. Even for more classical methods, nonnegative mass matrices can still be constructed for special sets of basis functions and finite elements [22].
The next lemma is a first step in analyzing the preconditioned spectrum (M, L(M )).
Lemma 4.8. Let A ∈ R n×n be symmetric positive definite. Then Proof. We prove each point below.
showing that 1 is an upper bound on the spectrum. 2. Note that if A is nonnegative, then (1, e) is an eigenpair of (A, L(A)), where e is the vector of all ones. Thus, the upper bound is always attained.
The next corollary is one of our main results.
Corollary 4.9. Let A, B ∈ R n×n be symmetric positive definite matrices and let all eigenvalues be numbered in ascending algebraic order. Then Proof. The result follows immediately from Lemma 4.8 and Corollary 4.5.
In words, lumping the mass matrix leads to an underestimate of the eigenvalues whereas lumping the stiffness matrix leads to an overestimate of the eigenvalues. Clearly, the above inequality is not of any practical interest for our problem and it is only mentioned for completeness.
In our quest to ensure an underestimate of the eigenvalues, we will heavily rely on the sufficient conditions stated in Corollary 4.5. In fact, the positive semi-definiteness of the error is also the main building block for ad hoc mass scaling strategies [23,24,25]. The objectives of mass lumping and mass scaling are similar, however, mass scaling primarily focuses on improving the stability. The issue of solving linear systems with the scaled mass matrix is usually disregarded and at best only addressed in a later design phase. On the contrary, mass lumping addresses both issues simultaneously.
Nevertheless, although the conditions stated in Corollary 4.5 are sufficient, they are not necessary, as the following example shows.
Example 4.11. Consider the symmetric positive definite matrices The eigenvalues of (A, B), (A,B) and of the error E =B − B are Thus, the error E is indefinite and yet λ k (A,B) ≤ λ k (A, B) for k = 1, 2.
The theoretical results are now illustrated on a numerical example.
Example 4.12 (Preconditioner comparison). We illustrate the theoretical bounds with a first numerical example. We consider an isogeometric discretization of our model problem (2.1) on a stretched square already used by previous authors [26,27]. Quadratic B-splines of maximal smoothness are used with 20 subdivisions in each parametric direction. All our experiments are done with GeoPDEs [28], an open source Matlab/Octave package for isogeometric analysis. The perturbed spectrum Λ(K,M ) is compared for two different choices ofM , namely the lumped mass matrix introduced in Definition 4.6 and a very good preconditioner developed by Loli et al. [27], which is provably optimal under mesh refinement. The interested reader is referred to the original article for its definition and construction. This experiment has two purposes: firstly, it illustrates the close connection between the quality of a preconditioner and the approximation of the eigenvalues. Secondly, it supports our first theoretical finding regarding mass lumping and pinpoints some characteristic features so far not captured by our bounds. The eigenvalue ratios together with their lower and upper bounds from (4.2) are shown in Figure 4.1 for the lumped mass matrix and for the preconditioner of Loli et al. In Figure 4.2, the exact discrete spectrum is compared to the approximate spectrum obtained with our two choices of preconditioners. The excellent preconditioner of Loli et al. leads to a nearly perfect approximation of the entire spectrum while mass lumping leads to an underestimate of the eigenvalues, in agreement with Corollary 4.9.   Example 4.12 also illustrates a typical and widely appreciated feature of mass lumping: In Figure 4.1, the smallest eigenvalues are reasonably well approximated, whereas the larger eigenvalues are strongly underestimated. It explains why mass lumping does not lead to completely erroneous approximate solutions. This feature is explored more thoroughly in a second and simpler example.
Example 4.13 (h-refinement for 1D Laplace). We perform an h-refinement experiment for the isogeometric discretization of the 1D Laplace eigenvalue problem on the unit line with homogeneous Dirichlet boundary conditions. The relative eigenvalue error λi−λi λi for the first 20 discrete eigenvalues is shown in Figures 4.3 for increasingly fine meshes. The spline order is fixed at p = 6. As expected, for a given spline degree the relative error decreases for increasingly fine meshes but is not even across all eigenvalues: the smallest eigenvalues are always much better approximated.  The reasons why mass lumping gives a much better approximation of the smallest eigenvalues can be intuitively explained with the perturbation theory of generalized eigenvalue problems pioneered by Stewart in the early 1970s. For generalized eigenvalue problems, the eigenvalue error is often measured in the chordal metric [18,29,30,31,32] because its ability to treat small and large eigenvalues uniformly. Error bounds are also alternatively formulated in terms of eigenangles, which were first introduced by Stewart in 1979 [30] and are related to eigenvalues by Similarly to eigenvalues, eigenangles satisfy a min-max characterization, which is the basis for deriving error bounds on eigenangles. The classical perturbation theory uniformly bounds the quantity sin(|θ i − θ i |) by the ratio of some error measure over the definiteness of the pencil. Due to the highly nonlinear nature of the cotangent function, large eigenvalues undergo large perturbations while small eigenvalues undergo small perturbations, which is consistent with the results in Figure 4.2. Quantitatively speaking though, the perturbations bounds from [30] seemed sharper for the smallest eigenvalues than for the largest ones in our numerical experiments. The classical perturbation theory may not be satisfactory for our problem for two main reasons. Firstly, the error E is treated as a random symmetric perturbation whereas in our problem the lumped mass matrix is constructed from the consistent mass matrix and therefore the error inherits some structure, which is ignored in the classical perturbation theory. Secondly, all error bounds we are aware of are uniform, meaning the error (in whichever measure is used) is bounded independently of the eigenvalue/eigenangle number similarly to Weyl's theorem for standard eigenvalue problems [20, Equation (6.3.4.1)]. Such bounds cannot explain differences in relative error as we are witnessing in Figure 4.3. Instead, focus should be directed at understanding how the error is interacting with the eigenspace. Surprisingly little attention has been paid to this problem and the few results we are aware of are limited to standard eigenvalue problems [33]. Such error bounds are not only valuable for mass lumping but also for selective mass scaling, where specific eigenvalues are selectively driven down [23]. Perturbation bounds that take into account eigenspaces are established in Theorem 4.14.
denote the eigenpairs of (A, B) and (Ã,B), respectively, where the eigenvectors u i andũ i are assumed B-orthonormal andB-orthonormal, respectively. Then Proof. We prove statement (4.3a) by first noticing that for any given vector x ∈ R n and any given scalar IfB is positive definite, there exists a matrixŨ ∈ R n×n ofB-orthonormal eigenvectors that diagonalizes (Ã,B) such thatŨ TÃŨ =D andŨ TBŨ = I, whereD is the diagonal matrix of eigenvalues of (Ã,B) (Lemma 3.2). Then, by taking the norm, we obtain Noticing that Ũ 2 In particular, the inequality trivially holds for any σ ∈ Λ(Ã,B). After choosing σ = λ i and x = u i , we obtain a residual bound. Since Au i = λ i Bu i , we finally obtain Statement (4.3b) follows by swapping the roles of (A, B) and (Ã,B).
A direct comparison between (4.3a) and (4.4) is generally not possible since they do not bound the same quantities. Nevertheless, if the eigenvalue gap is large enough relative to the error, then the perturbed eigenvalue that attains the minimum in (4.3a) is actually the ith one [33]. In this case, Theorem 4.14 is a clear improvement.
In the context of mass lumping, the error component in the eigenspaces associated to the smallest eigenvalues tends to be quite small. Moreover, for sufficiently fine meshes, the discrete eigenvalues are sufficiently close to the exact eigenvalues, which are well separated for 1D problems. Thus, Theorem 4.14 might very well directly provide an upper bound on the relative error for the smallest eigenvalues. This assertion is tested in the next example.
Example 4.16. The quality of the upper bounds given in (4.3a) and (4.4) is assessed numerically for a first order isogeometric discretization of the 1D Laplace. Figure 4.4 indicates that for the first few eigenvalues,λ i is the perturbed eigenvalue closest to λ i . However, it is no longer the case for the upper 82 % of the spectrum. The error grows faster than the relative gap and the pairing is lost. From Figure 4.4, the reader might misleadingly think that the upper bound in (4.3a) is actually always an upper bound on the relative error between λ i andλ i , regardless of the eigenvalue gap. Unfortunately, other experiments showed the contrary. The results using the upper bound in (4.3b) were quite similar and are omitted.

Error
Relative eigenvalue error

A generalization of mass lumping
For applications in explicit dynamics, the different approximation quality at different branches of the spectrum is an attractive feature of mass lumping. A good approximation of the smallest eigenvalues allows us to recover reasonably good solutions whereas a rather poor and gross underestimate of the largest eigenvalues, which usually do not contribute much to the solution, improves the stability of the time integration scheme. We will now design an entire class of preconditioners preserving those nice properties and generalizing the simple row-sum lumped mass matrix to more complex structures. A small preparatory lemma is needed before describing this class. Thus, ∀u ∈ R n , u = 0 where e is the vector of all ones. We now extend the proof to the more general case. Let us assume that A contains m rows of zeros. Since A is symmetric, if the ith row is zero, so is the ith column. Without loss of generality, we may assume that the rows and columns of A are ordered such that In some sense, if A is nonnegative, Lemma 4.17 states that the "upper bound" A L(A) is "attained" in the Loewner ordering. Namely, the error L(A) − A is not only positive semidefinite but also singular. For the sake of eigenvalue approximation, we are especially interested in singular errors. In fact, positive definite errors will necessarily effect the smallest generalized eigenvalues, which we of course want to avoid. Definition 4.18 (Preconditioners P i ). Let A ∈ R n×n be a symmetric positive definite matrix and consider the matrix splitting A = D i + R i where D i consists of all super and sub-diagonals smaller than i and R i is the remainder. We define the sequence of preconditioners P i = D i + L(R i ) for i = 1, . . . , n. In particular, we observe that P 1 = L(A) and P n = A.
We are now ready to state an important theorem generalizing Lemma 4.8. Proof. We prove all properties below.
1. We first prove that P i is symmetric positive definite for all i = 1, . . . , n. Symmetry is a straightforward consequence of its definition and the symmetry of A. Positive definiteness immediately follows from Lemma 4.17 and the positive definiteness of A: From the equivalent characterizations in Corollary 4.5, Λ(A, P i ) ⊂ (0, 1]. 2. For the second property, we note that where we have denoted ∆D i = D i+1 − D i and used the definition of the matrices R i and R i+1 to deduce that We note that ∆D i is obviously symmetric, but it may be singular. Nevertheless, Lemma 4.17 covers for such a situation and states that ∆D i − L(∆D i ) 0. Since P i and P i+1 are symmetric positive definite and P i P i+1 , Corollary 4.5 concludes that The same argument shows that more generally, P i P j for all j ≥ i. Thus, λ n (P j , P i ) ≤ 1 for all j ≥ i.
3. By the first property, λ n (A, P i ) ≤ 1. If A is nonnegative, this upper bound is attained. Indeed, if A is nonnegative, all preconditioners P i for i = 1, . . . , n are also nonnegative and all matrices have the same rowsum by construction. Thus, Ae = (D i + R i )e = (D i + L(R i ))e = P i e where e is the vector of all ones. Therefore, (1, e) is an eigenpair of (A, P i ) and (P j , P i ) for all i, j = 1, . . . , n. Consequently, λ n (A, P i ) = 1 and λ n (P j , P i ) = 1 for all j ≥ i.

Lemma 4.17 states that L(A)
A. Increasing the bandwidth of the preconditioner allows us to find increasingly better preconditioners that are "closer" to A in the Loewner ordering: L(A) = P 1 P 2 · · · P n−1 P n = A.
Theorem 4.19 then simply states that Λ(A, P i ) converges monotonically to 1 for increasing values of i. In particular, the theorem can be used to construct a sequence of preconditioners for the mass matrix and the first element of this sequence coincides with the usual row-sum lumped mass matrix. Since Λ(M, P i ) converges monotonically to 1 for increasing values of i, Λ(K, P i ) is expected to yield an increasingly good approximation of Λ(K, M ). Indeed, since P i P i+1 , Corollary 4.5 yields λ k (K, P i ) ≤ λ k (K, P i+1 ). By repeating the same argument for all i = 1, . . . , n − 1, we obtain In other words, Λ(K, P i ) converges monotonically from below to Λ(K, M ) for increasing values of i. In practice, we are only interested in small values of i such that the preconditioners P i have narrow bandwidths. Linear systems of size n with banded symmetric positive definite matrices of bandwidth b n can be solved in O(nb 2 ) floating point operations (flops), which represents a considerable saving over the O(n 3 ) complexity for general matrices using standard Gaussian elimination [35,36]. The tridiagonal case is particularly appealing since linear systems can be solved in O(n) flops using the Thomas algorithm. In general, choosing the right preconditioner is a compromise between accuracy and computational cost.
Example 4.20 (Lumped mass preconditioners P i -1D). We consider an isogeometric discretization of the 1D Laplace equation on Ω = (0, 1) with B-splines of order p = 3 and p = 5 with C p−1 continuity and 400 subdivisions. For a pth order discretization with C p−1 continuity, the mass and stiffness matrices are banded matrices of bandwidth p. Therefore, they have 2p + 1 bands. We construct the first three lumped mass preconditioners of the sequence (P 1 , P 2 and P 3 ). The preconditioner P i is a banded matrix of bandwidth i − 1 and has 2i − 1 bands. Thus, P 1 , P 2 and P 3 are diagonal, tridiagonal and pentadiagonal matrices, respectively. None of these preconditioners are able to represent the mass matrices exactly. Already for p = 3, the mass matrix has 7 bands. However, the further the bands are from the diagonal, the smaller the magnitude of their entries. The spectrum of (K, M ) and (K, P i ) for i = 1, 2, 3 are compared in Figures 4.5 and 4.6 for p = 3 and p = 5, respectively. The approximation quality improves tremendously already for small values of i. For p = 3, the preconditioner P 3 leads to a nearly perfect approximation. We notice that Λ(K, P i ) converges monotonically from below to Λ(K, M ), in perfect agreement with our theoretical results. Alert readers will also notice the "spikes" appearing towards the end of the spectrum. The largest eigenvalues of (K, M ) grossly overestimate those of the continuous operator and for this reason were coined outliers in [6] where they were first discussed. Unfortunately, in our case, Λ(K, P i ) approximates the entire spectrum of (K, M ), including the outliers. Fixing this issue will be the object of future work.   The performance of these preconditioners in approximating the spectrum of (K, M ) is directly related to the clustering on the eigenvalues of (M, P i ) around 1. The eigenvalues of the matrix pencils (M, P i ) for i = 1, 2, 3 are shown in Figures 4.7 and 4.8 for p = 3 and p = 5, respectively. As predicted theoretically, the spectrum of the preconditioned matrices converges monotonically to 1.    can also be used to construct a sequence of preconditioners for the stiffness matrix K. It can be easily shown that in this case the approximate spectrum converges monotonically from above to Λ(K, M ) for increasing values of i. However, the main computational bottleneck in explicit dynamics consists in solving linear systems with the mass matrix, not the stiffness matrix. Therefore, we will not explore this option any further.
Theorem 4.19 is based on purely algebraic considerations on the entries of the matrices. Thus, it holds regardless of how complicated the physical domain is or how large the dimension of the problem becomes. However, for multidimensional problems the mass matrix no longer has a banded structure. Therefore, if banded preconditioners are constructed, the preconditioned spectrum will most likely converge very slowly to 1. This foreboding is confirmed in the next two-dimensional example.      It is well-known in preconditioning practice that good preconditioners must preserve the essential structure of the matrices [37,38]. Therefore, it would be wise to further extend the family of lumped mass preconditioners to better account for the structure. For tensor product basis functions, there exists a natural generalization to multidimensional problems.

Multidimensional problems
We first discuss the case of a scalar PDE in 2D. The extension to 3D is straightforward and vector problems will be discussed at the end of this section. If the physical domain is a square, the density function is separable and tensor product basis functions are used, the mass matrix is expressed as M = M 1 ⊗ M 2 where M 1 and M 2 are mass matrices for univariate problems (see for instance [26]). We then consider the family of lumped mass preconditioners given by P ij = P 1,i ⊗ P 2,j where P 1,i and P 2,j are banded preconditioners defined in Definition 4.18 for some indices 1 ≤ i, j ≤ n. This definition is the natural extension to two-dimensional problems, as shown in the next theorem. Definition 4.23 (Preconditioners P ij ). Let A = A 1 ⊗ A 2 where A 1 , A 2 ∈ R n×n are symmetric positive definite matrices. We define the family of preconditioners P ij = P 1,i ⊗ P 2,j , where P 1,i and P 2,j are preconditioners constructed from A 1 and A 2 , respectively, for some indices 1 ≤ i, j ≤ n according to Definition 4.18. In particular, we observe that P nn = A.
Proof. We prove all properties below using the results of Theorem 4.19. Let us first define the product of two sets S 1 , S 2 ⊆ C as S 1 S 2 = {s 1 s 2 : s 1 ∈ S 1 , s 2 ∈ S 2 }.
1. We will first prove that A and P ij are symmetric positive definite. Our arguments are based on some basic properties of the Kronecker product. The reader may for instance refer to [39,Chapter 4] for their proof. For By the first and second property, since A 1 , A 2 and P 1,i , P 2,j are symmetric positive definite for all i, j = 1, . . . , n, then so are A 1 ⊗ A 2 and P 1,i ⊗ P 2,j . Therefore, λ k (A, P ij ) > 0 for all k = 1, . . . , n. Then, using Lemma 4.1 and the properties of the Kronecker product Thus, using the first property of Theorem 4.19, λ n (A, P ij ) = λ n (A 1 , P 1,i )λ n (A 2 , P 2,j ) ≤ 1 and Λ(A, P ij ) ⊂ (0, 1]. 2. Using Theorem 4.3, inequality (4.1a), we obtain λ k (A, P sq ) ≤ λ k (A, P ij )λ n (P ij , P sq ) 1 ≤ k ≤ n 2 .
Using again Lemma 4.1 and the properties of the Kronecker product, λ n (P ij , P sq ) = λ n (P i , P s )λ n (P j , P q ). Finally, when proving the third property of Theorem 4.19, we showed that λ n (P j , P i ) ≤ 1 for all j ≥ i. Therefore, λ n (P ij , P sq ) ≤ 1 if i ≥ s and j ≥ q and the result follows. 3. Moreover, if A 1 and A 2 are nonnegative, so are the preconditioners P 1,i and P 2,j for all i, j = 1, . . . , n constructed from A 1 and A 2 and Theorem 4.19 states that λ n (A 1 , P 1,i ) = 1 and λ n (A 2 , P 2,j ) = 1 for all i, j = 1, . . . , n. Therefore, λ n (A, P ij ) = λ n (A 1 , P 1,i )λ n (A 2 , P 2,j ) = 1.
Remark 4.25. The previous theorem is extended straightforwardly to matrices A = A 1 ⊗A 2 ⊗A 3 and the associated family of preconditioners P ijk = P 1,i ⊗ P 2,j ⊗ P 3,k where P 1,i , P 2,j and P 3,k are constructed from A 1 , A 2 and A 3 , respectively, for some indices 1 ≤ i, j, k ≤ n according to Definition 4.18.
Linear systems with a Kronecker product structure can be solved very efficiently and especially if the factor matrices themselves have good properties. It is the main motivation for designing Kronecker product preconditioners [26,27,40]. In the 2D case, if all factor matrices have size n and c 1 and c 2 denote the cost of solving a single linear system with P 1,i and P 2,j , respectively, then the cost of solving a linear system with P ij is O(n(c 1 + c 2 )). Thus, the preconditioner of largest bandwidth dictates the overall cost. If c 1 ≈ c 2 ≈ n, this cost amounts to O(n 2 ), which is linear in the size of the system. The analysis can be easily generalized to the 3D case.   Lumped mass preconditioners -2D). We consider exactly the same problem as in Example 4.22: the 2D Laplace model problem on the unit square Ω = (0, 1) 2 discretized with maximally smooth B-splines of order p = 3 and p = 5 and 20 subdivisions in each parametric direction. For this simple model problem, the mass matrix is expressed as a Kronecker product M = M 1 ⊗ M 2 , with M 1 = M 2 because of the same discretization parameters in each direction. We now compare the preconditioners P i for i = 1, 2, 3 from the previous example to the preconditioners P ii constructed following Definition 4.23. For i = 1, the two preconditioning strategies coincide. Hence, Figures 4.13 and 4.14 show the comparison for i = 2 and i = 3, respectively and p = 3. The same comparison for p = 5 is shown in Figures 4.15 and 4.16. In both cases, the preconditioners P ii lead to a significant improvement. However, it is not a fair comparison since they also usually have many more nonzero entries. Thus, let us select two different preconditioners with exactly the same sparsity pattern and the same number of nonzeros. For instance, P 2 and P 12 are both tridiagonal. Yet, Figure 4.17 reveals that the preconditioner P 12 performs much better.         Example 4.28 (Accuracy test -1D and 2D). As we have discussed in Section 1, the row-sum lumped mass matrix limits the convergence rate of the smallest eigenvalues to quadratic order (at least for 1D problems and for quadratic and cubic B-spline discretizations [6]). This property is now investigated numerically for our new class of lumped mass preconditioners for both 1D and 2D model problems. We consider a cubic isogeometric discretization of the 1D and 2D Laplace on the unit line and unit square, respectively, with mixed boundary conditions. The smallest eigenfrequency ω 1 = √ λ 1 for these model problems is given by We report the relative error in Figures 4.18 and 4.19 for our 1D and 2D problems, respectively. Unfortunately, our class of lumped mass preconditioners does not seem to improve the convergence rate but may improve the constant by several orders of magnitude. This fact is interesting given how close P 3 (or P 33 ) is to the consistent mass matrix M . Small differences are enough to ruin the convergence rate, which we expect is due to the consistency error of the method. The same trend was also observed for other types of boundary conditions.   If M S denotes the mass matrix of a multidimensional scalar PDE, and if the same approximation spaces are used for each component of the approximate solution of the vector PDE, then the mass matrix of a vector problem M V is either I d ⊗ M S or M S ⊗ I d depending on the ordering of the basis functions [41]. With the first definition M V is a block diagonal matrix. This ordering is used in GeoPDEs for instance [28]. In order to preserve the structure, the preconditioner is naturally defined as P V = I d ⊗ P S or P V = P S ⊗ I d . Again, due to the Kronecker product structure, all the properties of P S carry over to P V . Numerical testing for elasticity problems produced figures similar to those reported in Example 4.27 and are omitted.

Nearest Kronecker product approximation
As we have seen in Section 4.3, the performance of the lumped mass preconditioners can be drastically improved for tensor product basis functions if the mass matrix is expressed as a Kronecker product. Unfortunately, this structure cannot be guaranteed for most problems of practical interest. In several cases, the entries of the mass matrix are computed from The integration is done in the parametric domainΩ whereB i (x) are tensor product basis functions defined in the parametric domain (for instance B-splines in isogeometric analysis),ρ(x) = ρ(F (x)) is the density function, F :Ω → Ω is the mapping from parametric to physical domain and J F (x) is its Jacobian matrix. The mass matrix looses its Kronecker product structure when the integrand is no longer a separable function. A nontrivial density function, geometry or the use of NURBS basis functions are some of the causes. Nevertheless, for tensor product basis functions, the mass (and stiffness) matrix can be generally very well approximated by a sum of Kronecker products, a fact which has already been exploited for designing fast assemblage algorithms [26,42,43,44,45]. For 2D problems, we are interested in finding the best factor matrices B and C such that M ≈M = B ⊗ C. Several different approaches have been proposed in the literature. The one proposed in [42] and later simplified in [43], is tailored to PDE problems and attempts to directly approximate all non-separable functions in the integral (5.1) by a sum of separable functions. In other words, it finds a low-rank approximation for the integrand itself.
Another approach, used in [44] and based on earlier work from [46,47], is more general and does not require any knowledge of the origin of the matrix. Of course, the first approach seems better suited to our problem. However, its implementation, although simplified, remains technical. The implementation of the second approach, in comparison, is rather straightforward. The method relies entirely on standard matrix operations. For this reason, we will assess if the easiest method can deliver useful results. We will initially restrict the discussion to 2D problems by closely following the presentation in [46]. The extension to 3D problems will be discussed subsequently.

Approximation for 2D problems
Given a matrix M ∈ R n×n with n = n 1 n 2 , our goal is to find the best factor matrices B ∈ R n1×n1 and C ∈ R n2×n2 such that φ M (B, C) = M − B ⊗ C F is minimized. This problem was first investigated by Van Loan and Pitsianis [46]. They observed that since both the Kronecker product B ⊗ C and vec(B) vec(C) T form all the products b ij c kl for i, j = 1, . . . , n 1 and k, l = 1, . . . , n 2 but at different locations, there must clearly exist a linear mapping R : R n1n2×n1n2 → R n 2 1 ×n 2 2 such that R(B ⊗ C) = vec(B) vec(C) T . This mapping can be defined explicitly. Consider the following block matrix where A ij ∈ R n2×n2 for i, j = 1, . . . , n 1 and define the mapping as Then, by construction, for a matrix A = B ⊗ C, More generally, since the mapping R inherits the linearity of the vectorization operator, it transforms a sum of r Kronecker products into a rank r matrix. A matrix is said to have Kronecker rank r if it can be expressed as the sum of r Kronecker products. Thus, A has Kronecker rank r if and only if R(A) has (matrix) rank r. The mapping R was originally called rearrangement in [46] but later authors have also called it reordering. In any case, we must stress that R(A) is generally not a permutation of the rows and columns of A or any change of basis. The minimization of M − B ⊗ C F in the Frobenius norm allows to rewrite the problem in a more familiar way. The Frobenius norm is defined as the square root of the sum of squares of the matrix entries and does not depend on the location of these entries. Therefore, the mapping R can be applied to a matrix without affecting its Frobenius norm and the minimization problem becomes Thus, finding the best factor matrices B and C is equivalent to finding the best rank-1 approximation of R(M ). The transformation from M to R(M ) requires knowing the block-partitioning of the matrix M , which in the context of isogeometric analysis is always clear from the underlying discretization. The next theorem summarizes a very useful result in this context. Thus, the properties of M are inherited by its approximation B ⊗C. Once the Kronecker product approximation has been computed, the lumped mass preconditioners can be constructed from the factor matrices B and C following Theorem 4.24. The overall construction consists in a two-level approximation of the mass matrix. Algorithm 5.1 is a practical algorithm for computing the preconditioner. For simplicity, we will assume that the preconditioners of the factor matrices have equal bandwidth. More specifically, a practical implementation in Line 2 would extract a dense submatrix from R(M ) containing all its nonzero entries and then compute a low-rank approximation of that submatrix only. The truncated singular value decomposition is an obvious choice for computing a low-rank approximation. However, there exists a host of other black-box algorithms for computing good but not optimal low-rank approximations. Some of them do not require the matrix to be explicitly available. It was the main motivation in [44] for using Adaptive Cross Approximation (ACA). In our context, it means we could potentially compute a good Kronecker product approximation of the mass matrix without ever assembling it. In Line 5, the preconditionerP ii is not necessarily formed explicitly. Only the factorsP 1,i andP 2,i are stored.
Clearly, we can expect the approximation error of M by B ⊗ C to be related to the low-rank approximation error of R(M ) by vec(B) vec(C) T . The next argument attempts to formalize this intuition. Assuming R(M ) has rank r, its reduced singular value decomposition (SVD) is Our ability to approximate the spectrum of M also depends on the low-rank approximation error of R(M ). This fact is a direct consequence of [20, Corollary 6.3.8], applied to the symmetric matrices M andM The error we are committing will surely affect all eigenvalues. However, our main concern when approximating M byM is to guarantee that λ k (K,M ) remains close to λ k (K, M ) for the smallest eigenvalues while underestimating the largest eigenvalues. We have already noted in Corollary 4.5 that if the error E =M − M is symmetric positive semidefinite, then λ k (K,M ) ≤ λ k (K, M ). Unfortunately, this property is not satisfied by the nearest Kronecker product preconditioner, the error matrix E being in general indefinite, and the eigenvalues of (K,M ) might be greater than those of (K, M ). However, as we have seen in Example 4.11, the positive semidefiniteness of the error only guarantees an underestimate of the eigenvalues, and it is not a necessary condition. Moreover, our strategy relies on a two-level approximation of M . It is first approximated by a Kronecker structured matrixM and theñ M is approximated by a lumped mass preconditionerP ii . The relation between λ k (K,M ) and λ k (K,P ii ) is well understood thanks to Theorem 4.24, guaranteeing that λ k (K,P ii ) ≤ λ k (K,M ). The problem of understanding the relation between λ k (K, M ) and λ k (K,M ) once again amounts to assessing the quality of the preconditionerM . From (4.2), Although λ n (M,M ) may be greater than 1 due to the indefiniteness of the error, we will show numerically that it remains very close to 1. The use of Kronecker product approximations as preconditioners is not a new idea. It has already been explored for various applications including image processing [48], Markov chains [49], stochastic Galerkin finite element discretizations [50] and isogeometric analysis [26,27,40]. However, only few authors have undertaken the challenge of obtaining bounds for the preconditioned spectrum and most attempts are application specific. The following theorem is a variation of a general result presented in [50,Lemma 5.3] obtained by exhibiting the singular values and giving a more explicit bound on the condition number of the preconditioned matrix.
where κ denotes the spectral condition number.
Proof. By reverting the ordering, M is expressed as a sum of r Kronecker products where u i = vec(U i ) and v i = vec(V i ) for i = 1, . . . , r. Similarly, by construction, R(M ) = σ 1 u 1 v T 1 and thus M = σ 1 (U 1 ⊗ V 1 ). The symmetry and positive definiteness of U 1 and V 1 has already been established in [46,Theorem 5.8]. We now establish symmetry of all U i and V i for i = 1, . . . , r using proof arguments similar to those presented in [49,Theorem 4.1]. By the singular value decomposition, R(M )v k = σ k u k and R(M ) T u k = σ k v k for k = 1, . . . , r. Since By reshaping u k and v k to matrices, we obtain Since all B i and C i are symmetric by assumption, so are U k and V k for all k = 1, . . . , r. We now establish the bound on the condition number by noting that Assuming that δ < 1, sinceM −1/2 MM −1/2 is symmetric positive definite, The previous theorem is not entirely satisfactory since δ still depends on the unknown quantities max k |λ k (U i , U 1 )| and max k |λ k (V i , V 1 )|, which we were unable to bound in any convenient way. Nevertheless, our numerical experiments indicated that these quantities were O(1). Most importantly, δ depends on the singular value ratios σi σ1 for i = 2, . . . , r. Therefore, a large singular value gap (σ 2 σ 1 ) is important for constructing a good preconditioner and the smallest singular values will not contribute much to the sum. The quality of the approximation is now assessed numerically on a few examples. We consider the same discretization parameters in all examples, namely cubic B-splines of maximal smoothness and 20 subdivisions in each parametric direction. The factor matrices B and C are obtained by truncating the singular value decomposition of a dense submatrix extracted from R(M ) and containing all its nonzero entries. The rank of R(M ) is exactly 2 for Example 1 since | det(J F (x)| is the sum of two separable functions for the stretched square domain. The second example is more challenging with R(M ) having numerical rank 22. For each case, the discrete eigenvalues using the consistent mass matrix are compared to the approximate eigenvalues using the nearest Kronecker product (NKP) preconditioner with and without mass lumping. The curve usingP 33 nearly overlaps with the one for the NKP preconditioner and is omitted.   The condition number of the preconditioned mass matrix with the NKP preconditioner is 2.18 and 39.1 for Examples 1 and 2, respectively. The assumptions of Theorem 5.2 are satisfied for the first example and the upper bound on the condition number is 2.73. The assumptions are not satisfied for the second example due to the much slower decay of the singular values. In Figure 5.2, we notice that the largest eigenvalues of (K,M ) are actually slightly greater than those of (K, M ). This situation cannot be prevented but our second level approximation involving mass lumping usually fixes it. On the contrary, in Figure 5.4, the largest eigenvalues of (K, M ) are greatly underestimated when using the NKP preconditioner. This situation is very favorable and was generally experienced in our experiments with complicated problems. Any mass lumping done subsequently will bring down the eigenvalues further. However, the bandwidth must still be large enough for the smallest eigenvalues to be well approximated. The selective damping of the higher modes will be addressed in future work.

Approximation for 3D problems
The natural analog of the minimization problem for 3D discretizations is min φ M (B, C, D) = min M − B ⊗ C ⊗ D F where M ∈ R n×n with n = n 1 n 2 n 3 and the factor matrices B ∈ R n1×n1 , C ∈ R n2×n2 and D ∈ R n3×n3 . This extension was not originally investigated by Van Loan and Pitsianis [46] but was later analyzed in [49]. The definition of the mapping R must be adjusted such that R(B ⊗ C ⊗ D) = vec(B) • vec(C) • vec(D). More generally, it transforms a Kronecker product of p factor matrices to a pth order rank-1 tensor obtained by taking the outer product of the vectorization of the factor matrices. The symbol • is used to denote the outer product. Note that this definition is consistent with the one previously introduced for the 2D case since b • c = bc T . The minimization problem then translates to finding the best rank-1 tensor approximation of R(M ). Unfortunately, contrary to the matrix case, we cannot hope to compute this best rank-1 approximation and we must content ourselves with good rank-1 approximations such as those computed with the high order SVD (HOSVD) [51]. Another major shortcoming of the extension to higher dimensions is the lack of theoretical results related to the computed factor matrices. Some proof arguments used in [46] for the two dimensional case cannot be straightforwardly extended to higher dimensions. On the other hand, Algorithm 5.1 can be straightforwardly adapted to higher dimensional problems. The computational results were quite similar to those reported in Example 5.3 and are omitted.

Conclusion
Mass lumping is used on a daily basis for industrial applications in structural dynamics and consists in replacing the consistent mass matrix with some kind of easily invertible (typically diagonal) approximation. In this article, we have unraveled some attractive mathematical features of mass lumping partly responsible for its success. Our results hold under rather broad assumptions, indicating that some properties are not specific to any discretization method. More generally, we have highlighted the close connection between the quality of a preconditioner and the approximation of the eigenvalues. Based on this insight, we have generalized the concept of mass lumping and its properties to matrices with more complex structures including Kronecker products, whose structure can be used to solve linear systems very efficiently. We have then proposed a two-level approximation of the mass matrix in case it is not already expressed as a Kronecker product. Numerical experiments have highlighted the practical usefulness of this approximation in the context of isogeometric analysis by increasing the critical time step of explicit time integration schemes. However, our generalization of mass lumping improves the approximation of all discrete eigenvalues, including the outliers. Outlier eigenvalues are a distinctive feature of maximally smooth isogeometric discretizations. Apart from being completely meaningless, they significantly constrain the critical time step of explicit time integration schemes. Understanding and removing these outliers is therefore of high interest and several follow-up papers have been written on this topic; see e.g., [17,45,52,53,54]. Future work will attempt to selectively approximate the accurate low-frequency part of the spectrum while undermining the outliers. One possible research direction consists in combining the mass lumping techniques presented here with the optimal spline spaces in [55] which were shown to be outlier-free in [53] and [54].