Abstract

Computing the pseudoinverse of a matrix is an essential component of many computational methods. 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.

1. 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 [68] 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 [1012]. 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, 1417]. 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.

2. 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 3Dhp 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., ).

Given symmetric fine mesh element stiffness matrices assembled into , a nonzero fine mesh load vector , an interpolation matrix , a positive-definite symmetric weight matrix , and a dimensionless small parameter find a symmetric matrix for which reaches minimum, where

and 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 that minimizes (1) is the pseudoinverse of the local effective stiffness matrix . Derivation of (1) is presented in [6, 8]. Imposing the symmetry conditions on ) 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 a pseudoinverse of matrix , 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 ;(ii)for the pseudoinverse of the effective stiffness matrix (, 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 (). 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].

3. 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 can be decomposed into , where and are orthogonal matrices and 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., , 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].

4. 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 obtainwhere 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.

5. Computing Pseudoinverse by Iterative Methods

5.1. Higher Order Iterative Methods

In [19] Ben-Israel presents a group of methods leading to the pseudoinverse computation. For a nonzero matrix with the initial approximation and its residual satisfying ) and , respectively, for any integer the sequenceconverges to with .

The crucial assumption is the initial approximation of . Usually, it is assumed aswith a scalar being bounded as follows:

It was demonstrated in [19] that the approximate optimum order is the integer minimizingwhere and are dimensions of matrix .

5.2. 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 a given matrix as well as a positive scalar we compute iteratively , where is the maximum iteration number, for where , is the identity matrix of dimensions , , and . Then we compute aswith and . If stopping criteria (i.e., and ) are satisfied, then the pseudoinverse matrix is computed as follows:

Subsequently, the authors eliminated from the afore described algorithm a time step , which is matrix size dependent and thus convergence deteriorating in the case of large matrices. The obtained accelerated algorithm is as follows.

For and given matrix as well as positive scalar we compute iteratively ()where and . If stopping criteria (i.e., 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.

6. The qrginv Method

In [16] the authors proposed an efficient and elegant method (called ginv) for the pseudoinverse of matrix 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.

Authors in [16] consider for any () the mappingwith and being the sets of orthonormal and linearly independent vectors of . denotes the standard inner product in . Every rank- operator can be expressed asAn operator is linear and its matrix representation can be written in the following form: with being first elements of the standard basis and constituting -th column of the matrix , .

Recalling Theorem 3.2 from [32] the authors proposed the procedure for the Moore-Penrose inverse of the corresponding tensor-product matrix.

Theorem 1. Let be a Hilbert space. If 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 we look for in the following expression:where for 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.

Theorem 2. Let be the factorization of . Then, .

Using theorem 3.3.11 from [33] they proposed the formula for using the properties of matrix .

Theorem 3. Let be a matrix with rank(. Then, there exist matrices , , such that is obtained by by permuting its columns, is orthogonal, , is nonsingular and upper triangular, and .

Since , where is a permutation matrix from the decomposition, it can be written that . Finally, the authors observed using Theorem 2 and properties of thatwhere the matrix is computed using ginv method. Following [34], authors compute the rank of 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 .

7. 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 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 positive-definite. 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 ). Further details on the FEM and higher order elements can be found, e.g., in [18, 35].

7.1. 1D 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:

7.1.1. Test 1

Two elements of equal lengths with linear shape functions were used for ). The -th entry of the element stiffness matrix is , where is a given basis function. The assembled matrix is thenwith 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 , even though the squared minimal nonzero singular value equals . 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 and for the initial matrix were used. The stopping criteria were used as follows: the maximum number of iterations cannot exceed 100 or is less than .(iv)In the accelerated scheme . The tolerance for the error measured as proposed in [15] was set to .

7.1.2. Test 2

The next test is a repetition of the previous with higher approximation order, for the first element and 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.

7.1.3. Test 3

