Improving Newton–Schulz Method for Approximating Matrix Generalized Inverse by Using Schemes with Memory

: Some iterative schemes with memory were designed for approximating the inverse of a nonsingular square complex matrix and the Moore–Penrose inverse of a singular square matrix or an arbitrary m × n complex matrix. A Kurchatov-type scheme and Steffensen’s method with memory were developed for estimating these types of inverses, improving, in the second case, the order of convergence of the Newton–Schulz scheme. The convergence and its order were studied in the four cases, and their stability was checked as discrete dynamical systems. With large matrices, some numerical examples are presented to conﬁrm the theoretical results and to compare the results obtained with the proposed methods with those provided by other known ones.


Introduction
Calculating the inverse of a matrix, especially of a large size, is a very difficult task with a high computational cost.The alternative is the use of iterative algorithms to estimate it.In a vast range of fields, such as image and signal processing [1][2][3][4], encryption [5,6], control system analysis [7,8], etc., it is necessary to calculate the inverse or different generalized inverses in order to solve the problems posed.
In recent years, many iterative schemes of different orders of convergence have been designed to estimate the inverse of a complex matrix A or some generalized inverse (Moore-Penrose inverse, Drazin inverse, etc.).In 2013, Weiguo et al. in [9] constructed a sequence of third-order iterations converging to the Moore-Penrose matrix.In the same year as well, Toutounian and Soleymani in [10] presented a high-order method for approximating inverses and pseudo-inverses of complex matrices based on Homeier's scheme, with a derivative-free composition.More recently, Stanimirovic et al. in [11] designed efficient transformations of the hyperpower iterative method for computing the generalized inverses T,S with the aim of minimizing the number of required matrices by matrix products per cycle.In 2020, Kaur et al. established in [12] new formulations of the fifth-order hyperpower method to compute the weighted Moore-Penrose inverse, improving the efficiency indices.Such approximations were found to be robust and effective when implemented as a preconditioned matrix to solve the linear systems.All these schemes were designed with the starting point of iterative procedures without memory, that is where each new iteration is calculated using only the information provided by the previous iterate.
Iterative procedures that use more than one previous iterate to calculate the next one are called methods with memory.In the context of matrix inverse approximations, the secant scheme was proposed by the authors in [13].For a nonsingular matrix, the secant method gives an estimation of the inverse and, when the matrix is singular, an approximation of the pseudo-inverse and Drazin inverse.Furthermore, the superlinear convergence was demonstrated in all cases.
In this manuscript, we focused on constructing several iterative methods with memory, free of inverse operators and with different orders of convergence, for finding the inverse of a nonsingular complex matrix.We also analyzed the proposed schemes for computing the Moore-Penrose inverse of complex rectangular matrices.On the other hand, these procedures allow approximating other generalized inverses such as the Drazin inverse, the group inverse, etc., although this is beyond the scope of this work.
Let A be a complex n × n nonsingular matrix.The design of iterative algorithms without memory for estimating the inverse matrix, which we call Newton-Schulz-type methods, is mostly based on iterative solvers for the scalar equation f (x) = 0 applied to the nonlinear matrix equation: where F : C n×n → C n×n is a nonlinear matrix function.
The most-used iterative scheme to approximate A −1 is the Newton-Schulz scheme [14]: where I is an identity matrix of size n × n.
On the other hand, in the context of iterative procedures with memory, the authors presented in [13] a secant-type method (SM), whose iterative expression is: where X 0 , X −1 are the starting guesses.With a particular choice of the initial approximations, the authors proved that sequence {X k }, obtained by (3), converges to A −1 with order of convergence 1.618.
To analyze the convergence order of the iterative methods with memory for solving nonlinear equations f (x) = 0, the R-order is used (see [15]), which we summarize below.
Theorem 1.Let ψ be an iterative method with memory that generates a sequence {x k } of approximations to the root α of f (x) = 0, and let this sequence converge to α.If there exists a nonzero constant η and nonnegative numbers t i , i = 0, 1, . . ., m, such that the inequality: holds, then the R-order of convergence of the iterative method ψ satisfies the inequality: where s * is the unique positive root of the equation: Here, e k = x k − α denotes the error of the approximation in the kth iterative step.Continuing with the idea of designing iterative procedures with memory, in this work, we used the scalar iterative methods of Kurchatov and Steffensen with memory, which we adapted to the context of matrix equations.In the scalar case, this kind of procedure can reach a higher order of convergence than schemes without memory with the same number of functional evaluations and usually has better stability properties.However, as far as we know, only the secant scheme has been extended to matrix equations with good results, in [13].
Kurchatov's scheme is an iterative algorithm with memory with quadratic convergence to solve scalar equations f (x) = 0 (see [16]), which is deduced from Newton's method by substituting f (x) by Kurchatov's divided difference, that is: Regarding Steffensen's scheme, it was initially developed in [17], where the derivative f (x) was replaced by the divided difference f [x k , x k + f (x k )].It still holds the original quadratic convergence, but in practice, the set of converging initial estimations is quite smaller than that of Newton's scheme.It was Traub who, in [18], introduced an accelerating parameter γ in the divided difference such that f ( ] derived an iterative method whose error equation was: With the aim of increasing the order of convergence of the scheme, Traub defined an approximation of , obtaining the procedure: which, given initial values x 0 and γ 0 , has order of convergence 1 + √ 2. The stability of these schemes with memory for scalar problems was analyzed firstly in [19,20], showing very stable performance in both cases.In what follows, some definitions and properties about vectorial discrete dynamical systems are introduced.These tools are shown to be useful in the following sections.

