Tensor Train Accelerated Solvers for Nonsmooth Rigid Body Dynamics

The discrete element method (DEM) [1] is one of the most widely used approaches to simulate the multibody dynamics such as in granular materials. This method considers the granular medium as a collection of discrete particles; each responding to body forces such as gravity, inertia or drag, as well as repulsive or dissipative forces caused by contact. Materials in granular form are omnipresent in industry and understanding their dynamics is crucial for a broad range of application fields, including terramechanics, active media, additive manufacturing, nanoparticle self-assembly, avalanche dynamics, composite materials, pyroclastic flows, etc.

Distinct instances of the DEM are defined essentially by their modeling of contact forces, and thus in how collisions are handled. We distinguish two main types of discrete element methods in the literature: force penalty methods (DEM-P) and complementarity formulations (DEM-C). Penalty methods introduce one or multiple layers of spring-like forces between objects in contact, and may introduce additional fields to represent friction. These methods are widely used due to being computationally inexpensive and easy to implement. However, in many cases they introduce high stiffness, requiring extremely small time steps to retain stability during collisions. An in-depth and informative comparison between the two can be found in Ref. [2].

DEM-C methods such as in Refs. [3] and [4] enforce contacts using complementarity constraints, leading to a differential variational inequality problem upon discretization. This allows for the use of larger time steps in its integration, as contacts are enforced geometrically. Given a time-stepping scheme for this problem, a nonlinear complementarity problem (NCP) must be solved to compute the corresponding contact forces at each time-step. One common way to solve this optimization problem is via a linearization of the constraints, producing a linear complementarity problem [5,6]. An alternative relaxation method produces an equivalent, convex cone complementarity problem (CCP). Efficient quadratic cone programming techniques such as those in Refs. [7–11], and [4] have been proposed to solve the CCP.

In Ref. [12], the authors performed a comparison of these quadratic programming techniques focused on determining which class of methods performed best: second-order optimization methods, i.e., those that use Hessian information (Interior Point), or first-order methods, i.e., methods using only gradient information (Jacobi, Gauss-Seidel, Projected Gradient Descent). Second-order methods display quadratic convergence in a neighborhood of the solution, and thus their convergence is much faster and more problem-independent than that of linearly convergent first-order methods, requiring at least 1 to 2 orders of magnitude less iterations to reach a desired accuracy across experiments in this work. However, each iteration in second-order methods requires the solution of a (generally sparse) linear system, which can become costly as the dimensions of the many-body problem increase. For large problems, the added computational effort ultimately eclipses the gains obtained from the reduction in iteration counts.

The preferred solution method for large sparse systems is often a Krylov subspace iterative method. The iteration counts, and thus the performance of these methods are known to be directly affected by the eigenspectrum of the associated matrices. Preconditioning techniques can be used to cluster the eigenvalues away from zero and drastically reduce iteration counts; we refer the reader to Ref. [13] for a general review. Generic sparse preconditioners are most typically based on incomplete or sparsified factorizations, such as the well-known Incomplete LU and Cholesky methods [14]. In Ref. [12], a fast, parallel SaP (split and paralellize) [15] preconditioner was used to accelerate the primal dual interior point (PDIP) method, garnering reductions in iteration counts and execution times.

Most general-purpose preconditioners suffer a trade-off between precomputation costs and the resulting reduction in iteration counts, and their performance is often problem-dependent. Moreover, in the context of optimization problems such as the ones we are interested in, system matrices change every iteration, and it is often impossible or expensive to update the associated preconditioners. Finally, improving their scaling with the number of degrees-of-freedom and the level of sparsity is extremely important, as it is central to remaining competitive in large-scale problems.

In the last decades, a continued effort has been made to produce direct solvers for structured linear systems arising from differential and integral equations. These solvers entirely side-step the challenges related to convergence speed of iterative solvers. They can also lead to dramatic improvements in speed, in particular in situations where a large number of linear systems with coefficient matrices that stay fixed or can be updated via low rank modifications. Additionally, they also provide a methodology to produce robust preconditioners: low accuracy direct solvers may be used in conjunction with iterative refinement or the Krylov subspace method of choice; the use of an approximate low accuracy inverse requires less memory and precomputation time than a direct solver, at the expense of a slight increase in the number of iterations. These features have led to the adoption of these solvers in other areas of scientific computing and statistics.

In Ref. [16], an effective and memory-efficient solver based on the quantized tensor train decomposition (QTT) was presented. By recasting system matrices as tensors, the tensor train (TT) compression and inversion routines were used to produce direct solvers and robust preconditioners for integral equations in complex geometries in three dimensions. Key properties of this solver that differentiate it from other hierarchical matrix approaches feature sublinear computational costs and memory requirements with problem size N, as well as techniques to produce economic updates for a matrix and its inverse across time-steps, even when matrix size changes.

