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

Authors: Marek Klimczak; Witold Cecot
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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.


---

## Body

## 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 [6–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–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–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.


## 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 (Figure1). 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 the3Dhp 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.Figure 1
Scheme of the local numerical homogenization routine.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 withhat symbol (e.g., x^).Given symmetric fine mesh element stiffness matrices assembled intoK∈RN×N, a nonzero fine mesh load vector f∈RN, an interpolation matrix A∈RN×M, a positive-definite symmetric weight matrix B∈RN×N, and a dimensionless small parameter ϵ>0 find a symmetric matrix K^†∈RM×M for which E reaches minimum, where(1)EK^†=12K†-AK^†ATfB2+ϵ2K†-AK^†ATF,B2f22andx2=trace(xTx) denotes the Euclidean norm, xB=trace(xTBx) denotes the Euclidean norm weighted with B, and XF,B=trace(XTBX) denotes the Frobenius norm weighted with B. x and X are arbitrary column vector and matrix, respectively. The numbers of degrees of freedom for the fine and coarse meshes are denoted with N and M.The expression forE 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^†)T=0) into the Lagrangian, it was proved that problem (1) is equivalent to the continuous-time Lyapunov equation in a following form:(2)CX+XCT=Dwith auxiliary matricesC≔V-1U,D≔V-1(W+WT)V-1,U≔ATBA,V≔AT(ffT+ϵf22I)A,W≔ATBK†ffT+ϵf22IA.In (1) X† denotes the Moore-Penrose pseudoinverse of matrix X [20, 21]. In the case of the real matrices we call a matrix X†∈RN×M a pseudoinverse of matrix X∈RM×N, if it fulfills the following conditions:(3)XX†X=X(4)X†XX†=X†(5)XX†T=XX†(6)X†XT=X†XA reliable computation ofX† 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 matrixK∈RN×N;(ii)
for the pseudoinverse of the effective stiffness matrix (K^†∈RM×M, where M≪N).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 (A∈RM×N). 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 matrixA∈RM×N can be decomposed into A=UΣVT, where U∈RM×M and V∈RN×N are orthogonal matrices and Σ∈RM×N is a diagonal matrix with nonnegative entries being singular values of A. The columns of matrices U and V are called, respectively, the left- and right-singular vectors of A.Using the SVD ofA, the matrix A† is computed as A†=VΣ†UT, 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., TOL=ϵ·max(M,N)·max(Σ), where ϵ is a machine precision, whereas M and N 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 A. Using the Tikhonov regularized matrix [31] in the limit we obtain(7)A†=limδ→0ATA+δI-1AT=limδ→0ATAAT+δI-1where I is the identity matrix of the size compatible with ATA or AAT. Assuming δ as a small positive number one can compute an approximation of A†. Usually, δ is assumed as the small multiple of the squared smallest nonzero singular value of A. 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 A∈RM×N with the initial approximation X0 and its residual R0 satisfying X0=ATBAT(B∈RM×N) and ρ(R0)<1, respectively, for any integer p≥2 the sequence(8)Xk+1=XkI+I-AXk+I-AXk2+⋯+I-AXkp-1converges to A† with k→∞.The crucial assumption is the initial approximation ofX0. Usually, it is assumed as(9)X0=βAT,with a scalar β being bounded as follows:(10)0<β<2λmaxAAT.It was demonstrated in [19] that the approximate optimum order p is the integer p≥2 minimizing(11)fp=p-2+2N/Mlnp,p≠2q,q=1,2,…2q-1+N/Mqln2,p=2q,q=1,2,…where M and N are dimensions of matrix A.


### 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 matricesM≥N and M<N, where M and N are dimensions of the input matrix A. For the sake of brevity, only the first case is presented.ForM≥N and a given matrix Y(0)∈RN×M as well as a positive scalar ϵ we compute iteratively k=0,1,…,L, where L is the maximum iteration number, for i=1,…,M(12)hik=aiT-B1yik2IN+B1aiT-B1yik2,where B1=(ATA)2, IN is the identity matrix of dimensions N×N, AT=[a1,…,aM], and Y(k)=[y1(k),…,yM(k)]. Then we compute Y(k+1) as(13)Yk+1=Yk+WNAT-B1YkHk,with H(k)=diag(h1(k),…,hM(k)) and WN=(IN+B1). If stopping criteria (i.e., AT-B1Y(k+1)F2≤ϵ and k≥L) are satisfied, then the pseudoinverse matrix is computed as follows:(14)A†=ATAYk+1Subsequently, the authors eliminated from the afore described algorithm a time stephi(k), which is matrix A size dependent and thus convergence deteriorating in the case of large matrices. The obtained accelerated algorithm is as follows.ForM≥N and given matrix Y(0)∈RN×M as well as positive scalar ϵ we compute iteratively (k=0,1,…,L)(15)Yk+1=Yk+WNAT-B1Yk,where B1=(ATA)2 and WN=(IN+(ATA)2)-1. If stopping criteria (i.e., AT-B1Y(k+1)F2≤ϵ and k≥L) are satisfied, then the pseudoinverse matrix is computed as follows:(16)A†=ATAYk+1The inverse matrix that has to be computed in the afore presented algorithm is evaluated only once because it does not contain any iterative term.


