A graph-theoretic approach to sparse matrix inversion for implicit differential algebraic equations

In this paper, we propose an e fficient numerical scheme to compute sparse matrix inversions for Implicit Differential Algebraic Equations of large-scale nonlinear mechanical systems. We first formulate mechanical systems with constraints by Dirac structures and associated Lagrangian systems. Second, we show how to allocateinput-output relationsto the variables in kinematical and dynamical relations appearing in DAEs by introducing an oriented bipartite graph. Then, we also show that the matrix inversion of Jacobian matrix associated to the kinematical and dynamical relations can be carried out by using the input-output relations and we explain solvability of the sparse Jacobian matrix inversion by using the bipartite graph. Finally, we propose an e fficient symbolic generation algorithm to compute the sparse matrix inversion of the Jacobian matrix, and we demonstrate the validity in numerical e fficiency by an example of the stanford manipulator.


Introduction
Multibody systems such as space structures, manipulators, etc. are known to be represented as implicit mechanical systems with kinematical constraints, holonomic or nonholonomic, which may be eventually expressed by implicit nonlinear Differential-Algebraic Equations (DAEs).In particular, for the numerical integration of such DAEs, we need to employ stiffly stable implicit numerical integrators such as Backward Differentiation Formula (see, Hachtel et al., 1971;Brayton et al., 1972), since the DAEs are to be highly nonlinear and stiff in general.On the other hand, one may face at a serious problem in CPU time for solving the implicit nonlinear algebraic equations, especially, for the case of large-scale systems.Namely, increasing degrees of freedom of the system, it eventually requires much CPU time in computing the matrix inversion of Jacobian matrix of the implicit DAEs in Newton's iteration at each time, since the Jacobian matrix of discretized nonlinear algebraic equations may be random sparse in general.
A major stumbling block lies in the fact that the Jacobian matrix has the random sparseness as well as highly nonlinear in terms of generalized coordinates.So far, some nu-merical technique of sparse matrix inversions for VLSI circuits or networks has been developed by using the blocktriangularization of matrices (see, for instance, Orlandea et al., 1977a;Murata et al., 1985), where a structural analysis is effectively made by means of graph and matroid theory (for instance, refer to Murota, 2000).In these conventional approaches, one may properly find out pivots in the Gaussian elimination process at each time step in an ad hoc way (where we note that the choice of pivots is quite relevant with input-output relations as will be shown shortly).This eventually requires much CPU time to calculate the inversion of the Jacobian matrix unless utilizing some effective sparse matrix algorithms.Namely, it is almost impossible to figure out a prior fill-in and fill-out in Gaussian elimination since topological structure of such a sparse Jacobian matrix might be so much random and complicated.Thus, we need to develop an efficient numerical algorithm of sparse matrix inversion for solving large-scale implicit DAEs in a systematic way.
In this paper, we develop a graph-theoretic approach to computing sparse matrix inversion for large-scale nonlinear mechanical systems by using Dirac structures and associated implicit Lagrangian systems.The underlying idea is to regard a mechanical system as an interconnected system of Published by Copernicus Publications.
elements, in which systems are comprised of constitutive relations of physical elements, structural relations among the physical relations, and causal (input-output) relations among physical variables.In particular, focusing upon the inputoutput relations associated with all the kinematical and dynamical relations of original mechanical systems, we develop bipartite graphs and then we show how the sparse matrix inversion can be made by effectively using the input-output relations.Furthermore, we explain solvability for the sparse Jacobian matrix inversion associated to the DAEs by using the bipartite graph.Finally, we propose an efficient and systematic symbolic generation algorithm to compute the sparse matrix inversion of the Jacobian matrix and we demonstrate its validity in numerical efficiency by an example of the stanford manipulator.

Implicit Lagrangian systems
Let us review Dirac structures and associated implicit Lagrangian systems by following Yoshimura andMarsden (2006a,b, 2008).

Dirac structures
Let Q be an n-dimensional configuration manifold, whose kinematical constraints are given by a constraint distribution where ω a are m one-forms on Q. Define the distribution ∆ T * Q on T * Q by Let Ω be the canonical symplectic structure on T * Q and Ω :