The main goal of this work is to employ a QTT-based approach to provide a radical speed up to second-order optimization methods. We demonstrate its application to interior point methods such as the PDIP in the context of many-body dynamics problems; however, we expect our discussion to apply with little to no modification to a general class of Newton and quasi-Newton type methods. We first study compressibilty of the system matrices and their inverses in TT format for a range of target accuracies. We then demonstrate how factorization re-use can provide significant speed-ups to precomputation costs, reducing costs by orders of magnitude. For three validation tests of common soil mechanics phenomena—sedimentation, blade drafting and direct shear experiments—we show that a TT-based preconditioner displays efficient and robust performance for problems with > 10⁴ bodies, garnering up to an order of magnitude speed-up and greatly improved iteration counts when compared with state-of-the-art ILU-preconditioned methods. We also confirm an extremely significant gain in scaling of precomputation costs: precomputation for the TT preconditioner is sublinear (for all practical purposes, constant) as the number of collisions and matrix size increase.

The QTT decomposition is one of several approaches for approximate solution of linear systems based in hierarchical compression. In the context of sparse structured systems and related factorizations, work has been done for a number of hierarchical matrix formats: hierarchically semi separable [17–23], $H$matrices: [24–26], and fast multipole methods [27,28]. Producing efficient, global factorization updates and dealing with high storage costs is an ongoing challenge in these alternate formats.

In order to ensure clarity of exposition, we include list of notation and terminology employed throughout this work. The tensor train acceleration proposed involves the manipulation of vectors, matrices and tensors of different dimensions, as well as the associated indexing and operators (Table 1).

We include a glossary and notation for important quantities and parameters used throughout the text in Table 2.

We consider a granular material comprised of M rigid particles B_i in $ℝ3$⁠; the position of each body is uniquely described by the coordinates $xi$ for its center of mass and the rotation $Qi$ of a reference frame fixed to the body, represented by a unimodular quaternion. Quaternion representation of rotations in three dimensions is usually favored due to being compact and numerically stable; since we require the four component quaternions to be length one, it introduces only three independent degrees-of-freedom per particle. Let $qi=(xi,Qi)$ then be the six generalized coordinates for the ith body, and $q=(q1,q2,…,qM)∈ℝ6M$⁠.

where $L(q)$ is a linear operator that relates velocities to the rate of change in generalized coordinates, $M$ is the mass matrix and $v,fB,fC∈ℝ6M$ contain each body's translational and rotational velocities, the total body forces and torques applied to them and the reaction forces and torques due to contact dynamics, respectively.

We impose the following contact constraints: no two bodies should penetrate, and if there is contact, a normal force and a tangential frictional force act at the interface. Consider two bodies $Bi1$ and $Bi2$⁠, and let $Φi(q)$ be an unsigned distance function (also known as a gap function) for the pair $i=(i1,i2)$ satisfying $Φi(q)>0$ if the two bodies are separated, $Φi(q)=0$ if they are touching, $Φi(q)<0$ otherwise.

where μ is the static friction coefficient. The feasible set of forces $f=fN+fT$ for the minimization problem in Eq. (4) are known as a friction cone$ϒi={(γi,n,γi,1,γi,2) | ||(γi,1,γi,2)||≤μγi,n}$⁠. The complementarity condition in Eq. (3) is typically abbreviated using the notation $0≤γi,n⊥Φi(q)≥0$⁠. We note that an alternate complementarity formulation for Eq. (3) allowing for partially elastic collisions can be obtained by replacing $Φi$ with the normal velocity as discussed in Refs. [5] and [29].

For a system with $Nc=|A(q,δ)|$ contacts, this adds a number of constraints to the equations of motion proportional to N_c. We note that for dense granular flows, N_c itself scales as O(M) for M bodies.

We note that Eq. (12) is obtained from Eq. (8) via a linearization, dividing by Δt (which does not affect complementarity, but is numerically desirable). This makes the future velocity $vk+1$ the sole variable needed to enforce the complementarity condition. We also note that in this discretization, $(γi,n,γi,1,γi,2)$ constitute contact impulses, i.e., force magnitudes multiplied by the step length Δt.

While this relaxed problem may introduce artifacts when step size Δt, sliding velocity or friction are large, it enables the use of a wide range of quadratic programming methods for its solution:

• Projected Jacobi and Gauss–Seidel methods [4].

• Projected gradient descent methods like accelerated projected gradient descent [10], Barzilai and Borwein [30] and the Kucera and Preconditioned spectral projected gradient with fallback methods in Ref. [8].

• Krylov subspace methods: Gradient projected minimum residual in Ref. [8].

• PDIP [7,11,31–34].

• Symmetric cone interior point [9].