## 5.1. Higher Order Iterative Methods

In [19] Ben-Israel presents a group of methods leading to the pseudoinverse computation. For a nonzero matrix A∈RM×N with the initial approximation X0 and its residual R0 satisfying X0=ATBAT(B∈RM×N) and ρ(R0)<1, respectively, for any integer p≥2 the sequence(8)Xk+1=XkI+I-AXk+I-AXk2+⋯+I-AXkp-1converges to A† with k→∞.The crucial assumption is the initial approximation ofX0. Usually, it is assumed as(9)X0=βAT,with a scalar β being bounded as follows:(10)0<β<2λmaxAAT.It was demonstrated in [19] that the approximate optimum order p is the integer p≥2 minimizing(11)fp=p-2+2N/Mlnp,p≠2q,q=1,2,…2q-1+N/Mqln2,p=2q,q=1,2,…where M and N are dimensions of matrix A.


## 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 matricesM≥N and M<N, where M and N are dimensions of the input matrix A. For the sake of brevity, only the first case is presented.ForM≥N and a given matrix Y(0)∈RN×M as well as a positive scalar ϵ we compute iteratively k=0,1,…,L, where L is the maximum iteration number, for i=1,…,M(12)hik=aiT-B1yik2IN+B1aiT-B1yik2,where B1=(ATA)2, IN is the identity matrix of dimensions N×N, AT=[a1,…,aM], and Y(k)=[y1(k),…,yM(k)]. Then we compute Y(k+1) as(13)Yk+1=Yk+WNAT-B1YkHk,with H(k)=diag(h1(k),…,hM(k)) and WN=(IN+B1). If stopping criteria (i.e., AT-B1Y(k+1)F2≤ϵ and k≥L) are satisfied, then the pseudoinverse matrix is computed as follows:(14)A†=ATAYk+1Subsequently, the authors eliminated from the afore described algorithm a time stephi(k), which is matrix A size dependent and thus convergence deteriorating in the case of large matrices. The obtained accelerated algorithm is as follows.ForM≥N and given matrix Y(0)∈RN×M as well as positive scalar ϵ we compute iteratively (k=0,1,…,L)(15)Yk+1=Yk+WNAT-B1Yk,where B1=(ATA)2 and WN=(IN+(ATA)2)-1. If stopping criteria (i.e., AT-B1Y(k+1)F2≤ϵ and k≥L) are satisfied, then the pseudoinverse matrix is computed as follows:(16)A†=ATAYk+1The 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. Theqrginv Method

In [16] the authors proposed an efficient and elegant method (called ginv) for the pseudoinverse of matrix A∈RM×N 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, QR factorization and the reverse order law for generalized inverse matrices were used.Authors in [16] consider for any x∈RK (N<K) the mapping(17)e⊗f:RK→RKwithe⊗fx=x,efwith {e1,…,eN} and {f1,…,fN} being the sets of orthonormal and linearly independent vectors of RK. ·,· denotes the standard inner product in RN. Every rank-N operator A can be expressed as(18)A=∑i=1Nei⊗fi.An operator A is linear and its matrix representation can be written in the following form: (19)A=f11f21…fN10…0f12f22…fN20…0⋮⋮…⋮⋮⋮f1Kf2K…fNK0…0.with {e1,…,eN} being first N elements of the standard RK basis and fi=(fi1,fi2,…,fiK) constituting i-th column of the matrix A, i=1,2,…,N.Recalling Theorem 3.2 from [32] the authors proposed the procedure for the Moore-Penrose inverse of the corresponding tensor-product matrix.Theorem 1.
LetH be a Hilbert space. If A=∑i=1Nei⊗fi is a rank-N operator, then its generalized inverse is also a rank-N operator and for each x∈H it is defined by the relation (20)A†x=∑i=1Nλixeiwhere the unknown functions λi are the solutions of an appropriately defined N×N linear system of equations.The system of equations resulting from the properties of pseudoinverse matrices has a following form:(21)x,fi=∑j=1Nλjxfi,fj,i=1,2,…,NThe matrix of coefficients G∈N×N is called the Gram matrix (gij=fi,fj). For each j=1,2,…,N we look for λi(ej) in the following expression:(22)A†ej=∑i=1Nλiejei,where λi(ej) for i=1,2,,…,N constitutes the i-th row of pseudoinverse matrix A† that has a form (23)A†=λ1e1λ1e2…λ1eKλ2e1λ2e2…λ2eK⋮⋮⋮⋮λNe1λNe2…λNeK00…0⋮⋮⋮⋮00…0.A generalization of theginv method to arbitrary singular matrices was proposed in [17]. The authors proved the following theorem.Theorem 2.
LetA=QR be the QR factorization of A. Then, A†=R†Q∗.Using theorem 3.3.11 from [33] they proposed the formula for A† using the properties of matrix R.Theorem 3.
LetA∈RN×M be a matrix with rank(A)=r>0. Then, there exist matrices A^, Q, R such that A^ is obtained by A by permuting its columns, Q∈RN×N is orthogonal, R=R11R1200∈RN×M, R11∈Rr×r is nonsingular and upper triangular, and A^=QR.SinceA^=AP, where P is a permutation matrix from the QR decomposition, it can be written that AP=QR. Finally, the authors observed using Theorem 2 and properties of P that(24)A†=PR†QT,where the matrix R† is computed using ginv method. Following [34], authors compute the rank of R11 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.


