On Moore-Penrose Pseudoinverse Computation for Stiffness Matrices Resulting from Higher Order Approximation

Computing the pseudoinverse of amatrix is an essential component ofmany computationalmethods. It arises in statistics, graphics, robotics, numerical modeling, and many more areas. Therefore, it is desirable to select reliable algorithms that can perform this operation efficiently and robustly. A demanding benchmark test for the pseudoinverse computation was introduced. The stiffness matrices for higher order approximation turned out to be such tough problems and therefore can serve as good benchmarks for algorithms of the pseudoinverse computation. It was found out that only one algorithm, out of five known from literature, enabled us to obtain acceptable results for the pseudoinverse of the proposed benchmark test.


Introduction
The Moore-Penrose pseudoinverse arises in many real life applications.It is commonly used in graphics (image recognition [1], image restoration [2], data compression [3]), robotics [4], and statistics [5] as well as diverse branches of engineering [6,7].Our examination of various algorithms for the pseudoinverse was motivated by the necessity of using such an operation in the local numerical homogenization in solid mechanics.This method, proposed by Jhurani and Demkowicz (see [6][7][8] and also [9]), is one of the methods of the numerical upscaling, i.e., such a coarse mesh discretization that incorporates a fine mesh information.Other discretization-based homogenization methods are presented, for example, in [10][11][12].One of the key steps of the local numerical homogenization is a computation of the pseudoinverse matrix.We observed that whenever the order of approximation at any level is greater than one, the behavior of the algorithms for the pseudoinverse was very poor.We have tested five approaches to computation of the pseudoinverse and concluded that the Demmel-Kahan algorithm [13] was the best one.It is worth mentioning that all the methods performed very well for many other benchmark problems [7,8,[14][15][16][17].Therefore, a stiffness matrix for elasticity problem is a really demanding test for the pseudoinverse algorithm and can be used to verify the competitiveness of the newly developed methods.
In the next section the local numerical homogenization for solid mechanics is recalled.Then, several methods for the Moore-Penrose pseudoinverse computation are briefly described.Finally, numerical results are presented and discussed.

Local Numerical Homogenization
In this paper, we deal with the finite element (FE) analysis of heterogeneous materials.Let us assume that due to limited computational resources fine mesh solutions to problems with rapidly varying material parameters are infeasible.Instead of that, we generate a coarse mesh, which is, in fact, the finest mesh we can computationally afford (Figure 1).At this level, homogeneity of the domain is assumed and the mesh is adapted at the macroscale using automatic mesh refinement.We use for this purpose the Dhp code [18].Subsequently, the mesh is refined independently for every coarse element in order to comply with the material heterogeneity.The core of the local numerical homogenization is the computation of effective stiffness matrices for the macroelements.After that, an equivalent coarse mesh problem can be solved.At the postprocessing level also a fine mesh solution for each of coarse elements can be obtained.
A mathematical formulation of the effective coarse mesh stiffness matrix computation is recapitulated below following [6,8].All of the quantities referring to the macroscale are marked with hat symbol (e.g., x).
Given symmetric fine mesh element stiffness matrices assembled into  ∈ R 푁×푁 , a nonzero fine mesh load vector  ∈ R 푁 , an interpolation matrix  ∈ R 푁×푀 , a positive-definite symmetric weight matrix  ∈ R 푁×푁 , and a dimensionless small parameter  > 0 find a symmetric matrix K † ∈ R 푀×푀 for which  reaches minimum, where and ‖‖ 2 = √ ( 푇 ) denotes the Euclidean norm, ‖‖ 퐵 = √ ( 푇 ) denotes the Euclidean norm weighted with , and ‖‖ 퐹,퐵 = √ ( 푇 ) denotes the Frobenius norm weighted with . and  are arbitrary column vector and matrix, respectively.The numbers of degrees of freedom for the fine and coarse meshes are denoted with  and .
The expression for  consists of two terms.The first of them measures the error introduced to the local solution vector by the numerical homogenization, whereas the second term is a regularization one and guaranties the solution uniqueness.The matrix K † that minimizes (1) is the pseudoinverse of the local effective stiffness matrix K. Derivation of (1) is presented in [6,8].Imposing the symmetry conditions on K † ( K † − ( K † ) 푇 = 0) into the Lagrangian, it was proved that problem (1) is equivalent to the continuous-time Lyapunov equation in a following form: with auxiliary matrices In (1)  † denotes the Moore-Penrose pseudoinverse of matrix  [20,21].In the case of the real matrices we call a matrix  † ∈ R 푁×푀 a pseudoinverse of matrix  ∈ R 푀×푁 , if it fulfills the following conditions: A reliable computation of  † is the key point of the local numerical homogenization and it needs to be computed twice (for each coarse element): (i) for the assembled stiffness matrix  ∈ R 푁×푁 ; (ii) for the pseudoinverse of the effective stiffness matrix ( K † ∈ R 푀×푀 , where  ≪ ).
A comprehensive review of the theory and computational aspects of the pseudoinverse computation can be found in [14], in its updated version [19], and in the most recent book [22].Also Jhurani [7,8] presented several methods dedicated to the local numerical homogenization.
In the following chapters selected effective algorithms for the Moore-Penrose pseudoinverse matrix are described and numerically verified.Due to our field of interest (pseudoinverses of stiffness matrices), we limit the forthcoming discussion to real matrices ( ∈ R 푀×푁 ).In this paper, a group of pseudoinverse computation methods based on neural networks has not been covered.This approach can be found in papers [23,24].

Computing Pseudoinverse Using Singular Value Decomposition (SVD)
A well established and the most reliable algorithm for the pseudoinverse matrix computation is the one based on the singular value decomposition.Every matrix  ∈ R 푀×푁 can be decomposed into  = UΣV 푇 , where  ∈ R 푀×푀 and  ∈ R 푁×푁 are orthogonal matrices and Σ ∈ R 푀×푁 is a diagonal matrix with nonnegative entries being singular values of .
The columns of matrices  and  are called, respectively, the left-and right-singular vectors of .
Using the SVD of , the matrix  † is computed as  † = Σ †  푇 , where Σ † is obtained taking the reciprocal of nonzero diagonal elements of Σ and transposing of such a matrix.Numerically, by zero one means value smaller than, e.g.,  = ⋅max(, )⋅max(Σ), where  is a machine precision, whereas  and  are the dimensions of the matrix Σ.
One of the first practical algorithms was the method proposed in 1958 by Hestenes [25], which was based on the Givens rotations of the input matrix.Later on, two fundamental SVD algorithms were presented: Golub-Kahan (1965) [26] and its modification, Golub-Reinsch (1970) [27].They are two step algorithms based on the Householder reflections transforming the input matrix to an upper bidiagonal form at the first step and applying the modified QR algorithm to compute the SVD of this perturbed matrix subsequently.The case of very small singular values was covered by the Demmel-Kahan algorithm [13].Other algorithms are, for instance, the one-sided Jacobi orthogonalization [28] and the divide-and-conquer methods [29].Their idea is to modify the second step of SVD computation using other eigenvalue algorithms.The state-of-the-art review with more SVD algorithms tailored for specific matrices can be found in [30].

Computing Pseudoinverse by Tikhonov Regularization
Jhurani [7,8] suggested this method for the pseudoinverse computation as it can easily account for the sparsity of the input matrix .Using the Tikhonov regularized matrix [31] in the limit we obtain where  is the identity matrix of the size compatible with  푇  or  푇 .Assuming  as a small positive number one can compute an approximation of  † .Usually,  is assumed as the small multiple of the squared smallest nonzero singular value of .Nevertheless, for a good approximation of , singular values have to be computed, slowing down the whole routine.
The crucial assumption is the initial approximation of  0 .Usually, it is assumed as with a scalar  being bounded as follows: It was demonstrated in [19] that the approximate optimum order  is the integer  ≥ 2 minimizing where  and  are dimensions of matrix .
. .Iterative Methods Based on the Karush-Kuhn-Tucker Equations.Most of the numerical algorithms for the Moore-Penrose pseudoinverse are based on factorization.Despite their high reliability, they are computationally costly.Moreover, in the case of iterative methods, a wrong initial approximation deteriorates the convergence rate or even leads to the lack of convergence.
In [15] the authors proposed a novel iterative method.They formulated the pseudoinverse matrix computation problem as the matrix norm optimization problem with the Karush-Kuhn-Tucker conditions.Subsequently, the acceleration scheme was demonstrated.They proved in [15] that their method is globally convergent for arbitrary initial conditions.
Authors presented two versions of the algorithm for rectangular matrices  ≥  and  < , where  and  are dimensions of the input matrix .For the sake of brevity, only the first case is presented.
For  ≥  and given matrix  (0) ∈ R 푁×푀 as well as positive scalar  we compute iteratively ( = 0, 1, . . ., ) where  1 = ( 푇 ) 2 and 퐹 ≤  and  ≥ ) are satisfied, then the pseudoinverse matrix is computed as follows: The inverse matrix that has to be computed in the afore presented algorithm is evaluated only once because it does not contain any iterative term.