Basic Concepts of Qualitative Studies of Schemes with Memory
Let us start with an iterative scheme with memory, which is used to solve the scalar problem f (x) = 0, using two previous iterates to calculate the next one: with x 0 and x 1 being its seeds.Therefore, a solution α is estimated whether x k+1 = x k or, equivalently, x k = ψ(x k−1 , x k ).This estimation can be obtained as a fixed point of Ψ : R 2 −→ R 2 by means of: with again x 0 and x 1 being the seeds.
To study the qualitative behavior of a fixed point iterative scheme with memory, it is applied on a low-degree polynomial p(x), generating a vectorial rational operator, G(x).Then, the orbit of (z, x) ∈ R 2 is the set: {(z, x), G((z, x)), . . . ,G m ((z, x)), . ..}.
Therefore, a point (z, x) is fixed if it is the only point belonging to its orbit, and it is a T-periodic point if is orbit is only defined by (z, x), G((z, x)), . . ., G T−1 ((z, x)) and G T ((z, x)) = (z, x).
On the other hand, the orbit of a point can be classified depending on its asymptotic behavior, according to the next result.
Moreover, we define a point x c ∈ R 2 as a critical point of G if det(G (x c )) = 0. Indeed, if x c is not a zero of p(x), it is named a free critical point.In a similar way, if a fixed point is not a zero of p(x), it is named a strange fixed point.
Denoting by x F ∈ R 2 an attracting fixed point of: G, its basin of attraction B(x F ) is defined as: The union of all the basins of attraction of G is defined as the Fatou set, F (G), and its complementary set in R 2 is the Julia set J (G).It holds all the repelling fixed points and sets the boundary among the basins of attraction.
In Sections 2 and 3, we build a Kurchatov-type and a Steffensen-type iterative method with memory, respectively, to estimate the inverse of a nonsingular complex matrix.Both schemes are free of inverse operators, and we prove their order of convergence and stability.In Section 4, we extend these schemes to calculate the Moore-Penrose inverse of a complex rectangular matrix.Section 5 is devoted to the numerical tests to analyze their performance and confirm the theoretical results.We end the work with some conclusions.

Kurchatov-Type Method
The Kurchatov divided difference was defined by [16] as: .
For an equation f (x) = 0, Kurchatov's method has order or convergence two and iterative expression: .
If we apply this method to F(X) = X −1 − A, where F : C n×n → C n×n and A ∈ C n×n is an n × n nonsingular matrix, not necessarily diagonalizable, we obtain: Next, we show a technical lemma whose result we use later.
Lemma 1.Let A, B ∈ C n×n be invertible and commutative, AB = BA, then: Proof.Using algebraic manipulations, and the lemma is proven.
It is known that, for any nonsingular n × n complex matrix A, there exist n × n unitary matrices U and V, such that U * AV = Σ, Σ being the diagonal matrix of the singular values of , and we apply them in Equation ( 5). where then all the matrices D k are diagonal and D i D j = D j D i , for all i, j.Applying several algebraic manipulations on the last equation, we can ensure: Therefore, Using Lemma 1, we obtain: Then, Finally, In this iterative expression, still, an inverse matrix appears.Therefore, Kurchatov's scheme is not directly transferable for the calculation of the matrix inverses.
Therefore, we propose a slight change in the Kurchatov divided difference scheme: obtaining a new iterative scheme, . Now, we apply this iterative procedure to the matrix equation F Using the transformations defined above and Lemma 1, we obtain: and then, From this expression, by using again the transformations defined above, we have where three matrix products appear, but there are no inverse operators.Equation ( 8) corresponds to the iterative scheme of the modified Kurchatov-type method for matrix inversion, which we call MKTM.