## 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 functionv 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 ·1013). 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 lengthL, the Young modulus E, and cross-sectional area A for which the governing differential equation is the following:(25)-AEu′′=q,x∈0,L


#### 7.1.1. Test 1

Two elements of equal lengths with linear shape functions were used forL=2(AE=1). The ij-th entry of the element stiffness matrix is kijel=∫elφi′AEφj′dx, where φ is a given basis function. The assembled matrix is then(26)K=1-10-12-10-11with the pseudoinverse (27)K†=59-19-49-1929-19-49-1959Results for this test are shown in Table1. 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.Table 1
CPU time and relative errors, Test 1.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0009 0.1994 ⋅ 1 0 - 15 0.4083 ⋅ 1 0 - 15 0.2289 ⋅ 1 0 - 15 0.3284⋅10-15 [8] 0.0005 0.1131⋅10-06 0.0377⋅10-06 0.1460⋅10-06 0 ∗ [19] 0.0012 0.0142⋅10-14 0.1776⋅10-14 0.0324⋅10-14 0.0096 ⋅ 1 0 - 14 [15] 0.0008 0.1930⋅10-04 0.5790⋅10-04 0.8899⋅10-15 0.5382 ⋅ 1 0 - 15 [17] 0.0007 0.2538⋅10-15 0.3187⋅10-15 0.2289⋅10-15 0.7065 ⋅ 1 0 - 15 U p   ∗ to the machine precision.Remarks(i)
For Matlab’spinv default TOL 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 algorithmp=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 Xk+1-XkF/XkF is less than 1·10-8. (iv)
In the accelerated schemeY0=1·10-9K. The tolerance for the error measured as proposed in [15] was set to 1·10-8.


#### 7.1.2. Test 2

The next test is a repetition of the previous with higher approximation order,p1=2 for the first element and p2=3 for the second one. The Peano shape functions were used. The assembled stiffness matrix is as follows: (28)K=1-10000-120-10000130000-1010000001300000015and its exact pseudoinverse is (29)K†=59-190-4900-19290-1900003000-49-1905900000030000005.A comparison of pseudoinverse algorithms is shown in Table 2.Table 2
CPU time and relative errors, Test 2.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0014 0.1968 ⋅ 1 0 - 15 0.0648 ⋅ 1 0 - 15 0.1448 ⋅ 1 0 - 15 0.2077⋅10-15 [8] 0.0010 0.1786⋅10-05 0.0096⋅10-05 0.1477⋅10-05 0 ∗ [19] 0.0022 0.0011⋅10-13 0.1189⋅10-13 0.0016⋅10-13 0.0045 ⋅ 1 0 - 13 [15] 0.0979 0.0312⋅10-03 0.3759⋅10-03 0.0311⋅10-13 0.0348 ⋅ 1 0 - 13 [17] 0.0004 0.2505⋅10-15 0.5060⋅10-15 0.1448⋅10-15 0.4469 ⋅ 1 0 - 15 U p   ∗ to the machine precision.


