Analysis of Fixing Nodes Used in Generalized Inverse Computation

In various fields of numerical mathematics, there arises the need to compute a generalized inverse of a symmetric positive semidefinite matrix, for example in the solution of contact problems. Systems with semidefinite matrices can be solved by standard direct methods for the solution of systems with positive definite matrices adapted to the solution of systems with only positive semidefinite matrix. One of the possibilities is a modification of Cholesky decomposition using so called fixing nodes, which is presented in this paper with particular emphasise on proper definition of fixing nodes. The generalised inverse algorithm consisting in Cholesky decomposition with usage of fixing nodes is adopted from paper [1]. In [1], authors choose the fixing nodes using Perron vector of an adjacency matrix of the graph which is only a sub-optimal choice. Their choice is discussed in this paper together with other possible candidates on fixing node. Several numerical experiments including all candidates have been done. Based on these results, it turns out that using eigenvectors of Laplacian matrix provides better choice of fixing node than using Perron vector.


Introduction
The motivation of this paper is a computation of a generalized inverse of a symmetric positive semidefinite (SPS) matrix which arises for example in solution of contact problems.One of the possibilities is a modification of Cholesky decomposition.
For the purposes of this paper, we will consider the meshes arising during discretization of numerical prob-lems from graph theory point of view and we will analyse these meshes as graphs.
The particular purpose of this analysis is to find certain nodes (further defined as "fixing nodes") in meshes to reduce numerical instability in computation of a generalized inverse of semidefinite stiffness matrix.
Algorithm presented in this paper is taken from [1].In [1], authors use the Perron vector of an adjacency matrix for the computation of fixing nodes, which is only a sub-optimal choice.In this paper, the proper definition of fixing nodes is presented.Several candidates to fixing nodes are provided based on (spectral) graph techniques together with experimental results.Based on these results, it turns out that using eigenvectors of Laplacian matrix provides better choice of fixing node than using Perron vector (eigenvector of the adjacency matrix).The best choice of fixing node is analysed at the end of this paper.

Generalized Inverse Algorithm
Let us consider a system of consistent linear equations: with an SPS matrix A. In case of FETI method, matrix A is called the stiffness matrix.
If A ∈ R n×n and b ∈ Im A, where Im A denotes the range of A, then a solution x = A + b+Rc of the system of linear equations Eq. ( 1) can be expressed by means of a (left) generalized inverse matrix, A + ∈ R n×n , R is the matrix of rigid body modes, and c is a vector of coefficients.
There are several known algorithms how to compute generalized inverse matrix.There can be used either standard direct methods or some of iterative methods such as Cholesky decomposition, singular value decomposition, and their combinations.
In paper Cholesky decomposition with fixing nodes to stable computation of a generalized inverse of the stiffness matrix of a floating structure [1], authors Dostal, Kozubek, Kovar, Markopoulos, and Brzobohaty present their algorithm of computation of generalized inverse matrix based on Cholesky decomposition, adapted to the solution of systems with SPS matrix.Generalized inverse algorithm [1] consists in a decomposition of the SPS matrix A ∈ R n×n : where r is the number of displacements corresponding to the fixing nodes, and P is a permutation matrix.
Then the generalised inverse A + is computed as: where S † denotes the Moore-Penrose generalized inverse.Permutation matrix P, P = P 2 P 1 , is computed in two steps.First, matrix P 1 is found in the form: where A JJ is nonsingular and A II corresponds to the degrees of freedom of the M fixing nodes.
Second, reordering algorithm on P 1 AP T 1 is applied to get a permutation matrix P 2 which leaves the part A II without changes and enables the sparse Cholesky factorization of A JJ .

3.
Fixing Nodes and Center-like Points