Local representations
Let us choose local coordinates q i on Q so that Q is locally represented by an open set W ⊂ R n .The constraint set ∆ Q defines a subspace of T Q, which we denote by ∆(q) ⊂ R n at each point q ∈ W. If the dimension of ∆(q) is n − m, then we can choose a basis e m+1 (q), e m+2 (q), . . ., e n (q) of ∆(q).Recall that the constraint sets can be also represented by the annihilator of ∆(q), which is denoted by ∆ • (q) spanned by such one-forms that we write as ω 1 , ω 2 , . . ., ω m .Using π Q : T * Q → Q locally denoted by z = (q, p) → q and T π Q : T T * Q → T Q; (q, p, q, ṗ) → (q, q), it follows that Let points in T * T * Q be locally denoted by (q, p, β, u), where β is a covector and u is a vector.Then, the annihilator of ∆ T * Q is locally represented as Since we have the local formula Ω (q, p)•w (q,p) = (q, p, − ṗ, q), the condition α (q,p) − Ω (q, p) • w (q,p) ∈ ∆ • T * Q reads α + ṗ ∈ ∆ • (q), and w − q = 0, where α (q,p) = (q, p, α, w) and w (q,p) = (q, p, q, ṗ).Thus, the induced Dirac structure is locally represented by Representation (I): let us introduce a matrix representation of D ∆ Q given in Eq. ( 2).First, let N T (q) be an n × m matrix whose m-column vectors ω 1 (q), ..., ω m (q) span the basis of ∆ • (q), namely, N T (q) = [ω 1 (q), ..., ω m (q)] and the distribution ∆(q) ⊂ R n T q Q may be represented by Using Lagrange multipliers λ = (λ 1 , ..., λ m ) ∈ R m , one has Thus, the induced Dirac structure can be represented by Representation (II): as shown in Eq. ( 3), for Representation (I) for the induced Dirac structure, we utilized the Lagrange multipliers, which represent constraint forces in constrained mechanical systems.Here, we develop another representation of D ∆ Q on T * Q without using the Lagrange multipliers.
Let us choose an n × (n − m) matrix B(q) = [e m+1 (q), . . ., e n (q)], whose column vectors span the basis of ∆(q).Then, it follows that the distribution ∆(q) ⊂ R n T q Q can be also represented by where u = (u m+1 , ..., u n ) ∈ R n−m .Note that the orthogonality condition between N(q) and B(q) holds: The above condition naturally comes from the fact that ∆ • is the annihilator of the distribution ∆; namely, in other words, the basis e m+1 (q), ..., e n (q) is orthogonal to the dual basis ω 1 (q), ..., ω m (q) at each q ∈ Q.Therefore, one can read that Thus, the induced Dirac structure D ∆ Q ⊂ T T * Q ⊕ T * T * Q can be represented without using the Lagrange multipliers as

Ehresmann connection and structural relations
We briefly review an Ehresmann connection associated with nonholonomic mechanical systems; as to the details, for example, refer to Yoshimura and Marsden (2006b).
Assume that there is a bundle structure with a projection π : Q → R for Q; that is, there exists another manifold R called the base.We call the kernel of T q π at any point q ∈ Q the vertical space denoted by V q .An Ehresmann connection A is a vertical vector-valued one-form on Q, which satisfies Thus, we can split the tangent space at q such that T q Q = H q ⊕ V q , where H q = Ker A q is the horizontal space at q. Suppose there exist nonholonomic constraints ∆ Q ⊂ T Q, which are given by m (< n) algebraic equations for n generalized velocity vector v = q = ( q1 , ..., qn ) ∈ ∆(q) ⊂ T q Q as in Eq. ( 1).Let us choose an Ehresmann connection A such that H q = ∆ Q (q) or we assume that the connection is chosen such that the constraints are written as A • v q = 0, where the constraint distribution ∆ Q is spanned by a set of m independent one-forms, which is given, in local coordinates q i = (r α , s a ) for Q, by ω a = ds a − J a α (r, s)dr α .
In a matrix representation, where v is locally split into dependent velocity v • = q• = ( q1 , ..., qm ) and v * = q * = ( qm+1 , ..., qn ) independent velocity and J is a submatrix associated with the constraints.Geometrically speaking, this splitting corresponds to a choice of Ehresmann connections for the given constraints (see Yoshimura and Marsden, 2006b).
Corresponding to the annihilator, one has the dynamical relations associated to the generalized force vector where ∆ • Q denotes the annihilator of ∆ Q , and the generalized force vectors dual to v • and v * respectively.On the other hand, the input-output relation between Q • and Q * is reverse to v • and v * ; namely, Q • is the input and Q * the output.In the above, the orthogonality condition holds: The matrices N(q) and B(q) are called connection matrices (see Yoshimura, 1995).This orthogonality condition denotes principle of virtual work, which is given by where , denotes a duality pairing.
The dual set of constraints given by Eqs. ( 4) and ( 5) indicates structural relations, namely, it represents how physical elements are interconnected.Thus, we sometimes call the structural relation an interconnection among the physical elements.In circuit theory, it is known that Eqs. ( 4) and (5) correspond to KCL and KVL and also that the virtual work principle is known as Tellegen's theorem.Furthermore, there exists a relation called duality principle as in Fig. 1, which is known as Planck-Okada-Arsove principle in circuit theory (see Yoshimura, 1995).