#### 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 toE1=5.1·102 and E2=5.1·103. The finite element mesh consists of 10 elements with the Young modulus equal to E1 and E2 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.Table 3
CPU time and relative errors, Test 3.
Method Time [s] K K † K - K F K F K † K K † - K † F K † F K K † T - K K † F K K † F K † K T - K † K F K † K F SVD 0.0059 0.0005 ⋅ 1 0 - 11 0.0434 ⋅ 1 0 - 11 0.0644 ⋅ 1 0 - 11 0.2611 ⋅ 1 0 - 11 [8] 0.0025 0.6177 ⋅ 1 0 - 06 0.0095 ⋅ 1 0 - 00 0.3379 ⋅ 1 0 - 06 0.8813 ⋅ 1 0 - 07 [19] 0.0170 0.0001 ⋅ 1 0 - 10 0.3446 ⋅ 1 0 - 10 0.0002 ⋅ 1 0 - 10 0.0150 ⋅ 1 0 - 10 [15] 3.4590 0.2453 ⋅ 1 0 - 06 0.1265 ⋅ 1 0 - 03 0.3349 ⋅ 1 0 - 05 0.0016 ⋅ 1 0 - 00 [17] 0.0026 0.0030 ⋅ 1 0 - 12 0.0538 ⋅ 1 0 - 12 0.0429 ⋅ 1 0 - 12 0.1366 ⋅ 1 0 - 12In the case of the 1D tests all of the tested algorithms provided reliable results in a short time as the matricesK were very small. Concluding 1D tests part,(i)
the SVD was used with default Matlab’spinv 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 ofY0;(v)
the methodqrginv 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 dimensions1×1×1. The whole domain was discretized with 8 equal finite elements. Linear elasticity case was considered. The Young modulus E=1·106 and the Poisson ratio ν=0.3.


#### 7.2.1. Test 4a

In this test, linear shape functions were used. The matrixK has dimensions 81×81. Results of Moore-Penrose pseudoinverse computations are shown in Table 4.Table 4
CPU time and relative errors, Test 4a.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0078 0.4808 ⋅ 1 0 - 03 0.2340 ⋅ 1 0 - 02 0.4283 ⋅ 1 0 - 02 0.4977⋅10-02 [19] 0.0280 0.7028⋅10-13 0.3091⋅10-09 0.1318⋅10-14 0.1547 ⋅ 1 0 - 14 [17] 0.0069 0.7028⋅10-13 0.1572⋅10-14 0.1809⋅10-11 0.5138 ⋅ 1 0 - 14


#### 7.2.2. Test 4b

In the previous test, the SVD-based algorithm gave the worst results. It was due to the defaultTOL 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:Algorithm 1:Given K compute its pseudoinverse K† using pinv function.
SetTOLpinv=1, TOL=1⋅10-6 while  icond<4  do i c o n d = 0 T O L p i n v = 0.1 ⋅ T O L p i n v if  TOLpinv<1⋅10-20  then T O L p i n v = 1 ⋅ 1 0 - 1 and T O L = 10 ⋅ T O L end if K †=pinv(K,TOLpinv if  KK†K-KF/KF<TOL  then i c o n d = i c o n d + 1 end if if  K†KK†-K†F/K†F<TOL  then i c o n d = i c o n d + 1 end if if  KK†T-KK†F/KK†F<TOL  then i c o n d = i c o n d + 1 end if if  (K†K)T-K†KF/K†KF<TOL  then i c o n d = i c o n d + 1 end if end while return  K†Updated results for the SVD-based algorithm are shown in Table5.Table 5
CPU time and relative errors, Test 4b.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0064 0.7039 ⋅ 1 0 - 13 0.2704 ⋅ 1 0 - 14 0.6579 ⋅ 1 0 - 14 0.5456⋅10-14


#### 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 matrixK 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.Table 6
CPU time and relative errors, Test 5 (Peano shape functions).
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 1.5543 0.7596 ⋅ 1 0 - 07 0.1946 ⋅ 1 0 - 10 0.6767 ⋅ 1 0 - 09 0.2657⋅10-10 [19] 175.44 0.0022⋅10+00 0.0018⋅10+00 0.0008⋅10+00 1.4142 ⋅ 1 0 + 00 [17] 1.4129 0.0001⋅10-09 0.0004⋅10-09 0.5837⋅10-09 0.1063 ⋅ 1 0 - 09Looking at Tables5 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 usingpinv and qrginv functions, respectively, K1† and K2†. Consequently, their pseudoinverses are K1 and K2. One may compute the following error indicators:(i)
In Test 4b:K-K1F/KF=1.0434·10-13 but K-K2F/KF=1, when K1-K2F/K1F=5.2436·10-12.(ii)
In Test 5:K-K1F/KF=1, K-K2F/KF=1 and K1-K2F/K1F=12.7059(!).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)(K†)†=K. In Test 4b, matrices K1 and K2 have almost the same norms, however, only the error indicator K-KcomputedF/KF computed for K1 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 (K-KcomputedF/KF) of order 1·10-12 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.Table 7
CPU time and relative errors, Test 5 (Legendre shape functions).
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 1.4112 0.0047 ⋅ 1 0 - 11 0.0344 ⋅ 1 0 - 11 0.1863 ⋅ 1 0 - 11 0.0671⋅10-11 [19] 101.35 0.0017⋅10+00 0.0018⋅10+00 0.0009⋅10+00 1.4142 ⋅ 1 0 + 00 [17] 1.8285 0.0005⋅10-10 0.0005⋅10-10 0.2127⋅10-10 0.0054 ⋅ 1 0 - 10Table 8
CPU time and relative errors, Test 5 (Lobatto shape functions).
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 1.4259 0.0357 ⋅ 1 0 - 12 0.1948 ⋅ 1 0 - 12 0.5772 ⋅ 1 0 - 12 0.2635⋅10-12 [19] 107.22 0.0013⋅10+00 0.0017⋅10+00 0.0009⋅10+00 1.4142 ⋅ 1 0 + 00 [17] 1.8454 0.0004⋅10-10 0.0005⋅10-10 0.1100⋅10-10 0.0020 ⋅ 1 0 - 10