Projected Jacobi and Gauss–Seidel methods, while requiring only fast matrix applies of $N$⁠, have slow, linear convergence, often requiring thousands of iterations to obtain a significant reduction for the objective function. The projected gradient descent and Krylov subspace methods have since been proposed, providing considerable iteration count reductions while retaining cost-efficiency per time-step. However, they remain linearly convergent, and they require an increasing number of iterations as problem size increases for a variety of problems of interest.

Interior point methods are often based on a modified Newton or quasi-Newton step, displaying quadratic convergence near minima and iteration counts which are less dependent on problem size. However, they require the approximate solution of large linear systems in order to produce the Newton step, thus losing this competitive advantage due to computational cost.

The parameter τ controls the strength of the barrier, and as $τ→∞, Bτ(z)→I(z)$ with I(z) the indicator function over (–∞, 0], and the problem in Eq. (21) becomes equivalent to the original constrained program defined by Eqs. (19) and (20). Given the set of problems defined by Eq. (21), interior point methods proceed by following along the corresponding central path${γ*(τ):τ>0}$ of optima, which are located inside of the feasible set, and converge to the solution of the original problem as $τ→∞$⁠.

We note that, as $τ→∞$⁠, Eq. (26) becomes a strict complementarity condition.

where $M̂$ is a diagonal matrix defined by $M̂=∑i=12Ncλi∇2fi(γ)=diag(m̂)$⁠, with $m̂=[μ12λ1,λ1,λ1,…,μ12λNc,λNc,λNc]$⁠. $B,C$ are banded rectangular matrices, and $E$ is diagonal. A typical sparsity pattern for a multibody dynamics problem is shown in Fig. 1.

We note that working with this Schur complement matrix $S$ is widely preferred, as it is symmetric positive definite and smaller in size. We additionally observe that its sparsity pattern across all PDIP iterations is always that of matrix $N$⁠. That is due to the fact that $M̂$ is diagonal and $BE−1C$ is block-diagonal (3 × 3 blocks). Finally, in all of our experiments, $S$ was relatively more compressible in the Tensor train format. For these reasons, our discussion and acceleration efforts below center around solving the associated Schur complement system.

Most state-of-the-art methods for the direct solution and preconditioning of this system, the Tensor train decomposition included, will attempt to exploit sparsity structure in $S$⁠. We recall that $N$ is the Hessian of the quadratic $q(γ)$ given by $DTM−1D$ (Eq. (15)), with $D∈ℝ6M×3Nc$ contact transformation matrix and $M$ diagonal mass matrix. In order to understand the sparsity pattern of $N$⁠, we partition it into 3 × 3 blocks corresponding to each contact. The (ℓ, i)th block of $D$ is nonzero if the ith contact involves body B_ℓ. As a result, the (i, j)th block of $N$ is nonzero if the ith and jth contacts share a body in common.

Let us consider a toy contact problem involving three spherical objects lying on a flat surface as shown in Fig. 2(a). Given the network of bodies at contact in Fig. 2(b), the block sparsity pattern in $N$ corresponds to the edge-adjacency graph in Fig. 2(c); two edge nodes (indexed by body pairs) are connected if they share a body in common. The resulting pattern for this matrix (and thus for $S$⁠) for this toy problem is shown in Fig. 3.

In the context of multibody dynamics, we know the system matrices for the Newton step in Eq. (31) to be sparse and highly structured, and dependent on the iterates $(γ,λ)$ as the PDIP iteration proceeds to the solution of problem (17). Further, from one time-step to the next, we expect the matrix required to obtain the Newton step to change in size as the active set of constraints (corresponding to pairs of objects in contact) evolves. It is because of these changes both within and between timesteps that producing an efficient and robust solution technique for the Newton system remains challenging.

Among currently available hierarchical compression techniques, the TT decomposition features compression and inversion algorithms that are applicable to a large set of structured matrices and that lend themselves to inexpensive global updates. In addition, they have shown to achieve sublinear precomputation times. We thus propose to use it as a framework for direct solution and preconditioning of iterative solvers for the linear systems in each PDIP iteration.

In this section, we give a cursory description of the QTT decomposition as a method to efficiently compress, invert and perform fast arithmetic with approximate representations of structured matrices. We then present a general discussion of its application to solving the linear systems associated with the PDIP for the CCP. This constitutes, to our knowledge, the first application of hierarchical compression solvers to the acceleration of second-order optimization methods. Although it is beyond the scope of this work, we expect the techniques laid out in this section to be readily applicable to a more general class of interior point and other Newton and quasi-Newton based methods for smooth convex problems.

The TT decomposition provides a powerful tensor compression technique [36] via low rank representation akin to that of a generalized singular value decomposition. We will focus on its application in the approximation of tensorized vectors and matrices given a hierarchical subdivision of their indices, known as QTT. Matrices in this setting are further interpreted as tensorized operators acting on such tensorized vectors. We outline how this interpretation allows us to effectively employ the TT as a tool for hierarchical compression and inversion of structured matrices.

