Generalized Inverses Estimations by Means of Iterative Methods with Memory

A secant-type method is designed for approximating the inverse and some generalized inverses of a complex matrix A. For a nonsingular matrix, the proposed method gives us an approximation of the inverse and, when the matrix is singular, an approximation of the Moore–Penrose inverse and Drazin inverse are obtained. The convergence and the order of convergence is presented in each case. Some numerical tests allowed us to confirm the theoretical results and to compare the performance of our method with other known ones. With these results, the iterative methods with memory appear for the first time for estimating the solution of a nonlinear matrix equations.


Introduction
Recently, many iterative methods without memory have been published for approximating the inverse or some generalized inverse of a complex matrix A of arbitrary order (see, for example, [1][2][3][4][5][6] and the references therein). This topic has a significant role to play in many areas in applied sciences and engineering, such as multivariate analysis, image and signal processing, approximation theory, cryptography, etc. (see [7]).
The discretization process of boundary problems or partial differential equations by means of divided difference technique or finite elements yields to an important number of linear systems being solved. This statement is applicable both in equations with integer derivatives and in the case of fractional derivatives (see, for example, [8,9]). In these linear problems, usually the matrix of coefficients is too big or ill-conditioned to be solved analytically. Thus, iterative methods can play a key role.
The main purpose of this manuscript is to design a secant-type iterative scheme with memory, free for inverse operators and efficient under the point of view of CPU-time, for estimating the inverse of a non-singular complex matrix. We also argue the generalization of the proposed scheme for approximating the Drazin inverse of singular square matrices and the Moore-Penrose inverse of complex rectangular matrices. As far as we know, this is the first time that this kind of methods with memory is applied to estimate generalized inverses. This might be the first step to develop higher-order methods with memory in the future. This kind of schemes has proven to be very stable for scalar equations; we expect a similar performance in the case of matrix equations.
Let us consider a non-singular complex matrix A of size n × n. The extension of the iterative methods for the real equation g(x) = ax − 1 = 0 to obtain the inverse of A, that is the zero of the matrix function G(X) = X −1 − A, gives us the so-called Schulz-type schemes.
The most known of these schemes to estimate A −1 is the Newton-Schulz method [10], whose iterative expression is X k+1 = X k (2I − AX k ), k = 0, 1, . . . , (1) where I denotes the identity matrix of order n. Schulz [11] demonstrated that the eigenvalues of matrix I − AX 0 must be less than 1 to assure the convergence of the scheme in Equation (1). Taking into account that the residuals Newton-Schulz method has quadratic convergence. In general, it is known that this scheme converges to A −1 with X 0 = αA * or X 0 = αA, where 0 < α < 2/ρ(A * A), ρ(·) denotes the spectral radius, and A * is the conjugate transpose of A. Such schemes are also used for sensitivity analysis when accurate approximate inverses are needed for both square and rectangular matrices. On the other hand, for a nonsingular matrix A ∈ C n×n , Li et al. [12] suggested the scheme with X 0 = αA * . They proved the convergence of m-order of {X k } to the inverse of A. This result was extended by Chen et al. [13] for computing the Moore-Penrose inverse. Other iterative schemes without memory have been designed for approximating the inverse or some generalized inverses.
In this paper, we construct an iterative method with memory (that is, k + 1 iterate is obtained not only from the iterate k but also from other previous iterates) for computing the inverse of a nonsingular matrix. In the iterative expression of the designed method, inverse operators do not appear. We prove the order of convergence of the proposed scheme and we extend it for approximating the Moore-Penrose inverse of rectangular matrices and the Drazin inverse of singular square matrices.
For analyzing the order of convergence of an iterative method with memory, we use the concept of R-order introduced in [14] by Ortega and Rheinboldt and the following result.
Let us consider an iterative method with memory (IM) that generates a sequence {X k } of estimations to the solution ξ, and let us also assume that this sequence converges to ξ. If there exists a nonzero constant η and nonnegative numbers t i , 0 ≤ i ≤ m such that the inequality holds, where e k is the error of iterate X k , then the R-order of convergence of (IM) satisfies where s * is the unique positive root of the polynomial The proof of this result can be found in [14]. From here, the work is organized as follows. In the next section, we describe how a secant-type method, free of inverse operators, is constructed for estimating the inverse of a nonsingular complex matrix, proving its order of convergence. In Sections 3 and 4, we study the generalization of the proposed methods for computing the Moore-Penrose inverse of a rectangular complex matrix and the Drazin inverse of a singular square matrix. Section 5 is devoted to the numerical test for analyzing the performance of the proposed schemes and to confirm the theoretical results. With a section of conclusions, the paper is finished.