## 7.1. 1D Tests

For the 1D tests we analyze a simple example of a bar of lengthL, the Young modulus E, and cross-sectional area A for which the governing differential equation is the following:(25)-AEu′′=q,x∈0,L


### 7.1.1. Test 1

Two elements of equal lengths with linear shape functions were used forL=2(AE=1). The ij-th entry of the element stiffness matrix is kijel=∫elφi′AEφj′dx, where φ is a given basis function. The assembled matrix is then(26)K=1-10-12-10-11with the pseudoinverse (27)K†=59-19-49-1929-19-49-1959Results for this test are shown in Table1. 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.Table 1
CPU time and relative errors, Test 1.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0009 0.1994 ⋅ 1 0 - 15 0.4083 ⋅ 1 0 - 15 0.2289 ⋅ 1 0 - 15 0.3284⋅10-15 [8] 0.0005 0.1131⋅10-06 0.0377⋅10-06 0.1460⋅10-06 0 ∗ [19] 0.0012 0.0142⋅10-14 0.1776⋅10-14 0.0324⋅10-14 0.0096 ⋅ 1 0 - 14 [15] 0.0008 0.1930⋅10-04 0.5790⋅10-04 0.8899⋅10-15 0.5382 ⋅ 1 0 - 15 [17] 0.0007 0.2538⋅10-15 0.3187⋅10-15 0.2289⋅10-15 0.7065 ⋅ 1 0 - 15 U p   ∗ to the machine precision.Remarks(i)
For Matlab’spinv default TOL 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 algorithmp=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 Xk+1-XkF/XkF is less than 1·10-8. (iv)
In the accelerated schemeY0=1·10-9K. The tolerance for the error measured as proposed in [15] was set to 1·10-8.


### 7.1.2. Test 2

The next test is a repetition of the previous with higher approximation order,p1=2 for the first element and p2=3 for the second one. The Peano shape functions were used. The assembled stiffness matrix is as follows: (28)K=1-10000-120-10000130000-1010000001300000015and its exact pseudoinverse is (29)K†=59-190-4900-19290-1900003000-49-1905900000030000005.A comparison of pseudoinverse algorithms is shown in Table 2.Table 2
CPU time and relative errors, Test 2.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0014 0.1968 ⋅ 1 0 - 15 0.0648 ⋅ 1 0 - 15 0.1448 ⋅ 1 0 - 15 0.2077⋅10-15 [8] 0.0010 0.1786⋅10-05 0.0096⋅10-05 0.1477⋅10-05 0 ∗ [19] 0.0022 0.0011⋅10-13 0.1189⋅10-13 0.0016⋅10-13 0.0045 ⋅ 1 0 - 13 [15] 0.0979 0.0312⋅10-03 0.3759⋅10-03 0.0311⋅10-13 0.0348 ⋅ 1 0 - 13 [17] 0.0004 0.2505⋅10-15 0.5060⋅10-15 0.1448⋅10-15 0.4469 ⋅ 1 0 - 15 U p   ∗ to the machine precision.


### 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 toE1=5.1·102 and E2=5.1·103. The finite element mesh consists of 10 elements with the Young modulus equal to E1 and E2 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.Table 3
CPU time and relative errors, Test 3.
Method Time [s] K K † K - K F K F K † K K † - K † F K † F K K † T - K K † F K K † F K † K T - K † K F K † K F SVD 0.0059 0.0005 ⋅ 1 0 - 11 0.0434 ⋅ 1 0 - 11 0.0644 ⋅ 1 0 - 11 0.2611 ⋅ 1 0 - 11 [8] 0.0025 0.6177 ⋅ 1 0 - 06 0.0095 ⋅ 1 0 - 00 0.3379 ⋅ 1 0 - 06 0.8813 ⋅ 1 0 - 07 [19] 0.0170 0.0001 ⋅ 1 0 - 10 0.3446 ⋅ 1 0 - 10 0.0002 ⋅ 1 0 - 10 0.0150 ⋅ 1 0 - 10 [15] 3.4590 0.2453 ⋅ 1 0 - 06 0.1265 ⋅ 1 0 - 03 0.3349 ⋅ 1 0 - 05 0.0016 ⋅ 1 0 - 00 [17] 0.0026 0.0030 ⋅ 1 0 - 12 0.0538 ⋅ 1 0 - 12 0.0429 ⋅ 1 0 - 12 0.1366 ⋅ 1 0 - 12In the case of the 1D tests all of the tested algorithms provided reliable results in a short time as the matricesK were very small. Concluding 1D tests part,(i)
the SVD was used with default Matlab’spinv 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 ofY0;(v)
the methodqrginv 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.1.1. Test 1