Implicit DAEs for Lagrangian systems
Here, we show how the notion of Dirac structures can be fit into the formulation of implicit Lagrangian systems (see Yoshimura andMarsden, 2006a,b, 2008).Let L be a Lagrangian on T Q, which is given by where we assume that L is only associated to q• , M is a mass matrix whose components are functions of q • , and U denotes a potential energy function of q • .This implies that L is possibly degenerate.
H. Yoshimura: A graph-theoretic approach to sparse matrix inversion The Lagrange-d'Alembert principle is given by b a d dt where δq satisfies the constraint So, one can obtain the dual dynamical relation where and it directly induces Notice that τ indicates the external forces.Of course, equations of motion can be written as Furthermore, one has the second-order condition (see Marsden and Ratiu, 1999): From Eqs. (4), ( 5), ( 6), ( 7) and ( 8), we can obtain the following local differential-algebraic equations where In the above, we locally set space which is a model space for Q.Thus, the mathematical model of the Lagrangian system is given by the implicit DAEs: 3 Sparse tableau approach From the viewpoint of numerical analysis for mechanical systems, there exist two kinds of dynamical problems; namely, the forward dynamics and inverse dynamics.
Recall that W = W m ×W n−m = R m ×R n−m is the model space for Q.The forward dynamics analysis is the case in which given a smooth input vector as a vector function of time t, numerically integrate Eq. ( 9) in terms of t to obtain as the output, where x ∈ W × W × W * .On the other hand, the inverse dynamics analysis is the case in which given a smooth input vector as a vector function of time t, then compute In this paper, we explore the case of the inverse dynamics analysis by the sparse tableau approach, where the state vector is given by x To do this, let us first discretize Eq. ( 9) at time t = t n .By using the Backward Differentiation Formula (BDF) developed by Gear (see, for instance, Brennen et al., 1995), the time-derivative term ẋn = ẋ(t n ) associated to the state vector x n = x(t n ) may be approximately discretized by the backwards where h = t n − t n−1 denotes a time step, k is a backward order and α i indicates the coefficient associated to the i-th backward order.Substituting Eq. ( 11) into Eq.( 9), we obtain the discretized nonlinear algebraic equations as follows: Recall the algorithm of the Sparse Tableau Approach is given in Fig. 2 (see Hachtel et al., 1971;Brayton et al., 1972), where we linearize Eq. ( 12) at each time step t = t n as and where J = [∂G i /∂x j ] t=t n is the Jacobian matrix and x (r) n denotes the r-th iterated corrector vector at t n .Then, it follows x n G( x n ) (r) x n J 1 G( x n ) x n (r) x n (r) x n (r) x n " In the Newton method, it is necessary to take initial values near from the solution x n and the k-th prediction formula x Pr n is given by where γ i is the i-th coefficient.
In the inverse dynamics analysis, the state vector is given by x ) and it follows from Eq. ( 10) that the Jacobian matrix is given as in Fig. 3, where I n stands for the n-th degree unit matrix.