Convergence Analysis
Let us consider now two n × n unitary matrices U and V, satisfying By using Equation ( 7), we obtain: Next, we demonstrate the order of convergence of the iterative procedure MKTM.Theorem 3. Let A ∈ C n×n be a nonsingular matrix and X −1 and X 0 be the initial approximation such that matrices V * X −1 U and V * X 0 U are diagonal.Then, sequence {X k } obtained by means of the iterative expression (8), where Σ = U * AV is the singular-value decomposition of A, converges to the inverse A −1 , with convergence order 1.6180.

Proof. Let us denote
. By means of component-by-component calculations, we obtain: Let us also denote e j k = d j k − 1/σ j ; by subtracting 1 σ j from the both sides of the last equation, we obtain: This result shows that, for each j = 1, 2, . . ., n, d k+1 in Equation ( 9) converges to 1 σ j with order of convergence (1 + √ 5)/2 ≈ 1.6180, which is the only positive root of λ 2 − λ − 1 = 0 (see Theorem 1).In this way, for each j = 1, 2, . . ., n, there exists {C Using this result, we obtain: Therefore, we can affirm that {X k } converges to A −1 .
To demonstrate the stability of the modified Kurchatov-type iterative scheme, we used the definition introduced by Higham in [14] for the stability of the iterative scheme Z k+1 = H(Z k ), with a fixed point Z * .Assuming that H is Frechét differentiable in fixed point Z * , then the process is stable in a neighborhood of Z * if there exists a positive constant Theorem 4. The modified Kurchatov-type iterative scheme to estimate the inverse of a nonsingular complex matrix defined by expression (8) is stable.

Qualitative Analysis of Kurchatov-Type Scheme
Now, we considered the Kurchatov-type scheme as a scalar iterative procedure for solving nonlinear equations, analyzed its convergence, and studied its stability.In this analysis, the stability is understood as the dependence on the initial estimations.It is made by applying the vectorial real dynamical concepts to the Kurchatov-type scheme, whose expression is: In the following result, we show that this scheme for solving nonlinear scalar equations f (x) = 0 has superlinear convergence.Theorem 5. Let us consider a sufficiently differentiable function f : I ⊆ R → R in an open neighborhood I of the simple root α of nonlinear equation f (x) = 0. Furthermore, let us assume that f (x) is continuous at α and that x 0 and α are near enough between them.Then, sequence x k , k = 0, 1, 2, . .., generated by the Kurchatov-type scheme, converges to α with order of convergence 2 , its error equation being: where O 2 (e k−1 , e k ) means that the following terms in the error equation depend on the powers of errors e k−1 and e k such that the sum of the exponents is at least two; moreover, e k = x k − α and Proof.It is known that the Taylor expansion of f (x k ) around α is: By using the Genocchi-Hermite formula (see [15]) and the expansion of f (x + th) in the Taylor series around x, with , the expansion of the divided difference is: Then, the expression of the difference divided as a function of the errors e k−1 and e k is: Therefore, By applying Theorem 1, the only positive real root of λ 2 − λ − 1 = 0 (where the coefficients correspond to the exponents of e k−1 and e k in the error equation) is the convergence order of the method, that is 1+ From now on, we denote by KT the fixed point operator associated with the Kurchatovtype method applied on p(z) = z 2 − 1.As it does not use derivatives, it is not possible to establish a scaling theorem.Therefore, we cannot make the analysis on generic seconddegree polynomials.
The fixed point operator depends on two variables: x k (denoted by x) and x k−1 (denoted by z).Therefore, Let us analyze now the qualitative performance of the rational operator by means of the asymptotic behavior of the fixed points and the existence of free critical points.
Theorem 6.The only fixed points of rational operator KT are the roots x = z = ±1, both being superattracting.Moreover, KT has no free critical points.
From this results, we concluded that no other performance than convergence to the roots is possible, when the Kurchatov-type is applied on p(x), as any other basin of attraction would need a free critical point inside.
In Figure 1, we show the dynamical plane of the KT operator, x and z being real, corresponding to the abscissa and ordinate axis, respectively (see [22] for the routines).The mesh used has 1000 × 1000 points, and the maximum amount of iterations is 40.The distance to the root lower than 10 −3 is the stopping criterion.
Then, if each point of the mesh is considered as a seed of the method, when it converges to one of the roots of p(z) (located in x = z), it is represented in orange or green color, depending on the root it has converged to.The color is brighter the lower the amount of iterations is.When the initial estimation, after the iterative process, reaches 40 iterations with no convergence, it is colored in black.
In Figure 1, we observe that the basins of attraction of both roots of p(z) have symmetric shapes.In spite of the only basins of attraction being those of the roots (Theorem 6), there are black areas in the dynamical planes.They correspond to areas of slow convergence, whose initial estimations converge with a higher number of iterations.In [19,20], the stability analysis of the secant and Steffensen with memory schemes was performed, among other iterative schemes with memory.In both cases, the dynamical planes plotted showed full convergence to the roots, without black areas of slow convergence or divergence.Later on, we checked how the the differences among the Kurchatov-type and secant or Steffensen with memory stability affected their numerical performance.

Steffensen with Memory
Now, we considered the Traub-Steffensen family [18], whose iterative expression is: where γ is an arbitrary constant different from zero.Furthermore, if γ = 1 in this iterative scheme, we obtain the well-known Steffensen method [17].It is easy to show that this iterative procedure has order of convergence two for any value of γ, its error equation being: where e k = x k − α represents the error in each iteration and c j = 1 j! f (j) (α)/ f (α), j ≥ 2. Let us also remark that this is a derivative-free scheme without memory, because it uses only the information of the current iteration.Studying the equation of the error (12) of family (11), we notice that, if γ = −1/ f (α), the iterative scheme (11) achieves the third order of convergence.This value of the parameter cannot be used, as α is unknown.Therefore, the approximation of f (α) as the first-order divided differences , which is the secant scheme, and ( 11) is an iterative procedure with memory.

Iterative Scheme Design for the Matrix Inverse
Our aim now was to adapt Steffensen's scheme with memory to solve the matrix Equation (1).In [23] Monsalve et al. adapted the secant method to estimate the solution of Equation (1) for A diagonalizable, and in [13], the authors extended this result to any nonsingular matrix by means of the expression: where W 0 and W −1 are the seeds.
In a similar way, Expression (11) for matrix equations takes the form: which includes the inverse calculations for the estimation of A −1 .This circumstance must be avoided.
As A is a nonsingular n × n complex matrix, there exist n × n unitary matrices U and Then, Equation ( 14) takes the form: If the initial approximations X −1 and X 0 satisfy D −1 = V * X −1 U and D 0 = V * X 0 U, then D k are diagonal matrices and D i D j = D j D i for all i, j.By applying several algebraic manipulations on Equation ( 15), we can ensure: and therefore, where I denotes the identity matrix of size n.Now, taking k V and after more algebraic manipulations, we obtain: given the initial approximations X 0 and X −1 .We denote this iterative scheme with memory as SMM.Let us remark that Expression (18) does not include inverse operators, as intended.

Convergence Analysis and Stability
We present in the next result the convergence analysis of the iterative method given by Expression (18).Theorem 7. Let U * AV = Σ be the singular-value decomposition of the a nonsingular matrix A ∈ C n×n .Furthermore, let X 0 and X −1 be the initial approximations to A −1 , such that V * X −1 U and V * X 0 U are diagonal matrices.Then, sequence {X k }, generated by (18), converges to A −1 with order of convergence 1 + √ 2 ≈ 2.4142.

Proof.
Let U and V be unitary matrices such that U * AV = Σ, Σ being the diagonal matrix of the singular values . ., from Equation ( 18), we obtain: Let us remark that Therefore, in Equation ( 19), using component-by-component calculations, we obtain: We denote e j k = d j k − 1/σ j .We subtract 1 σ j from both sides of Equation ( 20), and we obtain: Theorem 1 allows us to conclude that, for each j from 1 to n, d k+1 converges to 1 σ j with the convergence order the only positive root of [15]).In this way, for each j = 1, 2, . . ., n, there exist a sequence {C On the other hand, where Using this result, we obtain: Then, {X k } converges to A −1 , and the proof is finished.
In the following result, we prove the stability of this scheme.
Theorem 8. Steffensen's method with memory to estimate the inverse of a given matrix with the expression (18) is a stable iterative scheme.
Proof.Steffensen's scheme with memory can be written as a fixed point scheme as follows: Then, given P = (P 1 , P 2 ) T , we deduce that: Next, for Z = Z * = (A −1 , A −1 ) T , we have: S (Z * )P = 0 0 I 0 therefore, S (Z * ) is an idempotent matrix, and the iterative process ( 18) is stable.

Approximating the Moore-Penrose Inverse
For a complex matrix that is not necessarily square, the pseudo-inverse can be defined as an object similar to an inverse matrix.Given a complex matrix, it is possible to define many possible pseudo-inverses, but the most-frequent one is the Moore-Penrose inverse.
In cases when a system of equations has no solution, there is no inverse of the matrix of the coefficients that defines the system.In these cases, it may be useful to find a value that is an approximate solution in terms of error minimization, such as being able to find the best fit of a dataset using the pseudo-inverse (Moore-Penrose inverse) of the coefficient matrix.See, for example, the classical texts of Ben-Israel [24] and Higham [14].
We now extend the proposed iterative methods of this manuscript for calculating the Moore-Penrose inverse A † of an m × n complex matrix A. This is the unique n × m matrix X that satisfies the following conditions: A computationally simple and accurate way to compute the Moore-Penrose inverse is through the use of the singular-value decomposition of A. Then, if the rank of A is r ≤ min(m, n), we obtain: where Σ = diag(σ 1 , σ 2 , . . ., σ r ), σ 1 ≥ σ 2 ≥ . . .≥ σ r > 0, with U ∈ C m×m and V ∈ C n×n a unitary matrix.Furthermore, it is known that: where Theorem 9. Let A ∈ C m×n be a matrix with singular-value decomposition satisfying (22) whose rank(A) = r.If X −1 and X 0 are the initial estimations such that: given that Σ −1 and Σ 0 are diagonal matrices of size r × r, then the sequences {X k }, generated by Equation (18) using SSM and (8) using MKTM, converge to A † with order of convergence 1 + √ 2 and (1 + √ 5)/2, respectively.
Proof.Given the singular-value decomposition of A (see ( 22)), we define the matrix D k : for any fixed arbitrary value of k, with Σ k ∈ C r×r .Now, by using the iterative expression in Equation ( 18) and ( 8), we obtain: respectively.Because Σ −1 and Σ 0 are diagonal matrices, then all matrices Σ k are diagonal as well.Then, the expression with Just as was argued in (8), we have that: Therefore, we can say that {X k } converges to A † , with the desired order of convergence.

Numerical Experiments
In this section, we present the numerical tests of the behavior of the Steffensen method with memory (SMM) and Kurchatov type's (MKTM) designed to calculate the inverse and the Moore-Penrose inverse applied to different matrices.For comparison, we used the Newton-Schulz method (NS) [14], and the secant method (SM) [13].The numerical calculations were made with MATLAB 2022a (MathWorks, USA) using a 3 GHz 10-Core Intel Xeon W processor, 64 GB 2666 MHz DDR4 (iMac Pro).As stopping criteria for all numerical tests, we used ||X k+1 − X k || 2 < 10 −6 or ||F(X k+1 )|| 2 < 10 −6 , F(X) = 0 being the nonlinear matrix equation to be solved for estimating the inverse or pseudo-inverse of a complex matrix A.
In addition, when verifying the theoretical results numerically, we used the order of approximate computational convergence (COC) introduced by Jay [25], defined as: Another form of the numerical approximation of the order of theoretical convergence presented by the authors in [26], denoted as ACOC, is defined as: .
To show the order of convergence of the methods in the numerical tests, we used any of these estimates of computational order.In the tables, we write "-" when the COC (or ACOC) vector is unstable.Furthermore, the mean of elapsed time after 50 executions of the codes appears in the tables, calculated by using the CPU-time command.
Example 1.As a first example, we looked for the inverse random matrices of size n × n, where n = 10, 100, 200, 300, 400, and 500.Since Newton-Schulz's method needs one initial point, we took it as X 0 = A T /||A|| 2 .The rest of the methods need two initial points, which we took as X −1 = A T /||A|| 2 and X 0 = 0.5X −1 .
Table 1 shows the results obtained by approximating the inverses of nonsingular random matrices of sizes n = 10, 100, 200, 300, 400, and 500 using the Newton-Schulz, secant, Steffensen with memory, and Kurchatov-type methods.The number of iterations, residuals, and the COC value are shown.The results confirmed the theoretical order of convergence obtained from each method.Using all methods, the approximation of the inverse of A was obtained.In all cases, the Steffensen method with memory showed better results in terms of the number of iterations and the computational time.The graphs shown in Figure 2 represent the results presented in Table 1 for the cases of matrices of n = 10, 300, and 500.Graphs of the number of iterations vs. COC for three of the cases in Table 1.
Example 2. Now, we built square matrices of size n × n using different MATLAB functions, such as: (a) A = gallery( lehmer , n).Symmetric and positive definite matrix of size n × n, a i,j = i/j, ∀i, j.(b) A = gallery( riemann , n).Riemann matrix of size n × n.
A = gallery( ris , n).Hankel matrix of size n × n.(d) A = gallery( grcar , n).Toeplitz matrix of size n × n.(e) A = gallery( leslie , n).Leslie matrix of size n × n with application in problems of population models.(f) A = gallery( parter , n).Parter matrix of size n × n.
Here, we used as stopping criteria ||X k+1 − X k || 2 < 10 −10 or ||F(X k+1 )|| 2 < 10 −10 and the same initial approximations as in Example 1.The numerical results obtained are shown in Table 2 and Figure 3.As in the previous example, the proposed methods showed a good performance in terms of stability, precision, and number of iterations required.The results obtained for the number of iterations, the residuals, and the ACOC value of Example 3 are shown in Table 3 and Figure 4.The methods gave us an approximation of the inverse Moore-Penrose and showed the same behavior as in the previous examples.3.

Conclusions
In this manuscript, we widened the set of iterative methods able to be applied for estimating generalized inverses of complex matrices, by using schemes with memory able to improve the Newton-Schulz scheme.Two procedures with memory were designed for approximating the inverses of nonsingular complex matrices or pseudo-inverses, in the case that the matrices are singular.The order of convergence and stability were proven, and in the case of the Steffensen with memory scheme, its order of convergence improved that of the Newton-Schulz method.
The method using Kurchatov's divided differences cannot be directly adapted to estimate generalized inverses of complex matrices.To overcome this difficulty, Kurchatov-type divided differences were used.In this process, there a decrease in the order of convergence of the starting method and a change in its behavior.
The technique used to adapt the iterative methods can be applied to the resolution of other types of matrix equations with an especially significant role in various areas such as control theory, dynamic programming, ladder networks, statistics, etc.
This research opens new ways in the design of iterative procedures for solving this kind of nonlinear matrix equation, with promising numerical performance, which showed agreement with theoretical results and stable results.

Figure 2 .
Figure 2.Graphs of the number of iterations vs. COC for three of the cases in Table1.

Figure 3 .Example 3 .
Figure 3. Graphs of the number of iterations vs. COC for three of the sizes in Table 2.Example 3. Finally, we tested the methods for computing the Moore-Penrose inverse of m × n random matrices for different values of m and n.The matrices of the initial approximations were Matrix size 20 × 10 (b) Matrix size 300 × 400 (c) 1000 × 900

Figure 4 .
Figure 4. Graphs of the number of iterations vs. ACOC for three of the sizes in Table3.

Table 1 .
Results obtained by approximating the inverse of random matrices of different sizes.

Table 2 .
Results obtained by approximating the inverse of classical square matrices of different sizes.

Table 3 .
Results obtained to approximate the Moore-Penrose inverse of a rectangular random matrix.