A Secant-Type Method for Matrix Inversion
Let us recall that, for an scalar nonlinear equation g(x) = 0, the secant method is an iterative scheme with memory such that , k ≥ 0, given x 0 and x −1 as initial approximations. For a nonlinear matrix equation G(X) = 0, where G : C n×n → C n×n , the secant method can be described as where X 0 and X −1 are initial estimations and being A k a suitable linear operator satisfying Thus, it is necessary to solve, at each iteration, the linear system A k+1 S k = Y k . It is proven in [15] that, with this formulation, secant method converges to the solution of G(X) = 0.
Let us consider an n × n nonsingular complex matrix A. We want to construct iterative schemes for computing the inverse A −1 of A, that is, iterative methods for solving the matrix equation The secant method was adapted by Monsalve et al. [15] to estimate the solution of Equation (5), that is the inverse of A, when the matrix is diagonalizable. The secant method applied to G(X) = X −1 − A (see [15]) gives us: Now, we extend the result presented in [15] to any nonsingular matrix, not necessarily diagonalizable. If A is a nonsingular complex matrix of size n × n, then there exist unitary matrices U and V, of size n × n, such that being Then, from Equation (6), Several algebraic manipulations allow us to assure that If we choose initial estimations, X −1 and X 0 , such that D −1 = V * X −1 U and D 0 = V * X 0 U are diagonal matrices, then all matrices D k are diagonal and therefore D i D j = D j D i , for all i, j. Thus, from Equation (8), we assure and, from this expression, we propose the secant-type method: being X 0 and X −1 initial approximations given. The analysis of the convergence of the iterative method with memory in Equation (9) is presented in the following result. Theorem 1. Let A ∈ C nxn be a nonsingular matrix, with singular value decomposition U * AV = Σ. Let X 0 and X −1 be such that V * X −1 U and V * X 0 U are diagonal matrices. Then, sequence {X k }, obtained by Equation (9), converges to A −1 with super-linear convergence.
Proof. Let us consider U and V unitary matrices such that the singular values decomposition in Equation (7) By subtracting 1 σ j from both sides of Equation (10) and denoting e j From Equation (11), we conclude that, for each value of j from 1 to n, d j k+1 in Equation (10) Thus, Therefore, which allows us to affirm that {X k } converges to A −1 .
On the other hand, Highan in [10] introduced the following definition for the stability of the iterative process Z k+1 = H(Z k ), with a fixed point Z * . If we assume that H is Frechét differentiable in Z * , the iteration is stable in a neighborhood of Z * if the Frechét derivative H (Z * ) has bounded powers, that is, there exists a positive constant C such that Therefore, the following result can be stated for the secant method. (9) for the estimation of inverse matrix is a stable iterative scheme.

Theorem 2. The secant method in Equation
Proof. The proof is made demonstrating that H (Z * ) is an idempotent matrix.
The secant-type method described as a fixed point scheme, can be written as It is easy to deduce that Thus, H (Z * ) is an idempotent matrix and the iteration is stable.

A Secant-Type Method for Approximating the Moore-Penrose Inverse
Now, we would like to extend the proposed iterative scheme for computing the Moore-Penrose inverse [7] of a m × n complex matrix A, denoted by A † . It is the unique n × m matrix X satisfying the equations AXA = A, XAX = X, (AX) * = AX, (XA) * = XA.
The convergence of the method in Equation (9) for Moore-Penrose inverse is established in the following result.