Input-output relations
The Jacobian matrix J(x (r)  n ) obtained in Eq. ( 13) apparently has the characteristic of random sparseness.So, we develop an efficient symbolic generation for computations of the sparse Jacobian matrix inversion for Newton's iterations.To do this, let us consider an input-output relation among state variables for every relation in ( 10).Now, we can uniquely allocate the input-output relation to the kinematical and dynamical relations in Eqs. ( 4) and ( 5) as follows: Similarly, for the second-order conditions between v = (v • , v * ) and q = ( q• , q * ), one has Furthermore, as to the equations of motion, it follows The input-output relations in the mentioned above can be determined by physical causality.Needless to say that the time derivative terms ẋ(t n ) are expressed in terms of the backwards x n−i (t n ), i = 0, ..., 6 by using the BDF as in equation ( 11).Corresponding to the Jacobian matrix J in (13), define the causal Jacobian matrix Ĵ by assigning −1 to the input and +1 to the output as to the j-th variable in the i-th relation associated to −1 : the j-th variable is input, +1 : the j-th variable is output, 0 : otherwise.
Thus, the causal Jacobian matrix Ĵ is given in Fig. 4.
In the above, I m indicates the m-th unit matrix.Note that there exists an element with +1 in each row, which plays a role of the pivot in the Gaussian elimination.Further, I m,n−m is the m × (n − m) matrix, in which +1 are allocated to non-zero components of the m × (n − m) matrix J.