Evaluation of Eqs. (33) and (34) on the tensor $A$ is depicted in Fig. 4. The first step (33) produces a rank 2 decomposition $A(i1,i2,i3)=G1(i1,α1)V1(α1,i2,i3)$⁠. In Eq. (34) we then decompose $V1$ to separate dependency of y and z; each of the two terms results in a rank 2 decomposition $V1(α1,i2,i3)=G2(α1,i2,α2)G3(α3,i3)$⁠. We note that each term G_k (depicted in solid yellow, blue and red in Fig. 4) depends on one dimension i_k of the original tensor $A$⁠, and storage has been reduced from n³ to 2n + 4n + 2n = 8n.

What we have presented in this example is an exact tensor train decomposition of $A$⁠. We present the analogous definition for the TT decomposition of a d dimensional tensor.

In order to understand the chain of “low rank” decompositions in the tensor train format in terms of matrix low rank decompositions, we introduce auxiliary objects known as unfolding matrices.

As is the case with low rank matrix decompositions, most often a tensor $A$ of interest will be approximately of low rank, in the sense that given a target accuracy ε, a TT decomposition with low TT ranks $A$ may be found such that $||A−A||F<ε$⁠. This decomposition can be obtained by a sequence of low-rank approximations to $Ak$⁠. A generic algorithm proceeds as in Algorithm 1.

Due to the cost of successive low rank factorizations, implementing Algorithm 1 will predictably lead to relatively high computational cost, O(N) or higher, which is exponential in the tensor dimension d. Instead, we employ TT rank revealing strategies based on the multipass alternating minimal energy (AMEN) cross algorithm of Ref. [36], in which a low-TT-rank approximation is initially computed with fixed ranks and is improved upon by a series of passes through all cores.

We note that, since these algorithms obtain the TT approximation based on an iterative, greedy rank detection procedure, they allow for inexpensive updates to a TT compressed representation of a tensor $A$⁠. If we wish to produce the TT representation of $Ã=A+E$⁠, the AMEN Cross algorithm may be started using the TT representation of $A$ as an initial guess. If $E$ has small TT ranks, this provides a significant speed-up for this compression algorithm, which often converges in a few iterations to the updated representation.

Because AMEN cross and related TT rank revealing approaches proceed by enriching low-TT-rank approximations, all computations are performed on matrices of size $rk−1nk×rk$ or less. In Ref. [16], it is shown that the complexity for this algorithm is thus bounded by O(r³d) or equivalently $O(r3 log N)$⁠, where $r=max(rk)$ is the maximal TT-rank that may be a function of sample size N and accuracy ε. For a large number of structured matrices as well as their inverses, r typically stays constant or grows logarithmically with N [16,37–39]. The overall complexity of computations is then sublinear in N.

Consider a block-sparse, structured matrix $S$ such as the ones obtained from the Newton step system in PDIP applied to contact dynamics. Following the motivating example presented above, our objective is to tensorize matrix $S$ in such a way that it is highly compressible in the TT format (i.e., the resulting TT ranks for a given accuracy ε are reliably small).

For this purpose, we propose a tensorization process obtained from a hierarchical partition of the unknowns corresponding to rows and columns of $S$ based on their location in space. These kind of partitions, usually encoded with tree data structures, have been extensively used in the last two decades to produce multilevel low rank factorizations of structured matrices and associated fast direct solvers. Following Ref. [16], we show that the Tensor train decomposition applied to this tensorized matrix yields an effective hierarchical low rank factorization technique.

In order to demonstrate how this tensorization procedure works, we again consider the small contact problem involving three spherical objects lying on a flat surface presented in Fig. 2. As discussed in Sec. 2.4, the unknowns of the Schur complement matrix $S$ correspond to one normal and two tangent impulses $γi=(γi,n,γi,1,γi,2)$ for each contact point located in space at $ci$⁠.

$S$ can thus be naturally partitioned as a $Nc×Nc$ block-sparse matrix, with the (i, j)th 3 × 3 block encoding the relationship between the change in impulse $Δγj$ and the resulting residual $rγ,i$⁠. We may additionally identify these 3 × 3 blocks as samples from a bivariate function of impulses at each contact point $S(i,j)=K(γi,γj)$⁠; we recall that for interior point methods, $K$ is derived from the augmented Hessian of $f(γ)$⁠.

In order to tensorize this matrix, we first subdivide the physical domain along each coordinate direction (two subdivisions in Figs. 5(b) and 5(c)) and upon each subdivision, we assign a binary index $ik∈{0,1}$ to contact points on each side. In this example, each contact is identified with the leaf node of a binary tree $T$ of depth 2 depicted in Fig. 5(d). A final index $i1∈{0,1,2}$ is required for each of the three coordinates of the impulse vector $γi$⁠. We note that indices are numbered from local to global, matching the order they would appear when tensorizing a column vector. Replacing a single index i with the multiindex $(i1,i2,i3)$ induces a natural tensorization for any vector of samples in our 12 unknowns $(γ1,n,γ1,1,…,γ4,n,γ4,1,γ4,2)$⁠. Unknowns are ordered according to a depth-first (e.g., morton ordering) traversal of the tree, and then the vector is reshaped as a 3 × 2 × 2 tensor.

