Generalized Inverses Estimations by Means of Iterative Methods with Memory

Abstract

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.

1. 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\times n$. The extension of the iterative methods for the real equation $g\left(x\right)=ax-1=0$ to obtain the inverse of A, that is the zero of the matrix function $G\left(X\right)=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-A{X}_k),k=0,1,\dots ,$$

(1)

where I denotes the identity matrix of order n. Schulz [11] demonstrated that the eigenvalues of matrix $I-A{X}_0$ must be less than 1 to assure the convergence of the scheme in Equation (1). Taking into account that the residuals $E_k=I-A{X}_k$ in each iteration of Equation (1) satisfy $\parallel E_{k+1}\parallel \le \parallel E_k{\parallel }^2$, Newton–Schulz method has quadratic convergence. In general, it is known that this scheme converges to $A^{-1}$ with $X_0=\alpha A^{*}$ or $X_0=\alpha A$, where $0<\alpha <2/\rho \left(A^{*}A\right)$, $\rho (·)$ 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\in {\mathbb{C}}^{n\times n}$, Li et al. [12] suggested the scheme

$$X_{k+1}=X_k\left(mI-{\displaystyle \frac{m(m-1)}{2}}A{X}_k+\frac{m(m-1)(m-2)}{3!}-\cdots +{(-1)}^{m-1}{\left(A{X}_k\right)}^{m-1}\right),m=2,2,\dots ,$$

with $X_0=\alpha A^{*}$. They proved the convergence of m-order of $\left\{X_k\right\}$ 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 $\left\{X_k\right\}$ of estimations to the solution $\xi$, and let us also assume that this sequence converges to $\xi$. If there exists a nonzero constant $\eta$ and nonnegative numbers $t_i$, $0\le i\le m$ such that the inequality

$$|e_{k+1}|\le \eta \prod _{i=0}^m{\left|e_{k-i}\right|}^{t_i}$$

(2)

holds, where $e_k$ is the error of iterate $X_k$, then the R-order of convergence of (IM) satisfies

$$O_R(\left(IM\right),\xi )\ge s^{*},$$

(3)

where $s^{*}$ is the unique positive root of the polynomial

$$s^{m+1}-\sum _{i=0}^m{t}_i{s}^{m-i}=0.$$

(4)

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 Section 3 and Section 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.

2. A Secant-Type Method for Matrix Inversion

Let us recall that, for an scalar nonlinear equation $g\left(x\right)=0$, the secant method is an iterative scheme with memory such that

$$x_{k+1}=x_k-\frac{g\left(x_k\right)}{{\alpha }_k},$$

with ${\alpha }_k$ satisfying $g\left(x_k\right)-g\left(x_{k-1}\right)={\alpha }_k(x_k-x_{k-1})$, $k\ge 0$, given $x_0$ and $x_{-1}$ as initial approximations.

For a nonlinear matrix equation $G\left(X\right)=0$, where $G:{\mathbb{C}}^{n\times n}\to {\mathbb{C}}^{n\times n}$, the secant method can be described as

$$X_{k+1}=X_k-A_k^{-1}G\left(X_k\right),k\ge 0,$$

where $X_0$ and $X_{-1}$ are initial estimations and being $A_k$ a suitable linear operator satisfying

$$A_{k+1}(X_{k+1}-X_k)=G\left(X_{k+1}\right)-G\left(X_k\right)\iff A_{k+1}{S}_k=Y_k,$$

where $S_k=X_{k+1}-X_k$ and $Y_k=G\left(X_{k+1}\right)-G\left(X_k\right)$. 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\left(X\right)=0$.

Let us consider an $n\times 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

$$G\left(X\right)=X^{-1}-A=0.$$

(5)

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\left(X\right)=X^{-1}-A$ (see [15]) gives us:

$$\begin{array}{ccc}\hfill X_{k+1}& =& X_k-S_{k-1}{\left[G\left(X_k\right)-G\left(X_{k-1}\right)\right]}^{-1}G\left(X_k\right)\hfill \\ & =& X_k-(X_k-X_{k-1}){\left[X_k^{-1}-X_{k-1}^{-1}\right]}^{-1}(X_k^{-1}-A).\hfill \end{array}$$

(6)

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\times n$, then there exist unitary matrices U and V, of size $n\times n$, such that