In this test we analyze the same domain as previously but the material is heterogeneous. The Young moduli are equal to and . The finite element mesh consists of 10 elements with the Young modulus equal to and interchangeably. Randomly selected approximation orders , 8, 9, 2, 6, 5, 1, 4, 2, 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 .

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 ;(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] are considered in the next chapter for the 3D tests.

7.2. 3D Tests

For the 3D tests, we analyze a homogeneous cube of dimensions . The whole domain was discretized with 8 equal finite elements. Linear elasticity case was considered. The Young modulus and the Poisson ratio .

7.2.1. Test 4a

In this test, linear shape functions were used. The matrix has dimensions . Results of Moore-Penrose pseudoinverse computations are shown in Table 4.

7.2.2. Test 4b

In the 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:

Set ,
while    do
if    then
and
end if
=pinv(
if    then
end if
if    then
end if
if    then
end if
if    then
end if
end while
return  

Updated results for the SVD-based algorithm are shown in Table 5.

7.2.3. Test 5

In this test, the Peano shape functions of the third order were used for the approximation at each direction. The matrix has dimensions . 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, and . Consequently, their pseudoinverses are and . One may compute the following error indicators:(i)In Test 4b: but , when .(ii)In Test 5: , and .

Let us remind that in the 4b and 5 tests the approximation orders equal to 1 and 3, respectively. The above presented numbers practically disable the application of the local numerical homogenization with Peano higher order shape functions. As it was mentioned at the beginning of this paper, one has to compute the pseudoinverse twice. First, the input matrix is the assembled stiffness matrix. Even if it is pseudoinverted accurately, then (after some further calculations) one has to pseudoinvert the pseudoinverse of the effective macroscale stiffness matrix. So far, we have not found the algorithm giving for matrices of our interest (using Peano higher order shape functions) . In Test 4b, matrices and have almost the same norms, however, only the error indicator computed for is satisfactory. In Test 5, all three error indicators were prohibitively large. Only the results of the local numerical homogenization using a linear approximation and the pinv function were correct.

In order to proceed with the local numerical homogenization analysis, we applied the integrated Legendre or Lobatto polynomials (see, e.g., [35]) of higher order. Such basis functions and, consequently, the resulting stiffness matrices exhibit some superior quantities (c.f. [35]) comparing to the Peano polynomials. They enabled us to obtain afore presented error indicators () of order for the pinv function. Thus, only the Demmel-Kahan method is reliable for pseudoinversion of stiffness matrices. Results for Test 5 with Legendre and Lobatto shape functions are shown in Tables 7 and 8, respectively.

8. Concluding Remarks

We proposed a very demanding test for the computation of the Moore-Penrose pseudoinverse matrices. The finite element stiffness matrices resulting from higher order approximation are such difficult data for all algorithms of the pseudoinverse computation. It could be concluded that, for elasticity stiffness matrices, the Demmel-Kahan algorithm was the only suitable method for our purposes. The proposed matrices can also serve, however, as benchmarks for new algorithms for the Moore-Penrose pseudoinverse.

The difficulty of the stiffness pseudoinverse computation for the higher order approximation can be explained by the increasing condition number (the ratio of the largest and the smallest nonzero singular values) of these matrices. Stiffness matrices obtained using Peano shape functions exhibit higher condition numbers comparing to the integrated Legendre or Lobatto polynomials. Therefore, only the application of those functions allows one to obtain for arbitrary higher order approximation with a very high accuracy. It was obtained, however, only using the pinv function implementing algorithm presented in [13]. The remaining methods did not pass the test.

Data Availability

The authors have attached to this submission the file “matrix_K”, which stores the stiffness matrix used in Test 5.

Conflicts of Interest

No conflicts of interest can be stated.

Supplementary Materials

The Supplementary Material file stores the global stiffness matrix resulting from the linear elasticity problem. The analyzed domain (a unit cube) was discretized with 8 equal elements. Material data can be found in “Numerical Tests”. The approximation order is equal to 3 and the Peano shape functions were used. The stiffness matrix has dimensions 1029x1029. It was used for the pseudoinversion tests. (Supplementary Materials)