The Schur matrix $S$ in our example is 12 × 12, with rows and columns dependent on this very set of unknowns; it is depicted in the left-most matrix in Fig. 3. At first glance, it might be tempting to split row and column indices in the order they appear, obtaining tensor of dimension 6 with entries $S(i1,i2,i3,j1,j2,j3)$⁠; this, however, will invariably produce suboptimal TT ranks. Instead, let $TΠ$ denote the product tree $T×T$⁠, corresponding to hierarchical subdivision of both row and column indices. The first bisection produces 4 matrix 6 × 6 blocks (Fig. 3, center), indexed by a block index $b3=i3j3¯$⁠. The second, 16 3 × 3 blocks indexed by $(b2,b3)=(i2j2¯,i3j3¯)$ (Fig. 3, left). Finally, entries on each individual block are identified with the local index $b1=i1j1¯$⁠. This corresponds to ordering matrix entries according to the multiindex $(i1j1¯,i2j2¯,i3j3¯)$ (depth-first traversal on $TΠ$⁠) and reshaping the matrix as a 9 × 4 × 4 tensor.

Given a set of N_c contact points with locations $ci∈ℝ3$⁠, we set a box containing all of them as the root of the hierarchy and proceed to subdivide all nodes along each coordinate direction until each leaf node has at most m points. As was the case for our toy problem, each contact point $ci$ may now be encoded by a set of d indices: a local index $i1∈0,1,…,3m−1$ for the 3 m unknowns and d − 1 binary indices $(i2,i3,…,id)$ corresponding to its location on the tree $T$ at each level. This process is shown on line 1 of Algorithm 2. We note that the binary encoding of points in n dimensional space induced by this coordinate partition induces a so-called Morton ordering (also known as Z-curve order). For efficient algorithms to construct and update such partitions, we refer the reader to Ref. [40].

In general, the tree $T$ will be a uniform binary tree of depth d − 1 with empty nodes corresponding to regions in space with no contacts, and nonempty leaf nodes will contain different numbers of contacts (unless m = 1). Equivalently, every unknown will have an associated multiindex $(i1,…,id)$⁠, but not all multiindices will correspond to a valid unknown. In order to apply TT compression algorithms, we will thus work on an augmented matrix $Ŝ$ of size $2d−1n1×2d−1n1$ (with n₁ = 3 m) such that $Ŝ(i,j)=S(i,j)$ if both i and j correspond to valid unknowns, and $Ŝ(i,j)=δi,j$ otherwise. In other words, $Ŝ$ is a permutation of a 2 × 2 block diagonal matrix with blocks $S$ and identity matrix $I$ corresponding to the original and dummy variables, respectively. This corresponds to the function call in line 2 of Algorithm 2.

and use this routine to obtain an approximate TT factorization S via a greedy, rank-revealing process. We note that each core of $S, Gk(αk−1,ikjk¯,αk)$ depends only on the pair of row and column indices at the corresponding level of the hierarchy. When performing matrix arithmetic, such as matrix-vector product or inversion, TT cores are often reshaped as n_k × n_k matrices parametrized by $αk−1$ and α_k.

In Fig. 3, we illustrate the application of Algorithm 1 to produce the TT decomposition for the Schur complement in our toy contact problem, corresponding to line 6 of Algorithm 2. This method produces a TT decomposition via a sequence of low rank decompositions of unfolding matrices of the tensorized $S$⁠. Because of the way we have indexed $S$⁠, the columns of the first unfolding matrix $S1$ are obtained by vectorizing each of the 16 leaf blocks in $TΠ$⁠. If we apply a row interpolatory decomposition to obtain the low rank factorization $G1V1$⁠, this step of the algorithm finds a low rank interpolation basis of block entries. The algorithm then proceeds by merging children nodes, forming the next unfolding matrix and applying low rank decompositions in an upward pass through the matrix block tree.

As evidenced by the experimental results in Sec. 4, following this procedure reliably yields approximate TT decompositions for $S$ with low TT ranks. This is, first, because it exploits sparsity at different levels of resolution, corresponding to different levels of the tree $TΠ$ (unfolding matrix columns are zero if the corresponding block is zero). As discussed in Sec. 2.4, block-sparsity of $S$ is tied to a network of contact points; a block of $S$ is nonzero if they share a body in common. Unless bodies are substantially different in size and diameter, this network is highly localized: two contacts sharing a body will generally be nearby in space. We note that an interesting avenue of future research might consider tensorizing $S$ following graph partition techniques based on the adjacency graph in Fig. 2(c).