Theorem 3.
Let A ∈ C m×n be a matrix with rank(A) = r, with singular value decomposition Let X −1 and X 0 be initial estimations such that being Σ −1 and Σ 0 diagonal matrices of size r × r. Then, sequence {X k }, obtained by Equation (9), converges to A † with super-linear order of convergence.
Proof. Given the singular value decomposition of A, for any fixed arbitrary value of k, we define matrix D k as being Σ k ∈ C r×r . Thus, by using the iterative expression in Equation (9), we obtain Therefore, as Σ −1 and Σ 0 are diagonal matrices, so are all matrices Σ k , and the expression represents r scalar uncoupled iterations converging to 1 with M k = max 1≤i≤r {c 2 k }, being c i k > 0 such that sequence {c i k } tends to zero for k tending to infinity. With an analogous argument as in Theorem 1, which allows us to affirm that {X k } converges to A † , with the desired order of convergence.

A Secant-Type Method for Approximating the Drazin Inverse
Drazin, in 1958 (see [10]), proposed a different kind of generalized inverse, in which some conditions of the Moore-Penrose inverse and the index of the matrix appeared. The importance of this inverse has motivated many researchers to propose algorithms for its calculation.
It is known (see [10]) that the smallest nonnegative integer l, such that rank(A l+1 ) = rank(A l ) is called the index of A and it is denoted by ind(A). If A is a complex matrix of size n × n, the Drazin inverse of A, denoted A D , is the unique matrix X satisfying where l is the index of A.
If ind(A) = 1, then X is called the g-inverse or group inverse of A, and, if ind(A) = 0, then A is nonsingular and A D = A −1 . Let us observe that the idempotent matrix AA D is the projector on R(A l ) along N (A l ), where R(A l ) and N (A l ) denote the range and null spaces of A l , respectively.
In [16], the following result is presented, which is used in the proof of the main result. Li and Wei [1] proved that the Newton-Schulz method in Equation (1) can be used for approximating the Drazin inverse, using as initial estimation X 0 = αA l , where parameter α is chosen so that condition I − AX 0 < 1 is satisfied. One way for selecting the initial matrix used by different authors is where tr(·) is the trace of a square matrix. Another fruitful initial matrix is Using two initial matrices of these form, αA l , with α a constant, we want to prove that the sequence obtained by the secant-type method in Equation (9) converges to the Drazin inverse A D . In this case, we use a different type of demonstration than those used in the previous cases.

Theorem 4.
Let A ∈ C n×n be a square nonsingular matrix. We choose as initial estimations X 0 = α 0 A l 0 and X 1 = α 1 A l 1 , with l 0 , l 1 ≥ ind(A). Then, sequence {X k } k≥0 generated by Equation (9) satisfies the following error equation Thus, {X k } k≥0 converges to A D with order of convergence 1.618, that is, with super-linear convergence.
Proof. Let us define E k = I − AX k , k = 0, 1, . . .. Then, Therefore, E k+1 ≤ E k−1 E k . In addition, it is easy to prove that, if we choose X 0 and X 1 such that E 0 < 1 and E 1 < 1, then E k < 1, ∀k ∈ N. Now, we denote e k = A D − X k the error of iterate k. From the selection of X 0 and X 1 and by applying Proposition 1, we establish From this identity, there exists k 0 ∈ N such that Thus, { e k } k≥0 tends to zero and therefore {X k } k≥0 tends to A D .
On the other hand, Now, we analyze Ae k+1 .
In the last equality, we use that In addition, (I − AA D )Ae k = 0 and Ae k−1 (I − AA D ) = 0. Therefore, Finally, by applying the theorem of convergence for iterative methods with memory, as mentioned in the Introduction, we assure that the order of convergence of secant-type method is the unique positive root of λ 2 − λ − 1 = 0, that is λ = 1.618.