Fixing Nodes
Finding fixing nodes effectively and accurately is the main ingredient of the generalized inverse algorithm.
The active choice is made in sense of permutation matrix P 1 in Eq. ( 4) such that the rows/columns corresponding to those vertices marked as fixing nodes are permuted at the end of the original stiffness matrix (bottom-right block).Also, the rigid body motions corresponding to degrees of freedom of selected fixing nodes are removed.Let us rewrite the definition of fixing node (first written in [2]).Definition 1. Fixing node: Let Ax = b be a system of linear equations arising from a finite element or finite difference discretization of the problem, such that A has one-dimensional kernel (i.e. the singular part A II in Eq. ( 4) is formed by one zero element).
The one-fixing node is the node that makes the regular part A JJ of the stiffness matrix A produced by the permutation P nonsingular and well conditioned, i.e. permutation of this node to the last row/column of the matrix A makes the condition number of the regular part A JJ finite and sufficiently small.
The best choice of one-fixing node is the node k for which the regular part A JJ of the stiffness matrix A produced by the permutation P has the minimal condition number over all A kk (Symbol A kk denotes the principal submatrix arising by removing k-th row and column from original matrix A): Condition number cond( A kk ) is computed as: where λ max ( A kk ) and λ min ( A kk ) denote the largest and the smallest eigenvalue of the non-singular matrix A kk , respectively.Definition of more fixing nodes is composed accordingly.
Definition 2. i-fixing nodes: The i-fixing nodes (if the number i of i-fixing nodes is not important, we use only the term "fixing nodes" instead of the term "i-fixing nodes") are the set of i nodes that make the regular part A JJ of the stiffness matrix A produced by the permutation P nonsingular and well conditioned, i.e., permutation of these nodes to the bottommost rows, rightmost columns respectively, of the matrix A makes the condition number of regular part A JJ finite and sufficiently small (see Eq. ( 4)).
The best choice of i-fixing nodes is the set of i-fixing nodes for which the regular part A JJ of the stiffness matrix A produced by the permutation P has the minimal condition number over all such permutations.Further, we show techniques how to find one fixing node.One fixing node can stabilize well the solution of problems with one-dimensional kernel.Problem of finding more fixing nodes can be reduced to the problem of finding one fixing node by decomposition of the mesh corresponding to given problem by a suitable software into several parts.The fixing nodes are then found one in each part.
In this case, it is guaranteed that we can obtain the best choice of one-fixing node in each part in sense of Definition 1 but the resulting set of i-fixing nodes of the whole problem do not need to be the best choice of i-fixing nodes in sense of Definition 2.
Definition 1 operates on condition numbers of residual matrices, thus it is not suitable to evaluate the fixing node in a reasonable time.Therefore, a fast heuristic has been looked for to provide a good approximation of the fixing node.
As the fixing nodes should lie near the "center" of the mesh [3], there are several possibilities how to define center-like points of mesh, or corresponding graph respectively.

Center-like Points
One of the approaches to find an approximation of the best choice of one-fixing node is based on the mechanical interpretation of the problem.It follows the idea of choice of the fixing node near the "center" of mesh.This translates directly into choosing one of the vertices in the center of the graph.Finding graph centers is not suitable for the numerical solution of large problems but it provides a good referential point to the other methods.This approach has been published in [3].
Another idea how to identify the fixing node is based on spectral approach, namely on eigenvector of the adjacency matrix of the graph.Definition 3. Capital vertex : Let the capital vertex of a graph G be a vertex v i that corresponds to the index i of the highest value in the Perron vector of the corresponding adjacency matrix of graph G.
The capital vertex does not identify the center of the graph, rather a vertex from which most walks of a given length can be realized.By choosing the capital vertex to be the fixing node we expect to achieve a stable numerical solution.Moreover, identifying the capital vertex is fast.
In Fig. 1 overtaken from [1], we can see that the position of the fixing node in the rightmost side subdomain differs from the position where we naturally expect it should be (i.e.closer to the "center" of given subdomain).One can see that the highest value of the Perron vector arises at the vertex with the higher degree rather than in one with the mean value of the degree.
The last presented approach how to obtain a good approximation of the best choice of one-fixing node mentioned in this text is based on the eigenvectors of the Laplacian matrix.
The Laplacian matrix can be set up in dependence on the boundary conditions.As standard Laplacian matrix corresponds to problem with the Neumann boundary, it is more suitable to represent the problem that we are interesting in.
According to M. Fiedler theory [4], [5], the eigenvector corresponding to the second smallest eigenvalue of the Laplacian matrix is used to graph partitioning [6], [7].It is known, that the so-called Fiedler cut defines the cut in certain coordinate direction.The interesting thing is that in two-dimensional case, the eigenvector corresponding to the third smallest eigenvalue defines the cut in the second coordinate direction.Because the Fiedler cut is known to be somehow "optimal", it is reasonable to expect a good approximation of one-fixing node near the crossing of both cuts.
For purposes of definition of a generalized eigenvector cut, the notation based on "Fiedler's tree theorem" is used, especially notation of the characteristic vertex and the characteristic edge.
Using the other Fiedler's results, the Fiedler cut is formed by the characteristic set (i.e.set of characteristic vertices or characteristic edges).In general, characteristic set can be defined for arbitrary eigenvector of the Laplacian matrix and for arbitrary type of graph.
The k-level cut is defined as a characteristic set, i.e. either set of characteristic vertices or set of characteristic edges, where: • characteristic vertices are the vertices i that satisfy • characteristic edges are the edges ij such that i and j are adjacent in Following this notation, the 2-level cut is exactly the Fiedler cut.As the spectrum of graph (or meshes respectively), respects the physical properties of given objects, i.e. the space-dimension, also the crosseigenvector center has to be defined with respect to this characteristic (dimension of the mesh).Definition 5. Cross-eigenvector center : Let G = (V, E) be a graph on n vertices: • For one-dimensional mesh, the cross-eigenvector center is the vertex, edge respectively, that lies on the 2-level cut.
• For two-dimensional mesh, the cross-eigenvector center is the vertex, edge or 2D element respectively, that lies on the crossing of the 2-level cut and the 3-level cut.
• For three-dimensional mesh, the cross-eigenvector center is the vertex, edge, face or 3D element respectively, that lies on the crossing of the 2-level cut, 3-level cut and 4-level cut.
In computational arithmetic, i.e. in presence of rounding errors, there is usually problem with zero identification.The cross-eigenvector center is always assigned to the vertex with the smallest value (in absolute value).The exact position of the cross-eigenvector center can differ no more than half edge, half element respectively.
In Section 5. there is presented the sketch of the proof that the cross-eigenvector center delivers the best choice of one-fixing node (instead of the other centerlike points), which is its biggest advantage.In opposite to the capital vertex, the eigenvectors corresponding to the second, third, and fourth smallest eigenvalue of the Laplacian matrix are harder to compute than the eigenvector corresponding to the largest eigenvalue of the adjacency matrix, which is a disadvantage.