Second, the TT automatically exploits smoothness and symmetries of the bivariate function $K(γ1,γ2)$⁠. The resulting TT decomposition and rank bounds may thus be analyzed in terms of the approximation of $K$ by a sum of separable functions of the form $g1(x1)g2(x2)⋯gd(xd)$ provided by this decomposition, as has been done for integral and differential operators and their inverses.

Fast methods for TT matrix arithmetic are available for a number of operations, including inversion and matrix-vector and matrix-matrix products. The TT solvers in this work find an approximate TT structure for the inverse and employ the corresponding TT matrix-vector product $γ=S−1b$⁠. This mat-vec algorithm proceeds by contracting one dimension of the tensorized matrix $S$^{– 1} at a time, applying the kth TT core. Its complexity can be shown to be $O(r2N log N)$⁠.

TT inversion methods compute an approximate TT decomposition of $S−1$ given the TT decomposition of $S$⁠. In these algorithms, matrix equations for each tensor core of the inverse is solved iteratively. Following an alternating least squares algorithm, given an initial guess for the inverse in TT form, it proceeds by iteratively cycling through the cores (freezing all cores but one) and solving a linear system to update the kth core of S^text--1. TT ranks of the inverse are not known a priori (and are distinct to those for S), and so strategies to increase core ranks are needed to ensure convergence of the alternating least squares procedure to an accurate inverse representation. Further details about variants of this approach can be found in Refs. [41–43]. The complexity of this algorithm for a maximum TT rank r for both S and S⁻¹ is bound by $O(r4 log N)$⁠, as shown in Ref. [16].

We note that, since they share the same iterative, greedy rank detection structure with the AMEN Cross compression algorithm, TT direct solver routines may also be significantly sped up using an approximate initial solution. In the context of the Newton step in interior point methods, we expect both forward and inverse operators to be obtainable via low TT rank updates.

We now discuss how to adapt the tensor train decomposition framework to accelerate computation of the Newton steps within the PDIP method applied to the CCP. We first show how at a given timestep, the tensor train provides an easy-to-update approximate inverse for the Newton system's Schur complement matrix. We then describe a procedure to hot-start the TT compression and inversion algorithms reusing information from the previous timestep, even when the corresponding set of contacts (and thus, matrix size and structure) change. We summarize the entire PDIP step acceleration approach in Algorithm 3.

At a given timestep t, the kth PDIP iteration involves the solution of the Schur complement system, for 3N_c × 3N_c matrix $Sk=N+M̂k−BkEk−1Ck$⁠. In order to accelerate the Newton system solve, we construct approximate TT decompositions S_k and $Sk−1$⁠; we employ the AMEN cross and solve algorithms given user prescribed target accuracy and maximum TT rank parameters (lines 3–8, 14, and 15 in Algorithm 3). We then have the choice to employ $Sk−1$ as a direct solver or couple it with an iterative procedure (e.g., preconditioned Krylov method) if a more accurate solve is needed. For the numerical experiments presented in Sec. 4, we show the latter approach to have the best performance in terms of reducing overall solution costs.

where L_k is block-diagonal. While L_k is generally of matrix rank O(N_c), across all experiments in Sec. 4 we observe it to be of approximate low TT rank, and the same is found for $Sk+1−1$ as a perturbation of $Sk−1$⁠. This fact allows us to use the TT decompositions $Sk,Sk−1$ to hot-start the corresponding AMEN compression and inversion algorithms, reducing precomputation times considerably (∼ 10 × faster in our experiments).

Employing information from the solution of the CCP at a timestep t to hot-start the PDIP iteration at the next timestep t + Δ t is notoriously hard; even if it can be used to produce a feasible point that is close to the optimum (which is nontrivial due to changes in the set of contacts), efforts by the PDIP algorithm to preserve centrality might cause it to take small steps and waste time “returning” to the central path. For this reason, we initialize each PDIP iteration by making the tangential force impulses $γi,1,γi,2$ equal to zero, and the normal force $γi,n$ equal to either a constant preset value (e.g., 1) or to the value computed in the previous timestep if the corresponding contact persists across timesteps. This ensures that our initial value is feasible and lies safely inside the friction cone.

Since pairs of bodies may phase in and out of contact, the set of contacts considered at each timestep (and thus the set of unknowns) as well as the corresponding sparsity pattern both change as the system evolves in time. In Fig. 6 we show an example based on our toy contact problem. While the matrix size for this problem remains the same (since one contact is added and one is removed), we note that the sparsity structure is different due to a change in graph connectivity. However, as long as Δt and relative velocities are sufficiently small, it is likely that the set of persistent contacts (and the corresponding submatrix of $S$⁠) from one timestep to the next will be relatively large. In all validation experiments in Sec. 4, in fact, well over 90% of contacts persist once objects have sedimented.