The qrginv Method
In [16] the authors proposed an efficient and elegant method (called ginv) for the pseudoinverse of matrix  ∈ R 푀×푁 computation.Their concept is based on a special type of the tensor product of two vectors.However, it was limited to full rank rectangular matrices and to square matrices with at least one zero row or column (with the rest of the matrix being full rank).
This method was extended to arbitrary singular matrices in [17] constituting the qrginv method.In order to obtain this,  factorization and the reverse order law for generalized inverse matrices were used.
Theorem 1.Let  be a Hilbert space.If  = ∑ 푁 푖=1  푖 ⊗  푖 is a rank-N operator, then its generalized inverse is also a rank-N operator and for each  ∈  it is defined by the relation where the unknown functions  푖 are the solutions of an appropriately defined  ×  linear system of equations.
The system of equations resulting from the properties of pseudoinverse matrices has a following form: The matrix of coefficients  ∈ × is called the Gram matrix ( 푖 푗 = ⟨ 푖 ,  푗 ⟩).For each  = 1, 2, . . .,  we look for  푖 ( 푗 ) in the following expression: where  푖 ( 푗 ) for  = 1, 2, , . . .,  constitutes the -th row of pseudoinverse matrix  † that has a form A generalization of the ginv method to arbitrary singular matrices was proposed in [17].The authors proved the following theorem.Using theorem 3.3.11from [33] they proposed the formula for  † using the properties of matrix .Theorem 3. Let  ∈ R 푁×푀 be a matrix with rank() =  > 0.
Since Â = , where  is a permutation matrix from the  decomposition, it can be written that  = .Finally, the authors observed using Theorem 2 and properties of  that where the matrix  † is computed using ginv method.Following [34], authors compute the rank of  11 as the number of columns, in which at least one coefficient with the absolute value greater than the set tolerance exists.As in [34], it is assumed to be equal to 10 −5 .