Experiments
The results of positioning of fixing nodes based on listed definitions are presented in this section.The criterion to measure quality of the approximation of one-fixing node is the condition number of the regular part A JJ of the matrix A of given problem Eq. ( 4) after permutation of the row and column corresponding to the fixing node at the end of the matrix.
Given examples are quite small (less than 121 vertices and 100 edges) because the finding of the best choice of one-fixing node according to Definition 1 yields to removing all vertices, computing the condition number of the regular part of the matrix (i.e. the largest and the smallest eigenvalue) for each removal, and choosing the minimal one which is a huge timeconsuming process.
In figures the following symbols appear: • The best choice of one-fixing node satisfying Definition 1 is drawn as a circle .
• The graph centers are drawn as a square .
• The capital vertex satisfying Definition 3 is drawn as a triangle .
• The cross-eigenvector center satisfying Definition 5 is drawn as a star * .
Let us have a look in the rightmost subdomain in Fig. 2. As the cross-eigenvector center (color cyan) almost the same vertex as the best choice of one-fixing node (color green) has been detected, meanwhile the capital vertex (color red) has been assigned to the vertex with higher degree which is far from the best choice of one-fixing node.Let us briefly present the result of this phenomenon.
At first, let us have a look on the behaviour of centerlike points of the meshes with regular elements, see Fig. 3(a).In two-dimensional space, the element corresponds to the graph face (region bounded by edges).Our examples consist of quadrilateral type elements only.As we can see from Fig. 3(a) and from the Tab. 1, almost all approaches deliver a good approximation of one-fixing node.
One can suppose that the behaviour of all center-like points is similar.But this holds only for meshes with uniform elements, or similarly, for meshes with inner vertices of the same degree.Let us have a look on the mesh with mixed elements (quadrilateral and triangle type) for further testing.This means that the degree of inner vertices differs from four to six.
For well understanding of the problem, several eigenvector of testing meshes are plotted (compare Fig. 3, Fig. 4).In Fig. 3  The testing mesh (Example 1) is similar to the Cartesian product of two paths P 11 × P 11 and also, the spectrum of this mesh is similar to the spectrum of Cartesian product P 11 × P 11 .The mesh Example 1 can be understood as an approximate Cartesian product with deleted both edges and vertices.The missing left upper corner has no essential influence on the behaviour of the eigenvectors.
Another situation arises if the edges are added instead of deleted (see Fig. 4).The behaviour of the eigenvector of the adjacency matrix corresponding to the highest eigenvalue does not correspond to the classical behaviour of the spectrum of Cartesian products.Therefore, the capital vertex (based on the adjacency matrix) do not approximate the best choice of onefixing node well.
As we can see in Fig. 4(c), Fig. 4(d), the eigenvectors of the Laplacian matrix are still indifferent to change of the structure of the mesh.The cross-eigenvector center still approximates the best choice of one-fixing node well.
For better comprehension, Tab. 1 is presented.The first column represents the name of the test example.The regular condition number of the matrix A in the second column is computed by: where λ max (A) and λ min (A) correspond to the largest and to the smallest non-zero eigenvalue of A. The regular condition number is used as a reference value, because it represents the smallest boundary to the condition number of the generalized inverse (the condition number of the generalized inverse can be never smaller than the condition number of the original matrix).The third column represents the condition number of the matrix A JJ considered in Definition 1, i.e. the minimal possible condition number of the regular part A JJ after removing the best choice of one-fixing node.The condition number of the regular part A JJ after removing the graph center (without definition) is written in the fourth column.
If more vertices satisfy the definition, the vertex that causes the minimal condition number has been chosen.The fifth column represents the condition number of the matrix A JJ when the capital vertex satisfying Definition 3 is removed.The last column represents the condition number of the matrix A JJ when the crosseigenvector center satisfying Definition 5 is removed.
We can see that choosing the capital vertex as the approximation of the best choice of one-fixing node instead the cross-eigenvector center debases the size of the condition number from 426.5 (best choice of onefixing node, cross-eigenvector center) to 537.1, i.e. by 26 %.