$$\begin{array}{c}\hfill U^{*}AV=\Sigma =diag({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_n),\end{array}$$

(7)

being ${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_n>0$ the singular values of A.

Let us define $D_k=V^{*}{X}_{k}U$, that is $X_k=V{D}_k{U}^{*}$. Then, from Equation (6),

$$V{D}_{k+1}{U}^{*}=V{D}_k{U}^{*}-(V{D}_k{U}^{*}-V{D}_{k-1}{U}^{*}){(U{D}_k^{-1}{V}^{*}-U{D}_{k-1}^{-1}{V}^{*})}^{-1}(U{D}_k^{-1}{V}^{*}-U\Sigma V^{*}).$$

Several algebraic manipulations allow us to assure that

$$\begin{array}{c}\hfill D_{k+1}=D_k-(D_k-D_{k-1}){(D_k^{-1}-D_{k-1}^{-1})}^{-1}(D_k^{-1}-\Sigma ).\end{array}$$

(8)

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

$$D_{k+1}=D_{k-1}+D_k-D_{k-1}\Sigma D_k,$$

and, from this expression, we propose the secant-type method:

$$\begin{array}{c}\hfill X_{k+1}=X_{k-1}+X_k-X_{k-1}A{X}_k,k=0,1,2,\dots \end{array}$$

(9)

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\in {\mathbb{C}}^{nxn}$ be a nonsingular matrix, with singular value decomposition $U^{*}AV=\Sigma$. Let $X_0$ and $X_{-1}$ be such that $V^{*}{X}_{-1}U$ and $V^{*}{X}_{0}U$ are diagonal matrices. Then, sequence $\left\{X_k\right\}$, 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) is satisfied, where ${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_n>0$ are the singular values of A.

We define $D_k=V^{*}{X}_{k}U$, that is $X_k=V{D}_k{U}^{*}$, for $k\ge -1$. From Equation (9), we have

$$V{D}_{k+1}{U}^{*}=V{D}_{k-1}{U}^{*}+V{D}_k{U}^{*}-V{D}_{k-1}{U}^{*}U\Sigma V^{*}V{D}_k{U}^{*},$$

then

$$V{D}_{k+1}{U}^{*}=V(D_{k-1}+D_k-D_{k-1}\Sigma D_k)U^{*}$$

and therefore

$$D_{k+1}=D_{k-1}+D_k-D_{k-1}\Sigma D_k,$$

where $D_k=diag(d_k^1,d_k^2,\dots ,d_k^n)$.

Then, component by component, we obtain

$$\begin{array}{c}\hfill d_{k+1}^j=d_{k-1}^j+d_k^j-d_{k-1}^j{d}_k^j{\sigma }_j,j=1,2,\dots ,n.\end{array}$$

(10)

By subtracting ${\displaystyle \frac{1}{{\sigma }_j}}$ from both sides of Equation (10) and denoting $e_k^j=d_k^j-1/{\sigma }_j$, we get

$$\begin{array}{ccc}\hfill e_{k+1}^j& =& d_{k-1}^j+d_k^j-d_{k-1}^j{d}_k^j{\sigma }_j-{\displaystyle \frac{1}{{\sigma }_j}}\hfill \\ \hfill & =& -{\sigma }_j{e}_k^j{e}_{k-1}^j\hfill \end{array}$$

(11)

From Equation (11), we conclude that, for each value of j from 1 to n, $d_{k+1}^j$ in Equation (10) converges to ${\displaystyle \frac{1}{{\sigma }_j}}$ with order of convergence of the unique positive root of ${\lambda }^2-\lambda -1=0$, that is, $\lambda \approx 1.618$ (by using the result of Ortega–Rheinboldt mentioned in the Introduction).

Then, for each j, $1\le j\le n$, there exist a ${\{c_k^j\}}_k$ satisfying $c_k^j>0$, $\forall k$ and ${(c_k^j)}_k$ tends to zero when k tends to infinity. Moreover,

$$\left|e_{k+1}^j\right|\le c_k^j\left|e_k^j\right|,1\le j\le n.$$

Thus,

$${∥D_{k+1}-{\Sigma }^{-1}∥}_2^2=\sum _{j=1}^n{\left(e_{k+1}^j\right)}^2\le \sum _{j=1}^n{\left(c_k^j\right)}^2{\left(e_k^j\right)}^2\le n{m}_k^2{∥D_k-{\Sigma }^{-1}∥}_2^2,$$

where

$$m_k=\underset{1\le j\le n}{max}\left\{c_k^j\right\}.$$

Therefore,

$$\begin{array}{ccc}\hfill {∥X_{k+1}-A^{-1}∥}_2& =& {∥V{D}_{k+1}{U}^{*}-V{\Sigma }^{-1}{U}^{*}∥}_2\hfill \\ & =& {∥V(D_{k+1}-{\Sigma }^{-1})U∥}_2\hfill \\ & \le & {\parallel V\parallel }_2{∥(D_{k+1}-{\Sigma }^{-1})∥}_2{\parallel U^{*}\parallel }_2={∥(D_{k+1}-{\Sigma }^{-1})∥}_2\hfill \\ & \le & \sqrt{n}{m}_k{∥(D_k-{\Sigma }^{-1})∥}_2\hfill \\ & \le & \sqrt{n}{m}_k{∥(X_k-A^{-1})∥}_2,\hfill \end{array}$$

which allows us to affirm that $\left\{X_k\right\}$ 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\left(Z_k\right)$, 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^{\prime }\left(Z_{*}\right)$ has bounded powers, that is, there exists a positive constant C such that

$$\parallel H^{\prime }{\left(Z_{*}\right)}^k\parallel \le C,\forall k>0.$$

Therefore, the following result can be stated for the secant method.

Theorem 2.

The secant method in Equation (9) for the estimation of inverse matrix is a stable iterative scheme.

Proof. 

The proof is made demonstrating that $H^{\prime }\left(Z_{*}\right)$ is an idempotent matrix.

The secant-type method described as a fixed point scheme, can be written as

$$H\left(Z_k\right)=H\left(\begin{array}{c}{X}_k\\ X_{k-1}\end{array}\right)=\left(\begin{array}{c}{X}_{k-1}+X_k-X_{k-1}A{X}_k\\ X_k\end{array}\right).$$

It is easy to deduce that

$$H^{\prime }\left(\begin{array}{c}{X}_k\\ X_{k-1}\end{array}\right)Q=\left(\begin{array}{c}{Q}_1+Q_2-X_{k-1}A{Q}_1-Q_{2}A{X}_k\\ Q_1\end{array}\right),$$

where $Q={(Q_1,Q_2)}^T$. Then, for $Z=Z*={(A^{-1},A^{-1})}^T$, we have

$$H^{\prime }\left(Z_{*}\right)Q=\left(\begin{array}{c}0\\ Q_1\end{array}\right)=\left(\begin{array}{cc}0& 0\\ I& 0\end{array}\right)\left(\begin{array}{c}{Q}_1\\ Q_2\end{array}\right).$$

Thus, $H^{\prime }\left(Z_{*}\right)$ is an idempotent matrix and the iteration is stable. □

3. 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\times n$ complex matrix A, denoted by $A^{†}$. It is the unique $n\times m$ matrix X satisfying the equations

$$AXA=A,XAX=X,{\left(AX\right)}^{*}=AX,{\left(XA\right)}^{*}=XA.$$

If $rank\left(A\right)=r\le min\{m,n\}$, by using the singular value decomposition of A, we obtain

$$A=U\left[\begin{array}{cc}\Sigma & 0\\ 0& 0\end{array}\right]V^{*},$$

being $\Sigma =diag({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_r)$, ${\sigma }_1\ge {\sigma }_2,\dots ,\ge {\sigma }_r>0$. U and V are unitary matrices with $U\in {\mathbb{C}}^{m\times m}$ and $V\in {\mathbb{C}}^{n\times n}$. It is also known that

$$A^{†}=V^{*}\left[\begin{array}{cc}{\Sigma }^{-1}& 0\\ 0& 0\end{array}\right]U,$$

where ${\Sigma }^{-1}=diag(1/{\sigma }_1,1/{\sigma }_2,\dots ,1/{\sigma }_r)$.

The convergence of the method in Equation (9) for Moore–Penrose inverse is established in the following result.

Theorem 3.

Let $A\in {\mathbb{C}}^{m\times n}$ be a matrix with $rank\left(A\right)=r$, with singular value decomposition

$$U^{*}AV=\left(\begin{array}{cc}\Sigma & 0\\ 0& 0\end{array}\right).$$

Let $X_{-1}$ and $X_0$ be initial estimations such that

$$V^{*}{X}_{-1}U=\left(\begin{array}{cc}{\Sigma }_{-1}& 0\\ 0& 0\end{array}\right)and{V}^{*}{X}_{0}U=\left(\begin{array}{cc}{\Sigma }_0& 0\\ 0& 0\end{array}\right),$$

being ${\Sigma }_{-1}$ and ${\Sigma }_0$ diagonal matrices of size $r\times r$. Then, sequence $\left\{X_k\right\}$, 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

$$D_k=V^{*}{X}_{k}U=\left(\begin{array}{cc}{\Sigma }_k& 0\\ 0& 0\end{array}\right),$$

being ${\Sigma }_k\in {\mathbb{C}}^{r\times r}$. Thus, by using the iterative expression in Equation (9), we obtain

$$\left(\begin{array}{cc}{\Sigma }_{k+1}& 0\\ 0& 0\end{array}\right)=\left(\begin{array}{cc}{\Sigma }_{k-1}+{\Sigma }_k-{\Sigma }_{k-1}\Sigma {\Sigma }_k& 0\\ 0& 0\end{array}\right).$$

Therefore, as ${\Sigma }_{-1}$ and ${\Sigma }_0$ are diagonal matrices, so are all matrices ${\Sigma }_k$, and the expression

$${\Sigma }_{k+1}={\Sigma }_{k-1}+{\Sigma }_k-{\Sigma }_{k-1}\Sigma {\Sigma }_k$$

represents r scalar uncoupled iterations converging to ${\displaystyle \frac{1}{{\sigma }_i}}$, $1\le i\le r$ with super-linear order, that is to say,

$$\parallel {\Sigma }_{k+1}-{\Sigma }^{-1}{\parallel }_2^2\le r{M}_k^2{\parallel {\Sigma }_k-{\Sigma }^{-1}\parallel }_2^2,$$

with $M_k={max}_{1\le i\le r}\left\{c_k^2\right\}$, being $c_k^i>0$ such that sequence $\left\{c_k^i\right\}$ tends to zero for k tending to infinity.

With an analogous argument as in Theorem 1,

$$\parallel X_{k+1}-A^{†}{\parallel }_2\le \sqrt{r}{m}_k{\parallel X_k-A^{†}\parallel }_2,$$

which allows us to affirm that $\left\{X_k\right\}$ converges to $A^{†}$, with the desired order of convergence. □

4. 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\left(A^{l+1}\right)=rank\left(A^l\right)$ is called the index of A and it is denoted by $ind\left(A\right)$. If A is a complex matrix of size $n\times n$, the Drazin inverse of A, denoted $A^D$, is the unique matrix X satisfying

$$A^{l+1}X=A^l,XAX=X,\left(AX\right)=XA,$$

where l is the index of A.

If $ind\left(A\right)=1$, then X is called the g-inverse or group inverse of A, and, if $ind\left(A\right)=0$, then A is nonsingular and $A^D=A^{-1}$. Let us observe that the idempotent matrix $A{A}^D$ is the projector on $\mathcal{R}\left(A^l\right)$ along $\mathcal{N}\left(A^l\right)$, where $\mathcal{R}\left(A^l\right)$ and $\mathcal{N}\left(A^l\right)$ 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.

Proposition 1.

If $P_{A,B}$ is the projector on a space A along a space B, the following statements hold:

(a) 

$P_{A,B}C=C$ if and only if $\mathcal{R}\left(C\right)\subseteq A$.

(b) 

$C{P}_{A,B}=C$ if and only if $\mathcal{N}\left(C\right)\supseteq B$.

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=\alpha A^l$, where parameter $\alpha$ is chosen so that condition $\parallel I-A{X}_0\parallel <1$ is satisfied. One way for selecting the initial matrix used by different authors is

$$X_0=\frac{2}{tr\left(A^{l+1}\right)}{A}^l,$$

where $tr(·)$ is the trace of a square matrix. Another fruitful initial matrix is

$$X_0=\frac{2}{{2\parallel A\parallel }_2^{l+1}}{A}^l.$$

Using two initial matrices of these form, $\alpha A^l$, with $\alpha$ 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\in {\mathbb{C}}^{n\times n}$ be a square nonsingular matrix. We choose as initial estimations $X_0={\alpha }_0{A}^{l_0}$ and $X_1={\alpha }_1{A}^{l_1}$, with $l_0,l_1\ge ind\left(A\right)$. Then, sequence ${\left\{X_k\right\}}_{k\ge 0}$ generated by Equation (9) satisfies the following error equation

$$\parallel A^D-X_{k+1}\parallel \le \parallel A^D{\parallel \parallel A\parallel }^2\parallel A^D-X_{k-1}\parallel \parallel A^D-X_k\parallel .$$

Thus, ${\left\{X_k\right\}}_{k\ge 0}$ converges to $A^D$ with order of convergence $1.618$, that is, with super-linear convergence.

Proof. 

Let us define $E_k=I-A{X}_k$, $k=0,1,\dots$. Then,

$$\begin{array}{ccc}\hfill E_{k+1}& =& I-A{X}_{k+1}=I-A(X_k+X_{k-1}(I-A{X}_k))\hfill \\ & =& I-A{X}_k-A{X}_{k-1}(I-A{X}_k)\hfill \\ & =& (I-A{X}_{k-1})(I-A{X}_k)\hfill \\ & =& E_{k-1}{E}_k.\hfill \end{array}$$

Therefore, $\parallel E_{k+1}\parallel \le \parallel E_{k-1}\parallel \parallel E_k\parallel$. In addition, it is easy to prove that, if we choose $X_0$ and $X_1$ such that $\parallel E_0\parallel <1$ and $\parallel E_1\parallel <1$, then $\parallel E_k\parallel <1$, $\forall k\in \mathbb{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

$$A^{D}A{X}_k=X_k=X_{k}A{A}^D,\forall k\ge 0.$$

Thus,

$$e_k=A^D-X_k=A^D-A^{D}A{X}_k=A^D(I-A{X}_k)=A^D{E}_k.$$

From this identity, there exists $k_0\in \mathbb{N}$ such that

$$\parallel e_k\parallel \le \parallel A^D\parallel \parallel E_k\parallel \le \parallel A^D\parallel \parallel E_0{\parallel }^k{\parallel E_1\parallel }^k,\forall k\ge k_0.$$

Thus, $\{\parallel e_k{\parallel \}}_{k\ge 0}$ tends to zero and therefore ${\left\{X_k\right\}}_{k\ge 0}$ tends to $A^D$.

On the other hand,

$$\begin{array}{ccc}\hfill \parallel e_{k+1}\parallel & =& \parallel X_{k+1}-A^D\parallel =\parallel A^{D}A{X}_{k+1}-A^{D}A{A}^D\parallel \hfill \\ & =& \parallel A^D(A{X}_{k+1}-A{A}^D)\parallel \le \parallel A^D\parallel \parallel A{e}_{k+1}\parallel \hfill \end{array}$$

(12)

Now, we analyze $A{e}_{k+1}$.

$$A{e}_{k+1}=A(A^D-X_{k+1})=A{A}^D-I+I-A{X}_{k+1}=A{A}^D-I+E_{k+1}=A{A}^D-I+E_{k-1}{E}_k,$$

but

$$\begin{array}{ccc}\hfill E_{k-1}{E}_k+A{A}^D-I& =& (I-A{X}_{k-1})(I-A{X}_k)+A{A}^D-I\hfill \\ & =& (I-A{A}^D+A{A}^D-A{X}_{k-1})(I-A{A}^D+A{A}^D-A{X}_k)+A{A}^D-I\hfill \\ & =& (I-A{A}^D+A{e}_{k-1})(I-A{A}^D+A{e}_k)+A{A}^D-I\hfill \\ & =& {(I-A{A}^D)}^2+(I-A{A}^D)A{e}_k+A{e}_{k-1}(I-A{A}^D)+A{e}_{k-1}A{e}_k+A{A}^D-I\hfill \\ & =& A{e}_{k-1}A{e}_k.\hfill \end{array}$$

In the last equality, we use that ${(I-A{A}^D)}^2=I-A{A}^D$, in fact ${(I-A{A}^D)}^m=I-A{A}^D$, $\forall m\in \mathbb{N}$. In addition, $(I-A{A}^D)A{e}_k=0$ and $A{e}_{k-1}(I-A{A}^D)=0$.

Therefore,

$$\begin{array}{ccc}\hfill \parallel e_{k+1}\parallel & \le & \parallel A^D\parallel \parallel A{e}_{k+1}\parallel =\parallel A^D\parallel \parallel E_{k-1}{E}_k+A{A}^D-I\parallel \hfill \\ & \le & \parallel A^D{\parallel \parallel A\parallel }^2\parallel e_{k-1}\parallel \parallel e_k\parallel .\hfill \end{array}$$

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 ${\lambda }^2-\lambda -1=0$, that is $\lambda =1.618$. □

5. 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 $\parallel X_{k+1}-X_k{\parallel }_2<{10}^{-6}$ or $\parallel F\left(X_{k+1}\right){\parallel }_2<{10}^{-6}$.

To numerically check the theoretical results, Jay [17] introduced the order of approximate computational convergence (COC), defined as

$$\mathrm{order}\approx COC=\frac{ln(\parallel F\left(X_{k+1}\right){\parallel }_2/\parallel F\left(X_k\right){\parallel }_2)}{ln(\parallel F\left(X_k\right){\parallel }_2/\parallel F\left(X_{k-1}\right){\parallel }_2)}.$$

In a similar way, the authors presented in [18] another numerical approximation of the theoretical order, denoted by ACOC, and defined as

$$\mathrm{order}\approx ACOC=\frac{ln(\parallel X_{k+1}-X_k{\parallel }_2/\parallel X_k-X_{k-1}{\parallel }_2)}{ln(\parallel X_k-X_{k-1}{\parallel }_2/\parallel X_{k-1}-X_{k-2}{\parallel }_2)}.$$

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\times 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/{\parallel A\parallel }^2$ and for the secant method $X_{-1}={\displaystyle \frac{A^T}{{\parallel A\parallel }^2}}$ and $X_0=0.5{\displaystyle \frac{A^T}{{\parallel A\parallel }^2}}$.

In

Table 1. Results for approximating the inverse of a random matrix (Example 1).

Method	n	Iter	$\parallel {\mathit{X}}_{\mathit{k}+1}-{\mathit{X}}_{\mathit{k}}{\parallel }_2$	$\parallel \mathit{F}\left({\mathit{X}}_{\mathit{k}+1}\right){\parallel }_2$	COC
Newton–Schulz	10	19	$5.2\times {10}^{-7}$	$1.12\times {10}^{-14}$	2.0005
Secant	10	22	$5.4\times {10}^{-5}$	$9.8\times {10}^{-7}$	1.8660
Newton–Schulz	100	26	$2.0\times {10}^{-8}$	$5.4\times {10}^{-13}$	1.9988
Secant	100	36	$1.3\times {10}^{-6}$	$2.0\times {10}^{-7}$	1.6645
Newton–Schulz	200	32	$2.5\times {10}^{-12}$	$4.6\times {10}^{-12}$	2.0012
Secant	200	40	$1.6\times {10}^{-6}$	$1.8\times {10}^{-8}$	1.8866
Newton–Schulz	300	34	$3.1\times {10}^{-12}$	$5.9\times {10}^{-12}$	1.8888
Secant	300	40	$3.7\times {10}^{-5}$	$3.6\times {10}^{-7}$	1.8865
Newton–Schulz	400	36	$2.2\times {10}^{-10}$	$1.9\times {10}^{-11}$	2.0001
Secant	400	43	$3.5\times {10}^{-7}$	$1.1\times {10}^{-8}$	1.8222
Newton–Schulz	500	33	$1.5\times {10}^{-7}$	$1.2\times {10}^{-11}$	1.9999
Secant	500	36	$9.0\times {10}^{-5}$	$2.6\times {10}^{-7}$	1.6666

, 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\times 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. Results for approximating the Moore–Penrose inverse of a rectangular random matrix (Example 2).

Method	m	n	Iter	$\parallel {\mathit{X}}_{\mathit{k}+1}-{\mathit{X}}_{\mathit{k}}{\parallel }_2$	ACOC
Newton–Schulz	20	10	14	$9.7\times {10}^{-12}$	2.0005
Secant	20	10	13	$9.9\times {10}^{-10}$	1.6199
Newton–Schulz	200	100	17	$4.1\times {10}^{-10}$	2.0018
Secant	200	100	17	$2.02\times {10}^{-9}$	1.6210
Newton–Schulz	300	400	21	$2.03\times {10}^{-11}$	2.0007
Secant	300	400	27	$1.4\times {10}^{-7}$	1.6267
Newton–Schulz	500	600	23	$3.8\times {10}^{-9}$	2.0028
Secant	500	600	31	$5.2\times {10}^{-10}$	1.6197
Newton–Schulz	1000	900	25	$4.5\times {10}^{-8}$	2.0055
Secant	1000	900	36	$2.5\times {10}^{-9}$	1.6205

, 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\times 6$ with $ind\left(A\right)=2$.

$$A=\left(\begin{array}{cccccc}\hfill 1& \hfill -1& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill -1& \hfill 1& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill -1& \hfill -1& \hfill 1& \hfill -1& \hfill 0& \hfill 0\\ \hfill -1& \hfill -1& \hfill -1& \hfill 1& \hfill 0& \hfill 0\\ \hfill -1& \hfill -1& \hfill -1& \hfill 0& \hfill 2& \hfill -1\\ \hfill -1& \hfill -1& \hfill 0& \hfill -1& \hfill -1& \hfill 2\end{array}\right)$$

Here, its Drazin inverse is expressed by

$$A^D=\left(\begin{array}{cccccc}1/4& -1/4& 0& 0& 0& 0\\ -1/4& 1/4& 0& 0& 0& 0\\ 0& 0& 1/4& -1/4& 0& 0\\ 0& 0& -1/4& 1/4& 0& 0\\ 0& 0& -5/12& -7/12& 2/3& 1/3\\ 0& 0& -7/12& -5/12& 1/3& 2/3\end{array}\right).$$

By using the initial matrix $X_0={\displaystyle \frac{0.5}{tr\left(A^3\right)}}$ and the same stopping criterion as in the previous examples, Newton–Schulz method gives us the following information:

• $ACOC=2.0009$;
• $iter=11$; and
• Exact error $\parallel A^D-X_{11}{\parallel }_2=7.7716\times {10}^{-16}$.

On the other hand, secant method is used with $X_{-1}={\displaystyle \frac{1}{tr\left(A^3\right)}}$ and $X_0={\displaystyle \frac{0.5}{tr\left(A^3\right)}}$, obtaining:

• $ACOC=1.6225$;
• $iter=15$; and
• Exact error $\parallel A^D-X_{15}{\parallel }_2=1.8539\times {10}^{-13}$.

Example 4.

This is another example for computing the Drazin inverse of the following matrix B of size $12\times 12$ (see [1]) with $ind\left(B\right)=3$.

$$B=\left(\begin{array}{cccccccccccc}\hfill 2& \hfill 0.4& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill -2& \hfill 0.4& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill -1& \hfill -1& \hfill 1& \hfill -1& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill -1& \hfill -1& \hfill -1& \hfill 1& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 1& \hfill 1& \hfill -1& \hfill -1& \hfill 0& \hfill 0& \hfill -1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 1& \hfill 1& \hfill -1& \hfill -1& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill -1& \hfill -2& \hfill 0.4& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 2& \hfill 0.4& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill -1& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 1& \hfill -1& \hfill -1& \hfill -1\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill -1& \hfill 1& \hfill -1& \hfill -1\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0.4& \hfill -2\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0.4.4& \hfill 2\end{array}\right).$$

Now, its Drazin inverse is expressed by

$$A^D=\left(\begin{array}{cccccccccccc}0.25& -0.25& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 1.25& 1.25& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ -1.6641& -0.9922& 0.25& -0.25& 0& 0& 0& 0& -0.0625& -0.0625& 0& 0.1563\\ -1.1953& -0.6797& -0.25& 0.25& 0& 0& 0& 0& -0.0625& 0.1875& 0.6875& 1.3438\\ -2.7637& -1.0449& -1.875& -1.25& -1.25& -1.25& -1.25& -1.25& 1.4844& 2.5781& 3.3203& 6.6406\\ -2.7637& -1.0449& -1.875& -1.25& -1.25& -1.25& -1.25& -1.25& 1.4844& 2.5781& 4.5703& 8.5156\\ 14.1094& 6.3008& 6.625& 3.375& 5& -3& -5& -5& -4.1875& -8.5& -10.5078& -22.4609\\ -19.3242& -8.5078& -9.75& -5.25& -7.5& 4.5& 7.5& 7.5& 6.375& 12.5625& 15.9766& 33.7891\\ -0.625& -0.3125& 0& 0& 0& 0& 0& 0& 0.25& -0.25& -0.875& -1.625\\ -1.25& -0.9375& 0& 0& 0& 0& 0& 0& -0.25& 0.25& -0.875& -1.625\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1.25& 1.25\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& -0.25& 0.25\end{array}\right).$$

By using the initial matrix $X_0={\displaystyle \frac{0.5}{tr\left(A^5\right)}}$ 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 $\parallel B^D-X_{14}\parallel =1.8354\times {10}^{-9}$.

On the other hand, secant method is used with $X_{-1}={\displaystyle \frac{1}{tr\left(A^5\right)}}$ and $X_0={\displaystyle \frac{0.5}{tr\left(A^5\right)}}$, obtaining:

• $ACOC=1.6201$;
• $iter=20$; and
• Exact error $\parallel B^D-X_{20}\parallel =1.8453\times {10}^{-9}$.

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\times n$, constructed by using different Matlab functions. Specifically, the test matrices are:

(a) 

$A=gallery{(}^{\prime }ri{s}^{\prime },n)$. Hankel matrix of size n × n.

(b) 

$A=gallery{(}^{\prime }grca{r}^{\prime },n)$. Toeplitz matrix of size n × n.

(c) 

$A=gallery{(}^{\prime }lehme{r}^{\prime },n)$. Symmetric and positive definite matrix of size n × n, $a_{i,j}=i/j,\forall i,j$.

(d) 

$A=gallery{(}^{\prime }lesli{e}^{\prime },n)$. Leslie matrix of size n × n. This type of matrices appears in problems of population models.

(e) 

$A=gallery{(}^{\prime }inv{o}^{\prime },n)$. Matrix ill-conditioned of size n × n, such that $A^2=I$.

By using the stopping criterion

$$\parallel X_{k+1}-X_k{\parallel }_2<{10}^{-10}or{\parallel F\left(X_{k+1}\right)\parallel }_2<{10}^{-10}$$

and the initial matrix $X_{-1}=\frac{A^T}{{\parallel A\parallel }^2}$ and $X_0=0.5\frac{A^T}{{\parallel A\parallel }^2}$, we obtain the numerical results that appear in

Table 3. Results for approximating the inverse of classical square matrices (Example 5).

Method	Matrix	n	Iter	$\parallel {\mathit{X}}_{\mathit{k}+1}-{\mathit{X}}_{\mathit{k}}{\parallel }_2$	$\parallel \mathit{F}\left({\mathit{X}}_{\mathit{k}+1}\right)\parallel$	COC
Newton–Schulz	Lehmer	10	18	$3.5\times {10}^{-7}$	$6.3\times {10}^{-15}$	-
Secant	Lehmer	10	20	$3.9\times {10}^{-9}$	$1.7\times {10}^{-11}$	1.6164
Newton–Schulz	Hankel	100	8	$1.1\times {10}^{-5}$	$1.2\times {10}^{-11}$	1.9993
Secant	Hankel	100	11	$1.9\times {10}^{-12}$	$4.4\times {10}^{-13}$	1.6180
Newton–Schulz	Toeplitz	200	9	$1.6\times {10}^{-9}$	$3.2\times {10}^{-15}$	1.9976
Secant	Toeplitz	200	11	$6.3\times {10}^{-11}$	$6.4\times {10}^{-11}$	1.6182
Newton–Schulz	Toeplitz	300	9	$1.7\times {10}^{-9}$	$2.5\times {10}^{-15}$	1.9975
Secant	Toeplitz	300	11	$6.3\times {10}^{-11}$	$6.4\times {10}^{-11}$	1.6182
Newton–Schulz	Leslie	400	22	$4.3\times {10}^{-5}$	$2.3\times {10}^{-13}$	1.9995
Secant	Leslie	400	33	$4.2\times {10}^{-12}$	$1.0\times {10}^{-14}$	1.6177
Newton–Schulz	Leslie	500	23	$1.6\times {10}^{-6}$	$1.4\times {10}^{-16}$	2.0001
Secant	Leslie	500	25	$1.7\times {10}^{-12}$	$3.8\times {10}^{-15}$	1.6070

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

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