Two elements of equal lengths with linear shape functions were used forL=2(AE=1). The ij-th entry of the element stiffness matrix is kijel=∫elφi′AEφj′dx, where φ is a given basis function. The assembled matrix is then(26)K=1-10-12-10-11with the pseudoinverse (27)K†=59-19-49-1929-19-49-1959Results for this test are shown in Table1. 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.Table 1
CPU time and relative errors, Test 1.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0009 0.1994 ⋅ 1 0 - 15 0.4083 ⋅ 1 0 - 15 0.2289 ⋅ 1 0 - 15 0.3284⋅10-15 [8] 0.0005 0.1131⋅10-06 0.0377⋅10-06 0.1460⋅10-06 0 ∗ [19] 0.0012 0.0142⋅10-14 0.1776⋅10-14 0.0324⋅10-14 0.0096 ⋅ 1 0 - 14 [15] 0.0008 0.1930⋅10-04 0.5790⋅10-04 0.8899⋅10-15 0.5382 ⋅ 1 0 - 15 [17] 0.0007 0.2538⋅10-15 0.3187⋅10-15 0.2289⋅10-15 0.7065 ⋅ 1 0 - 15 U p   ∗ to the machine precision.Remarks(i)
For Matlab’spinv default TOL 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 algorithmp=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 Xk+1-XkF/XkF is less than 1·10-8. (iv)
In the accelerated schemeY0=1·10-9K. The tolerance for the error measured as proposed in [15] was set to 1·10-8.


## 7.1.2. Test 2

The next test is a repetition of the previous with higher approximation order,p1=2 for the first element and p2=3 for the second one. The Peano shape functions were used. The assembled stiffness matrix is as follows: (28)K=1-10000-120-10000130000-1010000001300000015and its exact pseudoinverse is (29)K†=59-190-4900-19290-1900003000-49-1905900000030000005.A comparison of pseudoinverse algorithms is shown in Table 2.Table 2
CPU time and relative errors, Test 2.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0014 0.1968 ⋅ 1 0 - 15 0.0648 ⋅ 1 0 - 15 0.1448 ⋅ 1 0 - 15 0.2077⋅10-15 [8] 0.0010 0.1786⋅10-05 0.0096⋅10-05 0.1477⋅10-05 0 ∗ [19] 0.0022 0.0011⋅10-13 0.1189⋅10-13 0.0016⋅10-13 0.0045 ⋅ 1 0 - 13 [15] 0.0979 0.0312⋅10-03 0.3759⋅10-03 0.0311⋅10-13 0.0348 ⋅ 1 0 - 13 [17] 0.0004 0.2505⋅10-15 0.5060⋅10-15 0.1448⋅10-15 0.4469 ⋅ 1 0 - 15 U p   ∗ to the machine precision.


## 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 toE1=5.1·102 and E2=5.1·103. The finite element mesh consists of 10 elements with the Young modulus equal to E1 and E2 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.Table 3
CPU time and relative errors, Test 3.
Method Time [s] K K † K - K F K F K † K K † - K † F K † F K K † T - K K † F K K † F K † K T - K † K F K † K F SVD 0.0059 0.0005 ⋅ 1 0 - 11 0.0434 ⋅ 1 0 - 11 0.0644 ⋅ 1 0 - 11 0.2611 ⋅ 1 0 - 11 [8] 0.0025 0.6177 ⋅ 1 0 - 06 0.0095 ⋅ 1 0 - 00 0.3379 ⋅ 1 0 - 06 0.8813 ⋅ 1 0 - 07 [19] 0.0170 0.0001 ⋅ 1 0 - 10 0.3446 ⋅ 1 0 - 10 0.0002 ⋅ 1 0 - 10 0.0150 ⋅ 1 0 - 10 [15] 3.4590 0.2453 ⋅ 1 0 - 06 0.1265 ⋅ 1 0 - 03 0.3349 ⋅ 1 0 - 05 0.0016 ⋅ 1 0 - 00 [17] 0.0026 0.0030 ⋅ 1 0 - 12 0.0538 ⋅ 1 0 - 12 0.0429 ⋅ 1 0 - 12 0.1366 ⋅ 1 0 - 12In the case of the 1D tests all of the tested algorithms provided reliable results in a short time as the matricesK were very small. Concluding 1D tests part,(i)
the SVD was used with default Matlab’spinv 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 ofY0;(v)
the methodqrginv 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 dimensions1×1×1. The whole domain was discretized with 8 equal finite elements. Linear elasticity case was considered. The Young modulus E=1·106 and the Poisson ratio ν=0.3.