The Best Choice of One-Fixing Node
Based on the experiments, it is seemed that the best approximation of the best choice of one-fixing node is the cross-eigenvector center.In this section, the particular theorem is presented.
Theorem 6.The best choice of one-fixing node: The cross-eigenvector center (according to Definition 5) is the best choice of one-fixing node (according to Definition 1).
I.e. if we remove the row and column corresponding to the cross-eigenvector center from the original matrix A, the remaining principal submatrix has the best condition number over all principal submatrices.
Remark that we restrict on the matrices with onedimensional kernel, i.e. the resulting submatrix has deleted one row and one corresponding column.In above definitions, we have worked with some general stiffness matrix A. As we consider the Laplacian matrix of graph, we will use the symbol L hereinafter.Symbols with ∼ as L, λ, etc. will further denote variables of corresponding reduced problem.If we could emphasize that the i-th row/column is removed, we assign the corresponding variables as L i , λ i , etc.
The condition number is considered in the form: From one of the variant of Cauchy-like interlacing theorems (see i.e. [8]), removing the fixing node does not change the λ max so much as the λ min .Therefore, the minimization of the condition number corresponds to the maximization of the λ min .Theorem 6 can be paraphrased into following form: Theorem 7. Maximization of λ min : Removing the vertex k corresponding to the cross-eigenvector center in Laplacian matrix L, the maximal value of λ min of L is obtained over all λ i min of all principal submatrices L i , i.e., the condition number L i is minimized over all i.
Finding the maximal value of λ min can be written: As the whole proof of this theorem consists in comprehensive analysis of all cases, we restrict on the main ideas only.The complete proof can be found in [2].
Proof.In our considerations, we come out from the Rayleigh's principle [9].The smallest eigenvalue of the reduced matrix L can be computed as: As we would like to find which vertex i from the original matrix L should be removed to obtain the maximal value of λ min , we have to rewrite the above equation using original matrix L. Remark that the removing of the vertex i in original matrix L corresponds to setting the corresponding component of vector x to zero, which can be written as: In the last equation, we have substitute the Laplacian matrix by its spectral decomposition, where Λ is a diagonal matrix of eigenvalues of L ordered in ascending order and Q is the matrix, whose columns are eigenvectors of L ordered accordingly.
After some manipulations (and substituting y = Q T x, y T = x T Q respectively) we get: where r i (Q) denotes the i-th row of matrix Q (and it is not the eigenvector).
For given n we solve following system of equations: subject to: q y 1 + q i2 y 2 + • • + q in y n = 0.
Here, q ij denotes the ij-th entry of matrix Q, i.e. q i denotes again the i-th row of matrix Q.
Thanks to Cauchy Interlacing Theorem [8]: we get the upper bound for λ i min ≤ λ 2 .Considering the particular example of regular rectangle mesh with known eigenvalues and eigenvectors, it can be shown, that we obtain exactly the upper bound for i such that q i2 = 0, [2]: Notice that in this notation, the q i2 corresponds to the i-th entry of the second smallest eigenvector of the Laplacian matrix which corresponds exactly to 2-level cut according to Definition 4.

Conclusion
In this paper, we have presented a modification in stabilization of generalised algorithm used in [1].Our modification consists in another choice of fixing node.
Several numerical experiments have be done to observe the influence of different choice of fixing nodes on the condition number of resulting matrix.Based on these experiments, it been shown that using eigenof the Laplacian matrix provides better choice of fixing node than using eigenvector of the adjacency matrix (used in [1]).
We presented the main ideas of the analytical proof that the capital vertex (defined using eigenvectors of the Laplacian matrix) provides the best choice of one-fixing node.The whole proof is based on knowledge of eigenvalues and eigenvectors of the particular example of regular rectangle mesh and it can be found in [2].
(b), Fig. 4(b), the eigenvectors of the adjacency matrix corresponding to the highest eigenvalue are plotted.In Fig. 3(c), Fig. 4(c), the eigenvectors of the Laplacian matrix corresponding to the second smallest eigenvalue are plotted, and in Fig. 3(d), Fig. 4(d), the eigenvectors of the Laplacian matrix corresponding to the third smallest eigenvalue are plotted.
Tab. 1: Approximation of the best choice of one-fixing node.