Our goal is then to re-use the TT decompositions for the initial Schur complement matrix $S0$ from one timestep to the next. From one timestep to the next, contact locations evolve, new active contacts may be introduced and inactive contacts need to be eliminated. We recall from Sec. 3.2 that the current set of active unknowns is considered as a subset of variables corresponding to a uniform binary tree with at most m contacts in each leaf node, and the matrix $S0$ is compressed within an augmented matrix $Ŝ0$ of size $2d−1n1×2d−1n1$ with n₁ = 3 m.

We must then update the indices for the new set of active contacts at time t + Δt

• Indices for contacts that are no longer active must be set to dummy variables and corresponding entries of $Ŝ0$ set to the identity.

• New contacts must be added to the hierarchy and encoded according to their location. Blocks for these contacts and their neighbors (sharing a rigid body) are now nonzero.

• Since binary encoding of contact locations is computationally inexpensive, we may readily detect and update indices for persistently active contacts whose location has significantly changed.

In order to preserve as much structure as possible, we only change the encoding of persistent contacts if there is an indexing conflict, if their current leaf node has more than m contacts, or if their location has migrated significantly. In Fig. 7 we show a simple one-dimensional example of the evolution of tensor indices from one timestep to the next, tracking active nodes in the resulting binary tree structure and corresponding sparsity of matrix $S0$⁠. We note that while N_c changes from 5 to 6, the size of the augmented matrix, depth of the binary tree and number of contacts per leaf (m = 1) for this example remain the same.

If this re-indexing is compatible with the encompasing binary tree, we may use our TT decompositions for $S$₀ and $S0−1$ as initial guesses for the corresponding initial matrix and inverse at the next timestep. Otherwise, we must compute a new hierarchical subdivision and compress the TT factorizations from scratch.

We now demonstrate the performance of the TT-based solver when accelerating the solution of the Newton step system, and thus of second-order methods such as PDIP, in the context of dense, multiple rigid-body dynamics. As mentioned in Secs. 1 and 3, we know that given a desired target accuracy ε, if the associated TT ranks r_k are bounded or slowly growing as a function of problem size N, factorization costs and storage for the TT solver are sublinear in N, and that applying this matrix to a given right-hand-side is $O(N log N)$⁠. We wish to study if Newton system matrices are TT compressible in this sense, and to test their performance and scaling in their solution as problem size grows (determined by the number of collisions N_c). By harnessing the ability of the TT to produce economic updates both from one Newton iteration to the next, as well as across timesteps, we demonstrate significant improvement for the solution of the complementarity problem using second-order methods.

The complementarity method for frictional contact, as well as the solution methods for its CCP relaxation have been thoroughly validated and contrasted with experimental data [12,44]. We setup three experiments based on standard phenomena in terramechanics, focusing on how TT-based linear solvers perform in the PDIP iteration. Models and contact dynamics simulations are performed using open-source library Project Chrono [45], PDIP solver comparisons are performed with a serial matlab implementation and all TT methods are based on the matlab TT-Toolbox [46]. All experiments are run on the serial queue of University of Michigan's Flux computing cluster ().

We first carry out an assessment to determine how the general behavior of TT ranks for the Schur matrix and its inverse vary with target accuracy ε and maximum allowed TT ranks, varying these parameter from 10⁻² to 10⁻⁴ and $r≤10,100,1000$⁠, respectively. This gives us a general idea of the compressibility of these matrices upon re-ordering their nodes according to a spatial hierarchy. Since algorithmic constants for TT compression, inversion and matrix-vector apply all depend on rank, this also informs out choice of ε. We present average results for 1000 PDIP iterations for a sedimentation experiment with M = 35,939 rigid bodies in Sec. 4.2; we note these are largely replicated across experiments, and that ranks and compression times grow very slowly with the number of degrees-of-freedom.

In Table 3 we observe that the Schur matrices and their inverses are indeed highly compressible, and that as we increase the maximum allowed rank, individual and average TT ranks converge. This is replicated in all our experiments, regardless of number of particles and number of collisions. In Table 4 we record the combined matrix compression and inversion times; as indicated in Sec. 3.2, performance of these algorithms is bounded by terms of the form $rk log Nc$⁠. Given that rank growth results in a substantial increase in precomputation times, in our experiments we find that a TT preconditioner approach with ε = 10⁻² is most effective in reducing overall computation for the PDIP iterations. We note, however, that this parameter selection is generally problem and implementation dependent.

We then wish to quantify the impact of the strategies described in Sec. 3.3 to re-use information in TT factorizations for the matrix and its inverse. For this purpose, we run two versions of the TT preconditioned PDIP iteration with and without factorization re-use for 100 timesteps for the same sedimentation problem described above. We set target accuracy to ε = 10⁻² and bound maximum ranks at $r≤10$⁠.