### 7.2.1. Test 4a

In this test, linear shape functions were used. The matrixK has dimensions 81×81. Results of Moore-Penrose pseudoinverse computations are shown in Table 4.Table 4
CPU time and relative errors, Test 4a.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0078 0.4808 ⋅ 1 0 - 03 0.2340 ⋅ 1 0 - 02 0.4283 ⋅ 1 0 - 02 0.4977⋅10-02 [19] 0.0280 0.7028⋅10-13 0.3091⋅10-09 0.1318⋅10-14 0.1547 ⋅ 1 0 - 14 [17] 0.0069 0.7028⋅10-13 0.1572⋅10-14 0.1809⋅10-11 0.5138 ⋅ 1 0 - 14


### 7.2.2. Test 4b

In the previous test, the SVD-based algorithm gave the worst results. It was due to the defaultTOL 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:Algorithm 1:Given K compute its pseudoinverse K† using pinv function.
SetTOLpinv=1, TOL=1⋅10-6 while  icond<4  do i c o n d = 0 T O L p i n v = 0.1 ⋅ T O L p i n v if  TOLpinv<1⋅10-20  then T O L p i n v = 1 ⋅ 1 0 - 1 and T O L = 10 ⋅ T O L end if K †=pinv(K,TOLpinv if  KK†K-KF/KF<TOL  then i c o n d = i c o n d + 1 end if if  K†KK†-K†F/K†F<TOL  then i c o n d = i c o n d + 1 end if if  KK†T-KK†F/KK†F<TOL  then i c o n d = i c o n d + 1 end if if  (K†K)T-K†KF/K†KF<TOL  then i c o n d = i c o n d + 1 end if end while return  K†Updated results for the SVD-based algorithm are shown in Table5.Table 5
CPU time and relative errors, Test 4b.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0064 0.7039 ⋅ 1 0 - 13 0.2704 ⋅ 1 0 - 14 0.6579 ⋅ 1 0 - 14 0.5456⋅10-14


### 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 matrixK 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.Table 6
CPU time and relative errors, Test 5 (Peano shape functions).
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 1.5543 0.7596 ⋅ 1 0 - 07 0.1946 ⋅ 1 0 - 10 0.6767 ⋅ 1 0 - 09 0.2657⋅10-10 [19] 175.44 0.0022⋅10+00 0.0018⋅10+00 0.0008⋅10+00 1.4142 ⋅ 1 0 + 00 [17] 1.4129 0.0001⋅10-09 0.0004⋅10-09 0.5837⋅10-09 0.1063 ⋅ 1 0 - 09Looking at Tables5 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 usingpinv and qrginv functions, respectively, K1† and K2†. Consequently, their pseudoinverses are K1 and K2. One may compute the following error indicators:(i)
In Test 4b:K-K1F/KF=1.0434·10-13 but K-K2F/KF=1, when K1-K2F/K1F=5.2436·10-12.(ii)
In Test 5:K-K1F/KF=1, K-K2F/KF=1 and K1-K2F/K1F=12.7059(!).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)(K†)†=K. In Test 4b, matrices K1 and K2 have almost the same norms, however, only the error indicator K-KcomputedF/KF computed for K1 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 (K-KcomputedF/KF) of order 1·10-12 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.Table 7
CPU time and relative errors, Test 5 (Legendre shape functions).
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 1.4112 0.0047 ⋅ 1 0 - 11 0.0344 ⋅ 1 0 - 11 0.1863 ⋅ 1 0 - 11 0.0671⋅10-11 [19] 101.35 0.0017⋅10+00 0.0018⋅10+00 0.0009⋅10+00 1.4142 ⋅ 1 0 + 00 [17] 1.8285 0.0005⋅10-10 0.0005⋅10-10 0.2127⋅10-10 0.0054 ⋅ 1 0 - 10Table 8
CPU time and relative errors, Test 5 (Lobatto shape functions).
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 1.4259 0.0357 ⋅ 1 0 - 12 0.1948 ⋅ 1 0 - 12 0.5772 ⋅ 1 0 - 12 0.2635⋅10-12 [19] 107.22 0.0013⋅10+00 0.0017⋅10+00 0.0009⋅10+00 1.4142 ⋅ 1 0 + 00 [17] 1.8454 0.0004⋅10-10 0.0005⋅10-10 0.1100⋅10-10 0.0020 ⋅ 1 0 - 10