Oriented bipartite graphs
Let us illustrate the input-output relations as to the Ja- bipartite graphs.Recall that the i-th row of the Jacobian matrix J = [J i j (x)] = [∂G i /∂x j ] t=t n corresponds to the i-th vector function G i (x, ẋ; u) and the j-th column corresponds to j-th vector x j respectively, where i, j = 1, ..., N(= 3n).Here, G i (x, ẋ; u) = 0 are the N-th order DAEs and x = (x 1 , ..., Given the causal Jacobian matrix Ĵ associated to J, let V r = {G 1 ,G 2 , . . .,G N } be the row-set of Ĵ and let V c = {∆x 1 , ..., ∆x 6 } be the column-set of Ĵ .Here, we assign V r and V c to vertex set V := V r × V c .Now, there are nonzero elements (±1) in the i-th row of the causal Jacobian matrix and suppose that they are in the k-th and l-th columns.For instance, as to the first equation G 1 = 0, one has the kinematical relation between x 1 = q • and x 2 = v • as ẋ1 − x 2 = 0.By the BDF discretization, the corresponding linearized equation is given as In order to illustrate input-output relations among the state variables by graphs, let us define arc-set by where the arcs represent some relations between vertices in V c (state variables) as to every vertex (equation) in V r .Furthermore, the direction of an arc indicates causality or an input-output relation among state variables associated to the column vertices.Let a ∈ A be an arc.Let us introduce a map s : A → V, which is given by, for an arc a ∈ A, s(a) denotes the initial vertex of a.Similarly, define a map t : A → V, which is by, for an arc a ∈ A, t(a) indicates the terminal vertex of a. So, we can define the set of vertices incident to a by {s(a), t(a)}.Sometimes, s(a) is called the source of a and t(a) the target of a. Hence, we can define a directed bipartite graph by by which we can illustrate the input-output relations associated to the Jacobian matrix as shown in Fig. 5.

Perfect matching
A matching M ⊂ A of a bipartite graph B = (V r , V c , A, s, t) is defined by a set of arcs without common vertices.Especially, M is called a perfect matching if it is a matching which matches all vertices of the graph, namely, every vertex of the graph is incident to exactly one edge of the matching: s(M) = V r and t(M) = V c .In Fig. 6, the arcs that are drawn by the broken lines consist of the perfect matching.If a perfect matching exists in B, the corresponding Jacobian matrix J is a square and nonsingular matrix since |V r | = |V c |. Here, |V| means the size of V, namely, the number of elements in the set of vertices V.
Let us see that |V + | = |V − | is equivalent with the fact that the Jacobian matrix is diagonalizable.First, let us show the following relation: This is clear because a perfect matching can be detected by the elementary row and column operations in matrix; namely, (1) interchanging two rows or columns; (2) adding a multiple of one row or column to another; (3) multiplying any row or column by a nonzero element, although there might be several perfect matchings for a given bipartite graph.Next, let us show the following relation: That the Jacobian matrix J is diagonalizable means that one can properly select a pivot in each row or column from nonzero elements of the matrix.In other words, the rank of J is equal to N. Recall that a cover is defined as a pair ( Vr , Vc ) of Vr ⊂ V r and Vc ⊂ V c such that there exist no arcs between V r \ Vr and V c \ Vc .For the case in which Rank J = N = |V|/2, the number of the minimum cover for the vertex set V = (V r , V c ) of the bipartite graph B is to be N.It follows from the König-Egerváry theorem that the number of arcs in a maximum matching is equal to the number of vertices in a minimum vertex cover; namely, So, the size of the maximum matching is to be N = |V|/2 and therefore we have shown that the matching is a perfect matching.As shown in Fig. 6, there exists a perfect matching For each arc a ∈ M, one can choose the target x = t(a) as the pivot associated to the k-th equation, where k = s(a) is the source of a. Thus, the proof has been done.Therefore, we can conclude the following relation: Perfect matchings exist in ⇔ J is diagonalizable.
this way, the solvability of the linearized Eq. ( 13) is clarified in the context of the perfect matching by using the bipartite graph associated to the causal Jacobian matrix.Recall that nonzero elements of Ĵ correspond to arcs of the bipartite graph and hence the number of nonzero elements and the number of the arcs are the same.In our study, in each row, there is only one output variable, which implies that the arcs associated with the output variables never share the vertices of other output arcs in the bipartite graph.Therefore, the set of edges of output variables has perfect matching.

Symbolic generation
Next, we show how the sparse matrix inversion of the Jacobian matrix can be done by symbolic manipulation.By the information from the causal Jacobian matrix, selecting the pivot, the Jacobian matrix can be diagonalized as shown in Fig. 7, where G (k) is the k-th tableau error vector, A (k) represents a matrix in the Jacobian matrix for the k-th step of elementary operations.Next, after forward Gaussian elimination for the Jacobian matrix, we can easily obtain the reduced Jacobian matrix as in Fig. 8.Note that we can do this by symbolic manipulation.Furthermore, after backward elimination as to the reduced Jacobian matrix by symbolic manipulation,  we can obtain the following corrector vector consequently: 3 ), ∆v * = G (1)  1 .
Thus, we can develop the inversion of the Jacobian matrix, i.e., ∆x (r) n = −J −1 (x (r) n )G(x (r) n ) can be explicitly done by symbolic manipulation with the causal information of Ĵ.As a result, we can make symbolic generation of Eq. ( 14).

Numerical analysis
Let us demonstrate the validly of our symbolic generation scheme by an example of the Stanford manipulator with 6 degrees of freedom as shown in Fig. 9.One can set up 96, 150, 210, 390, 738 system equations of DAE models for the stanford manipulator, where the numbers of systems equations correspond to those of the unknown variables according to the elimination of redundant variables.
We examine to compare the CPU time required in the routine of the Jacobian matrix inversion in each time step for the three cases: (1) Ordinary Gaussian elimination method (without any sparse matrix algorithm), (2) Sparse Gaussian elimination method of the inner product algorithm (as to the opensource subroutine, see Murata et al., 1985), (3) our method.The comparison of the required CPU time is illustrated in has a great advantage in the CPU time efficiency in comparison with other approaches.

Conclusions
We have shown symbolic generation of sparse matrix inversion for large-scale mechanical systems.We have set up Implicit DAE models in the context of Lagrangian and we have shown the input-output relations as to the Jacobian matrix associated with linearized equations.Then, we have shown the solvability of the linearized equations by using the bipartite graph.Furthermore, we have proposed symbolic generation of the random sparse matrix inversion for the Jacobian matrix.Finally, we have demonstrated the validity of our approach by numerical analysis with an example of the Stanford manipulator comparing with the inner-product sparse matrix Gaussian elimination algorithm as well as the standard Gaussian elimination algorithm.