Numerical Tests
The local numerical homogenization routine uses pseudoinverses of the stiffness matrices.Thus, for our tests we have chosen several of such matrices.They are assembled in macroelement domains fine mesh matrices with no kinematic boundary conditions accounted for.We use the standard Bubnov-Galerkin approach to solve a given boundary value problem.First, we multiply the governing equation by an arbitrary test function V and integrate both sides over the analyzed domain.Then, we integrate by parts selected left hand side terms to decrease the derivative order of the searched function and to account for the prescribed boundary conditions.As a consequence, the derivative of the test function is introduced.We assume that both the solution and test function are approximated by the same basis functions.Solution to the obtained variational formulation over the domain discretized into a number of finite elements constitutes the set of degrees of freedom (coefficients) of the assumed basis functions.Typically, polynomials of various orders are used.In this manuscript, we understand the approximation order as the order of these functions.Stiffness matrices, that are the field of our interest, are the left hand side coefficient matrices.Due to the discretization and the local basis function supports, the global stiffness matrix is the assembly of the element-wise quantities.In order to solve the resulting system of linear algebraic equation, one needs to account for the kinematic boundary conditions.However, for our purposes, we operate on pure stiffness matrices without any modification.They are singular and positivedefinite.In the case of linear basis functions, these matrices are also banded.The bandwidth is the number of degrees of freedom of a single element.The rank of stiffness matrices is approximately equal (slightly smaller, obviously, due to their singularity) to the global number of degrees of freedom.The condition number depends on the approximation order, size of the element, and material parameters (e.g., for the most complex matrix analyzed in Test 5, it is equal to 6.4 ⋅10 13 ).Further details on the FEM and higher order elements can be found, e.g., in [18,35].
. .D Tests.For the 1D tests we analyze a simple example of a bar of length , the Young modulus , and cross-sectional area  for which the governing differential equation is the following: . . .Test .Two elements of equal lengths with linear shape functions were used for  = 2 ( = 1).The -th entry of the element stiffness matrix is where  is a given basis function.The assembled matrix is then with the pseudoinverse Results for this test are shown in Table 1.For every solution obtained using methods from Sections 3 to 6 we compare CPU time and verify if the conditions for the Moore-Penrose pseudoinverse (3÷6) hold.The values shown in the tables are the averages of 100 program runs.The SVD solution was obtained using the Matlab's pinv function that uses the algorithm from [13].In the case of the method presented in Section 5.2, only its accelerated version is considered in the tests.