Numerical Experiments
In this section, we check the behavior for the calculation of the inverse, Moore-Penrose inverse and Drazin inverse, of different test matrices A, using the secant method, which we compared with the Newton-Schulz scheme in Equation (1). Numerical computations were carried out in Matlab R2018b (MathWorks, Natick, USA) with a processor Intel(R) Xeon(R) CPU E5-2420 v2 at 2.20 GHz. As stopping criterion, we used X k+1 − X k 2 < 10 −6 or F(X k+1 ) 2 < 10 −6 .
To numerically check the theoretical results, Jay [17] introduced the order of approximate computational convergence (COC), defined as In a similar way, the authors presented in [18] another numerical approximation of the theoretical order, denoted by ACOC, and defined as We use indistinctly any of these computational order estimates, to show the accuracy of these approximations on the proposed method. In the case of vector COC (or ACOC) is not stable, we write "-" in the corresponding table.
Example 1. In this example, matrix A is a n × n random matrix with n = 10, 100, 200, 300, 400, 500. The initial estimation used for the Newton-Schulz scheme is X 0 = A T / A 2 and for the secant method In Table 1, we show the results obtained by Newton-Schulz and secant-type method for the different random matrices, the number of iterations, the residuals, and the value of COC. The results are in concordance with the order of convergence of each scheme. All obtained random matrices are nonsingular and both methods give us an approximation of the inverse of A. Newton method needs lower number of iterations than Secant scheme, as was expected, being the first one quadratic and the latter one super-linear. Example 2. In this example, matrix A is a m × n random matrix for different values of m and n. The initial matrices are calculated in the same way as in the previous example.
In Table 2, we show the results obtained by Newton-Schulz and secant-type method for the different random matrices, the number of iterations, the residuals, and the value of ACOC. The results are in concordance with the order of convergence of each scheme, despite being non-square matrices. Both methods give us an approximation of the Moore-Penrose inverse of A. Example 3. In this example, we want to analyze the performance of the secant method for computing the Drazin inverse of the following matrix A of size 6 × 6 with ind(A) = 2.
Here, its Drazin inverse is expressed by On the other hand, secant method is used with , obtaining: Example 4. This is another example for computing the Drazin inverse of the following matrix B of size 12 × 12 (see [1]) with ind(B) = 3.

Now, its Drazin inverse is expressed by
By using the initial matrix X 0 = 0.5 tr(A 5 ) and the same stoping criterion as in the previous examples, Newton-Schulz method gives us the following information: • ACOC = 2.0031; • iter = 14; and • Exact error B D − X 14 = 1.8354 × 10 −9 .
On the other hand, secant method is used with X −1 = 1 tr(A 5 ) and X 0 = 0.5 tr(A 5 ) , obtaining: Again, the numerical tests confirm the theoretical results.
Example 5. Finally, in this example, we test Newton-Schlutz and secant methods on several known square matrices of size n × n, constructed by using different Matlab functions. Specifically, the test matrices are: (a) A = gallery( ris , n). Hankel matrix of size n × n. (b) A = gallery( grcar , n). Toeplitz matrix of size n × n.
(c) A = gallery( lehmer , n). Symmetric and positive definite matrix of size n × n, a i,j = i/j, ∀i, j.
(d) A = gallery( leslie , n). Leslie matrix of size n × n. This type of matrices appears in problems of population models. (e) A = gallery( invo , n). Matrix ill-conditioned of size n × n, such that A 2 = I. By using the stopping criterion X k+1 − X k 2 < 10 −10 or F(X k+1 ) 2 < 10 −10 and the initial matrix X −1 = A T A 2 and X 0 = 0.5 A T A 2 , we obtain the numerical results that appear in Table 3. In this cases, as in the previous ones, the proposed method shows good performance in terms of stability, precision, and number of iterations needed. We must take into account that both schemes have different orders of convergence, which is displayed in Table 3. Table 3. Results for approximating the inverse of classical square matrices (Example 5)

Conclusions
An iterative method with memory for approximating the inverse of nonsingular square complex matrices, the Moore-Penrose inverse of rectangular complex matrices, and the Drazin inverse of square singular matrices is presented. As far as we know, it is the first time that a scheme with memory is employed to approximate the solution of nonlinear matrix equations. The proposed scheme is free of inverse operators and its iterative expression is simple; therefore, it is computationally efficient. From particular initial approximations, the convergence is guaranteed for all matrices, without conditions. Numerical tests allowed us to analyze the performance of the proposed scheme and confirm the theoretical results.