In Fig. 8, we plot precomputation times for the initial Newton steps and for the entire PDIP iteration, in order to measure how effective our information re-use strategies for the TT are in bringing down compression and inversion costs. For the initial PDIP iteration, we see that except for the first timestep (for which both methods have no prior information), information re-use always provides a speedup, between 5 and 15× for most timesteps shown here. The accumulated effect of both re-use strategies can be seen in the second plot on the right; precomputation is improved for all timesteps, resulting in a 10–15× overall speedup. Across all experiments, we observe both an improvement by about an order of magnitude and a drastic reduction in the variance of precomputation times, making the TT approach faster and more robust. We also note that predictably, this speedup is larger for higher target accuracies, as the iterations through all TT cores become more computationally burdensome.

A randomly pertubed cubic lattice of (2n + 1)³ rigid particles of spherical shape of radius 0.1 and friction coefficient μ = 0.25 are dropped and sediment under gravity into a fixed box with a rotating mixer, as shown in Fig. 9. We then run simulations for n = 8, 16, 32, with a total of rigid bodies M = 4915, 35,939 and 274,627 (including the box and mixer), scaling up the box size to keep particle density roughly constant. We set a timestep size Δt = 0.025 and target accuracy 1 × 10⁻⁴ for the PDIP solver, until such time as all objects are deposited in the box and undergoing mixing.

We compute the Newton step through the solution of the corresponding Schur complement system. In order to test performance of the TT-based preconditioner, we compare precomputation and solution times for a preconditioned biconjugate gradient stabilized (BICGSTAB) method against Incomplete LU and unpreconditioned versions of this iterative solver. We note that while use of the conjugate gradient method may be generally preferrable for these systems, we observe its performance may degrade due to ill-conditioning as iterates approach the feasible set boundary, and so for simplicity of presentation we exclude it from our comparison.

As discussed in Sec. 3.3, the Tensor train approach allows us to produce approximate direct and preconditioned iterative solvers by varying target accuracy and maximum TT rank in the compression and inversion processes. Following preliminary parametric studies, we find that setting target accuracy for TT inversion to 1 × 10⁻² and capping maximum TT ranks at $r≤10$ provides the best performance in terms of the trade-offs involved in precomputation and TT preconditioner apply costs for our experimental setup.

We note that for practical purposes, the maximum number of iterations for the unpreconditioned solve was set at 1000; this limit was often reached reducing the overall accuracy of the resulting linear system solution and thus the quality of the PDIP iteration. Through further testing, we consistently observe average BICG iteration counts of 3–5 × 10³ and a fivefold and sixfold increase in computational cost when removing this constraint. These estimates should be considered whenever comparing either of the preconditioned methods against unpreconditioned BICG.

In Fig. 10 we can observe how performance for each of these solvers scales with the number of collisions N_c. Two factors are contributing to increase problem complexity in this experiment: as objects sediment in the box, the number of collisions tends to increase and the Schur complement matrix becomes less sparse. From this plot, we can readily observe that the TT preconditioned solver generally shows superior scaling, outperforming the ILU preconditioner at about N_c ∼ 20,000, and gaining an order of magnitude speed-up against the ILU preconditioner by the end of experiments with M = 35,939, 274,267. We can observe that while the fraction of time spent in precomputation tends to a constant for ILU (indicating similar asymptotic scaling for precomputation and solve times, experimentally $O(Nc3)$⁠, it quickly goes down to zero for the TT.

In Fig. 11, we take a closer look at average precomputation times for each PDIP iteration for TT and ILU. We note that this involves computing roughly 50–100 factorizations, one per Newton iteration. In all experiments, we confirm that compression and inversion times for the Tensor train approach grow extremely slowly with N_c, staying on a range from 5 to 15 s.

Finally, we compare BICGSTAB average iteration counts in Fig. 12. We note that while iteration counts generally increase at the beginning of the experiment, those for the TT grow slower and settle sooner as particles sediment; for M = 35,939, 274,627, the number of iterations for the ILU becomes up to eight times larger. For the unpreconditioned case, the maximum iteration count of 1000 is reached for a significant number of linear system solves, limiting the accuracy of the resulting Newton steps.

We follow the validation experiment in Ref. [12]; as shown in Fig. 13, we setup a drafting test involving a rectangular blade of width 0.1 moving through a container filled with spherical rigid particles of radius 0.1 and friction coefficient μ = 0.25. This test may be used to compute the force that the blade experiences as it moves through the granular flow, which eventually reaches a steady-state. As in the sedimentation test above, we setup a perturbed lattice of (2n + 1)³ spherical particles inside the box, and perform experiments for n = 8, 16 for a total number of M = 4915, 35,939 rigid bodies. The blade moves from one end of the box to the other with a prescribed sinusoidal velocity with period 4 s.