## 7.2.1. Test 4a

In this test, linear shape functions were used. The matrixK has dimensions 81×81. Results of Moore-Penrose pseudoinverse computations are shown in Table 4.Table 4
CPU time and relative errors, Test 4a.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0078 0.4808 ⋅ 1 0 - 03 0.2340 ⋅ 1 0 - 02 0.4283 ⋅ 1 0 - 02 0.4977⋅10-02 [19] 0.0280 0.7028⋅10-13 0.3091⋅10-09 0.1318⋅10-14 0.1547 ⋅ 1 0 - 14 [17] 0.0069 0.7028⋅10-13 0.1572⋅10-14 0.1809⋅10-11 0.5138 ⋅ 1 0 - 14


## 7.2.2. Test 4b

In the previous test, the SVD-based algorithm gave the worst results. It was due to the defaultTOL 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:Algorithm 1:Given K compute its pseudoinverse K† using pinv function.
SetTOLpinv=1, TOL=1⋅10-6 while  icond<4  do i c o n d = 0 T O L p i n v = 0.1 ⋅ T O L p i n v if  TOLpinv<1⋅10-20  then T O L p i n v = 1 ⋅ 1 0 - 1 and T O L = 10 ⋅ T O L end if K †=pinv(K,TOLpinv if  KK†K-KF/KF<TOL  then i c o n d = i c o n d + 1 end if if  K†KK†-K†F/K†F<TOL  then i c o n d = i c o n d + 1 end if if  KK†T-KK†F/KK†F<TOL  then i c o n d = i c o n d + 1 end if if  (K†K)T-K†KF/K†KF<TOL  then i c o n d = i c o n d + 1 end if end while return  K†Updated results for the SVD-based algorithm are shown in Table5.Table 5
CPU time and relative errors, Test 4b.
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 0.0064 0.7039 ⋅ 1 0 - 13 0.2704 ⋅ 1 0 - 14 0.6579 ⋅ 1 0 - 14 0.5456⋅10-14


## 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 matrixK 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.Table 6
CPU time and relative errors, Test 5 (Peano shape functions).
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 1.5543 0.7596 ⋅ 1 0 - 07 0.1946 ⋅ 1 0 - 10 0.6767 ⋅ 1 0 - 09 0.2657⋅10-10 [19] 175.44 0.0022⋅10+00 0.0018⋅10+00 0.0008⋅10+00 1.4142 ⋅ 1 0 + 00 [17] 1.4129 0.0001⋅10-09 0.0004⋅10-09 0.5837⋅10-09 0.1063 ⋅ 1 0 - 09Looking at Tables5 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 usingpinv and qrginv functions, respectively, K1† and K2†. Consequently, their pseudoinverses are K1 and K2. One may compute the following error indicators:(i)
In Test 4b:K-K1F/KF=1.0434·10-13 but K-K2F/KF=1, when K1-K2F/K1F=5.2436·10-12.(ii)
In Test 5:K-K1F/KF=1, K-K2F/KF=1 and K1-K2F/K1F=12.7059(!).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)(K†)†=K. In Test 4b, matrices K1 and K2 have almost the same norms, however, only the error indicator K-KcomputedF/KF computed for K1 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 (K-KcomputedF/KF) of order 1·10-12 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.Table 7
CPU time and relative errors, Test 5 (Legendre shape functions).
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 1.4112 0.0047 ⋅ 1 0 - 11 0.0344 ⋅ 1 0 - 11 0.1863 ⋅ 1 0 - 11 0.0671⋅10-11 [19] 101.35 0.0017⋅10+00 0.0018⋅10+00 0.0009⋅10+00 1.4142 ⋅ 1 0 + 00 [17] 1.8285 0.0005⋅10-10 0.0005⋅10-10 0.2127⋅10-10 0.0054 ⋅ 1 0 - 10Table 8
CPU time and relative errors, Test 5 (Lobatto shape functions).
Method Time [s] KK†K-KFKF K † K K † - K † F K † F KK†T-KK†FKK†F K † K T - K † K F K † K F SVD 1.4259 0.0357 ⋅ 1 0 - 12 0.1948 ⋅ 1 0 - 12 0.5772 ⋅ 1 0 - 12 0.2635⋅10-12 [19] 107.22 0.0013⋅10+00 0.0017⋅10+00 0.0009⋅10+00 1.4142 ⋅ 1 0 + 00 [17] 1.8454 0.0004⋅10-10 0.0005⋅10-10 0.1100⋅10-10 0.0020 ⋅ 1 0 - 10


## 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(K†)†=K 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.



---

*Source: 5060397-2019-02-24.xml*