Remarks
(i) For Matlab's pinv default  was used.
(ii) The parameter  in the Tikhonov regularization was set to 1 ⋅ 10 −9 , even though the squared minimal nonzero singular value equals 8.1079 ⋅ 10 −17 .For such a small number pseudoinverse cannot be found, since the matrices to be inverted in (7) are singular to working precision.Thus, the next larger singular value was assumed for  approximation.
(iii) In Ben-Israel's algorithm  = 3 and  = 1/9 for the initial matrix were used.The stopping criteria were used as follows: the maximum number of iterations cannot exceed 100 or (iv) In the accelerated scheme  0 = 1 ⋅ 10 −9 .The tolerance for the error measured as proposed in [15] was set to 1 ⋅ 10 −8 .
. . .Test .The next test is a repetition of the previous with higher approximation order,  1 = 2 for the first element and   2 = 3 for the second one.The Peano shape functions were used.The assembled stiffness matrix is as follows: and its exact pseudoinverse is A comparison of pseudoinverse algorithms is shown in Table 2.
. . .Test .In this test we analyze the same domain as previously but the material is heterogeneous.The Young moduli are equal to  1 = 5.1 ⋅ 10 2 and  2 = 5.1 ⋅ 10 3 .The finite element mesh consists of 10 elements with the Young modulus equal to  1 and  2 interchangeably.Randomly selected approximation orders (5, 8, 9, 2, 6, 5, 1, 4, 2, 8) for the elements were used.For the sake of clarity only results of pseudoinverse computations are shown in a tabular form (Table 3).The parameter  in the Tikhonov regularization was set equal to 1.3 ⋅ 10 −3 .
In the case of the 1D tests all of the tested algorithms provided reliable results in a short time as the matrices  were very small.Concluding 1D tests part, (i) the SVD was used with default Matlab's pinv tolerance with a good performance; (ii) the Tikhonov regularization is very sensitive to parameter .Looking for 's giving reasonable results was very time consuming.Evaluating their values, performed in terms of singular values, is very costly and rather impractical in the case of large matrices; (iii) the method proposed in [19] turned out to be very competitive; (iv) the KKT based method was always convergent but its efficiency very highly depends on the initial approximation of  0 ; (v) the method qrginv showed the best performance.Frequently, it was even better in terms of time and accuracy than Matlab's pinv.
Only the SVD [13], Ben-Israel's [19] higher order iterative method and the qrginv algorithm [17]  . . .Test .In this test, the Peano shape functions of the third order were used for the approximation at each direction.The matrix  has dimensions 1029 × 1029.We provide this matrix in the supplementary file matrix K.mat that was created using Matlab.Results of pseudoinverse computations are shown in Table 6.
Looking at Tables 5 and 6, it can be observed that the method proposed in [19] did not converge.The algorithm was stopped only because of the exceeding of the maximum iteration number.Being unconditionally convergent, the method turns out to be rather impractical in real applications.The pinv and qrginv functions turned out to be very fast.Unfortunately, homogenization based on their application yields unreliable results.The higher the approximation order was, the worse the final results were.We dealt with a super problem with the pseudoinversion embedded thus.Without this super problem, one could consider the results shown in Table 5 (for the SVD and qrginv) as the reliable ones.Poor results of the trivial local numerical homogenization example (uniaxial test for a homogeneous material) led us to a thorough analysis.It ended with an interesting observation concerning the pseudoinverse matrix computation described below.
Let us denote pseudoinverse matrices obtained using pinv and qrginv functions, respectively,  † 1 and  † 2 .Consequently, their pseudoinverses are  1 and  2 .One may compute the following error indicators: Let us remind that in the 4b and 5 tests the approximation orders equal to 1 and 3, respectively.The above presented Mathematical Problems in Engineering

Figure 1 :
Figure 1: Scheme of the local numerical homogenization routine.
Upto the machine precision.
Upto the machine precision.

Table 1 :
CPU time and relative errors, Test 1.

Table 5 .
are considered in the next chapter for the 3D tests...D Tests.For the 3D tests, we analyze a homogeneous cube of dimensions 1 × 1 × 1.The whole domain was discretized with 8 equal finite elements.Linear elasticity case was considered.The Young modulus  = 1 ⋅ 10 6 and the Poisson ratio ] = 0.3....Test a.In this test, linear shape functions were used.The matrix  has dimensions 81 × 81.Results of Moore-Penrose pseudoinverse computations are shown in Table4.. ..Test b. Ithe previous test, the SVD-based algorithm gave the worst results.It was due to the default  used for pinv Matlab's function.This parameter determines how small singular value is treated as zero.In order to obtain results with higher accuracy we proposed Algorithm 1:Updated results for the SVD-based algorithm are shown in

Table 2 :
CPU time and relative errors, Test 2.

Table 3 :
CPU time and relative errors, Test 3.