After notes on Chebyshev’s iterative method

In this case, F″(x) is a constant bilinear operator that we denote by B.

Using Taylor expansions,

F(yn)=F(xn)+F'(xn)(yn−xn)+12F''(xn)(yn−xn)2=12B(yn−xn)2.$$\begin{array}{} \displaystyle F(y_n)=F(x_n)+F^{'}(x_n)(y_n-x_n)+\frac{1}{2}F^{''}(x_n)(y_n-x_n)^2 =\frac{1}{2} B(y_n-x_n)^2. \end{array}$$(8)

Thus, F'(xn)−1F(yn)=−12LF(xn)(yn−xn)$\begin{array}{} \displaystyle F^{'}(x_n)^{-1} F(y_n)=-\frac{1}{2}L_{F}(x_n)(y_n-x_n) \end{array}$, and the method becomes

yn=xn−F'(xn)−1F(xn),xn+1=yn−F'(xn)−1F(yn),$$\begin{array}{} \displaystyle \begin{array}{*{20}{r}} {{y_n}}& = &{{x_n} - F'{{({x_n})}^{ - 1}}F({x_n}),}\\ {{x_{n + 1}}}& = &{{y_n} - F'{{({x_n})}^{ - 1}}F({y_n}),} \end{array} \end{array}$$

equivalently,

F'(xn)(yn−xn)=−F(xn),F'(xn)(xn+1−yn)=−F(yn).$$\begin{array}{} \displaystyle \begin{array}{*{20}{r}} {F'\left( {{x_n}} \right)\left( {{y_n} - {x_n}} \right)}& = &{ - F({x_n}),}\\ {F'\left( {{x_n}} \right)\left( {{x_{n + 1}} - {y_n}} \right)}& = &{ - F({y_n}).} \end{array} \end{array}$$(9)

Notice that we only need a LU decomposition and three evaluations (F(x_n), F(y_n) and F′(x_n)) in each iteration. In particular, for this problem Chebyshev’s method is more efficient than Newton’s method.

Moreover, we can obtain a very simple semilocal convergence:

Approximating inverse operators is a very common task in several areas of interest, such as physics, chemistry, engineering, etc. In a general context, we can formulate the following problem: given g, we are interested in calculating f ε Ω such that H(f) = g, where H : Ω → Y is an operator defined in a domain Ω of a Banach space X with values in a Banach space Y . It is clear that we have to calculate or approximate the inverse operator H⁻¹ for solving the previous equation. In this case, if g is in the domain of H⁻¹, there is a solution f = H⁻¹(g).

To approximate the inverse operator H⁻¹, we use Newton-type methods, so that better successive approximations to H⁻¹ are constructed from an initial approximation. The methods considered in this paper, when they are applied to compute an inverse operator, are not be based on solving linear systems, if not on the product of matrices. Notice that the formulation of the problem in this manner is very interesting.

Let X and Y be two Banach spaces and GL(X,Y)={H∈ℒ(X,Y):H−1exists}$\begin{array}{} \displaystyle GL(X,Y)=\{H\in\mathcal{L}(X,Y):\, H^{-1}\,\text{exists}\} \end{array}$, where ℒ(X,Y)$\begin{array}{} \displaystyle \mathcal{L}(X,Y) \end{array}$ is the set of bounded linear operators from the Banach space X into the Banach space Y . The problem that we are thinking about is the following: given an operator H ε GL(X,Y), we approximate H⁻¹. To do this, we first consider

ℱ:GL(Y,X)→ℒ(X,Y)andℱ(G)=G−1−H,$$\begin{array}{} \displaystyle \mathcal{F}:GL(Y,X)\rightarrow\mathcal{L}(X,Y) \quad\text{and}\quad \mathcal{F}(G)=G^{-1}-H, \end{array}$$

so that H⁻¹ is the solution of the equation ℱ(G)=0$\begin{array}{} \displaystyle \mathcal{F}(G)=0 \end{array}$.

If we observe Newton’s method,

{G0given,Gn+1=Gn−[ℱ′(Gn)]−1ℱ(Gn),n≥0,$$\begin{array}{} \displaystyle \left\{ \begin{array}{*{20}{l}} {{G_0}\;{\rm{given}},}\\ {{G_{n + 1}} = {G_n} - {{[{\cal F}'({G_n})]}^{ - 1}}{\cal F}({G_n}),\quad n \ge 0,} \end{array}\right. \end{array}$$

it is clear that inverse operators are used, but if we take into account the definition of ℱ$\begin{array}{} \displaystyle \mathcal{F} \end{array}$, calculate ℱ′(Gn)$\begin{array}{} \displaystyle \mathcal{F}'(G_n) \end{array}$ and do

ℱ′(Gn)(Gn+1−Gn)=−ℱ(Gn),n≥0,$$\begin{array}{} \displaystyle \mathcal{F}'(G_n)(G_{n+1}-G_n)=-\mathcal{F}(G_n),\quad n\geq 0, \end{array}$$(12)

then we can avoid the use of inverse operators for approximating G_n+1.

Indeed, to obtain the corresponding algorithm, we only need to compute ℱ′(Gn)$\begin{array}{} \displaystyle \mathcal{F}'(G_n) \end{array}$. So, given G ε GL(Y,X), as G⁻¹ exists, if

0<ε<1‖α‖‖G−1‖,$$\begin{array}{} \displaystyle 0 < \varepsilon < \dfrac{1}{\|\alpha\|\|G^{-1}\|}, \end{array}$$

we have ‖εα‖<1‖G−1‖$\begin{array}{} \displaystyle \|\varepsilon\alpha\| < \frac{1}{\|G^{-1}\|} \end{array}$ for α ε GL(Y,X). Therefore, it is known that G + εα ε GL(Y,X) and then

ℱ′(G)α=limε→01ε[ℱ(G+εα)−ℱ(G)]=−G−1αG−1.$$\begin{array}{} \displaystyle \mathcal{F}'(G)\alpha= \lim_{\varepsilon\rightarrow 0}\frac{1}{\varepsilon} [\mathcal{F}(G+\varepsilon\alpha)-\mathcal{F}(G)]= -G^{-1}\alpha G^{-1}. \end{array}$$

In consequence, Newton’s method is now given by the following algorithm:

{G0given,Gn+1=2Gn−GnHGn,n≥0.$$\begin{array}{} \displaystyle \left\{ \begin{array}{*{20}{l}} {{G_0}\;{\rm{given}},}\\ {{G_{n + 1}} = 2{G_n} - {G_n}H{G_n},\quad n \ge 0.} \end{array}\right. \end{array}$$(13)

In addition, Newton’s method does not use inverse operators for approximating an inverse operator and the order of convergence is two.

If we now consider Chebyshev’s method,

{G0given,Gn+1=Gn−[I+12Lℱ(Gn)][ℱ′(Gn)]−1ℱ(Gn),n≥0,$$\begin{array}{} \displaystyle \left\{ \begin{array}{*{20}{l}} {{G_0}\;{\rm{given}},}\\ {{G_{n + 1}} = {G_n} - \left[ {I + \frac{1}{2}{L_{\cal F}}({G_n})} \right]{{[{\cal F}'({G_n})]}^{ - 1}}{\cal F}({G_n}),\quad n \ge 0,} \end{array}\right. \end{array}$$

where Lℱ(Gn)=[ℱ′(Gn)]−1ℱ″(Gn)[ℱ′(Gn)]−1ℱ(Gn)$\begin{array}{} \displaystyle L_{\mathcal{F}}(G_n)=[\mathcal{F}'(G_n)]^{-1}\mathcal{F}''(G_n)[\mathcal{F}'(G_n)]^{-1}\mathcal{F}(G_n) \end{array}$, we can think that inverse operators are used, but we can do the same as in Newton’s method to see that Chebyshev’s method does not use them:

{ℱ′(Gn)(Pn−Gn)=−ℱ(Gn),n≥0,ℱ′(Gn)(Gn+1−Pn)=−12ℱ″(Gn)(Pn−Gn)2,$$\begin{array}{} \displaystyle \left\{ \begin{array}{*{20}{l}} {{\cal F}'({G_n})({P_n} - {G_n}) = - {\cal F}({G_n}),\quad n \ge 0,}\\ {{\cal F}'({G_n})({G_{n + 1}} - {P_n}) = - \frac{1}{2}{\cal F}''({G_n}){{({P_n} - {G_n})}^2},} \end{array}\right. \end{array}$$(14)

so that we can also avoid the use of inverse operators for approximating G_n+1.

Let α,β ε GL(Y,X) and then

0<ε<1‖β‖‖G−1‖,$$\begin{array}{} \displaystyle 0 < \varepsilon < \dfrac{1}{\|\beta\|\|G^{-1}\|}, \end{array}$$

so that G + εβ ε GL(Y,X) and

ℱ″(G)αβ=limε→01ε[ℱ′(G+εβ)α−ℱ′(G)α]=G−1αG−1βG−1+G−1βG−1αG−1.$$\begin{array}{} \displaystyle \mathcal{F}''(G)\alpha\beta= \lim_{\varepsilon\rightarrow 0}\frac{1}{\varepsilon} [\mathcal{F}'(G+\varepsilon\beta)\alpha-\mathcal{F}'(G)\alpha]= G^{-1}\alpha G^{-1}\beta G^{-1} + G^{-1}\beta G^{-1}\alpha G^{-1}. \end{array}$$

In consequence, we write Chebyshev’s method as

{G0given,Gn+1=3Gn−3GnHGn+GnHGnHGn,n≥0.$$\begin{array}{} \displaystyle \left\{ \begin{array}{*{20}{l}} {{G_0}\;{\rm{given}},}\\ {{G_{n + 1}} = 3{G_n} - 3{G_n}H{G_n} + {G_n}H{G_n}H{G_n},\quad n \ge 0.} \end{array}\right. \end{array}$$(15)

Observe that Chebyshev’s method does not use inverse operators for approximating an inverse operator and the order of convergence is three.

If we consider the complex equation f(x) = x^p − a with a∈ℂ$\begin{array}{} \displaystyle a\in\mathbb{C} \end{array}$, then Lf(x)=(p−1p)xp−axp$\begin{array}{} \displaystyle L_{f}(x)=\left(\dfrac{p-1}{p}\right)\dfrac{x^p-a}{x^p} \end{array}$ and Chebyshev’s method is reduced to

x0∈D,xn+1=2p2−3p+12p2xn+2p−1p2axn1−p−p−12p2a2xn1−2p,n≥0.$$\begin{array}{} \displaystyle {x_0} \in D,{\rm{ }}{x_{n + 1}} = \frac{{2{p^2} - 3p + 1}}{{2{p^2}}}{\mkern 1mu} {x_n} + \frac{{2p - 1}}{{{p^2}}}{\mkern 1mu} a{\mkern 1mu} x_n^{1 - p} - \frac{{p - 1}}{{2{p^2}}}{\mkern 1mu} {a^2}{\mkern 1mu} x_n^{1 - 2p},\quad n \ge 0. \end{array}$$

If we extend Chebyshev’s method for the computation of the pth root of a matrix, we then consider the space

Θ={A∈ℂr×r s.t. A has no nonpositive real eigenvalues}$$\begin{array}{} \displaystyle \Theta = \left\{A \in \mathbb{C}^{r\times r}\,\text{ s.t. } A \text{ has no nonpositive real eigenvalues} \right\} \end{array}$$

and approximate A1p$\begin{array}{} \displaystyle A^{\frac{1}{p}} \end{array}$ for a given matrix X ε Θ. For this, we consider

ℱ:Θ→ℂr×randℱ(X)=Xp−A,$$\begin{array}{} \displaystyle \mathcal{F}:\Theta\rightarrow\mathbb{C}^{r \times r} \quad\text{and}\quad \mathcal{F}(X)=X^{p}-A, \end{array}$$

so that A1p$\begin{array}{} \displaystyle A^{\frac{1}{p}} \end{array}$ is the solution of the equation ℱ(X)=0$\begin{array}{} \displaystyle \mathcal{F}(X)=0 \end{array}$. So, Chebyshev’s method is reduced to

X0∈Θ,Xn+1=2p2−3p+12p2Xn+2p−1p2AXn1−p−p−12p2A2Xn1−2p,n≥0.$$\begin{array}{} \displaystyle X_0\in\Theta,\quad X_{n+1} = \dfrac{2p^2-3p+1}{2p^2}\,X_n + \dfrac{2p-1}{p^2}\, A X_n^{1-p} - \dfrac{p-1}{2p^2}\,A^2 X_n^{1-2p},\quad n\geq 0. \end{array}$$(16)

Chebyshev’s method, as Newton’s method, cannot be used directly to approximate the principal pth root. Using the same idea as in [15], one can prove that it is not stable in a neighborhood of A1p$\begin{array}{} \displaystyle A^{\frac{1}{p}} \end{array}$. A small perturbation on the value of X_n is amplified in the following steps and in a finite arithmetic computation the algorithm diverges. This problem can be overcome using another algorithm which provides the same sequence but it is stable in a neighborhood of A1p$\begin{array}{} \displaystyle A^{\frac{1}{p}} \end{array}$.

The following iteration, given by Denman and Beavers in [9],

{X0=A,Y0=I,Xn+1=12(Xn+Yn−1),Nn+1=12(Yn+Xn−1),n≥0.$$\begin{array}{} \displaystyle \left\{ \begin{array}{*{20}{l}} {{X_0} = A,\quad {Y_0} = I,}\\ {{X_{n + 1}} = \frac{1}{2}({X_n} + Y_n^{ - 1}),}\\ {{N_{n + 1}} = \frac{1}{2}({Y_n} + X_n^{ - 1}),\quad n \ge 0.} \end{array}\right. \end{array}$$

is a stable variant of the Newton method for the matrix square root.

The instability of the simplified Newton iterations Xn+1=(Xn+AXn−1)/2$\begin{array}{} \displaystyle X_{n+1}=(X_n+AX_n^{-1})/2 \end{array}$ and Xn+1=(Xn+Xn−1A)/2$\begin{array}{} \displaystyle X_{n+1}=(X_n+X_n^{-1} A)/2 \end{array}$, shown by Higham [15] for the matrix square root, is mainly due to the one-sided multiplication by A. On the other hand, since X_n and A commute if AX₀ = X₀A, the iteration can be rewritten as

Xn+1=Xn+A12Xn−1A122,n≥0,$$\begin{array}{} \displaystyle X_{n+1}=\frac{X_n+ A^\frac{1}{2} X_n^{-1} A^\frac{1}{2}}{2}, \quad n\geq 0, \end{array}$$

and the iteration becomes stable in this form. However, it is useless since it involves the square root of A, but it helps us to stabilize the iteration by introducing the variable Yn=A12Xn−1A12=A−1Xn=XnA−1$\begin{array}{} \displaystyle Y_n= A^\frac{1}{2} X_n^{-1} A^\frac{1}{2}=A^{-1}X_n=X_nA^{-1} \end{array}$. The resulting iteration is that of Denman and Beavers, see [16].

Following these ideas, Iannazo [17] proposes the following two stable versions of the Newton method for the matrix pth root:

Xn+1=(p−1)Xn+(A1pXn−1)p−1A1pp,n≥0,$$\begin{array}{} \displaystyle {X_{n + 1}} = \frac{{(p - 1){X_n} + {{\left( {{A^{\frac{1}{p}}}X_n^{ - 1}} \right)}^{p - 1}}{A^{\frac{1}{p}}}}}{p},{\rm{ }}n \ge 0, \end{array}$$

and

{X0=I,N0=A,Xn+1=Xn((p−1)I+Nnp),Nn+1=((p−1)I+Nnp)−pNn,n≥0.$$\begin{array}{} \displaystyle \left\{ \begin{array}{*{20}{l}} {{X_0} = I,\quad {N_0} = A,}\\ {{X_{n + 1}} = {X_n}\left( {\frac{{(p - 1)I + {N_n}}}{p}} \right),}\\ {{N_{n + 1}} = {{\left( {\frac{{(p - 1)I + {N_n}}}{p}} \right)}^{ - p}}{N_n},\quad n \ge 0.} \end{array}\right. \end{array}$$

Observe that in the second case, the matrix A does not explicitly appear in the iteration.

Similarly, for Chebyshev’s method, we propose

Xn+1=(p−32+12p)Xn+(2−1p)(A1pXn−1)p−1A1pp−(p−1)(A1pXn−1)2p−1A1p2p2,n≥0.$$\begin{array}{} \displaystyle X_{n+1}=\frac{\left(p-\frac{3}{2}+\frac{1}{2p}\right) X_n + \left(2-\frac{1}{p}\right) \left(A^{\frac{1}{p}} X_n^{-1} \right)^{p-1} A^{\frac{1}{p}}}{p} - \frac{ (p-1) \left(A^{\frac{1}{p}} X_n^{-1} \right)^{2p-1} A^{\frac{1}{p}}}{2 p^2},\quad n\geq 0. \end{array}$$(17)

and

{X0=I,N0=A,Xn+1=Xn((p−32+12p)I+(2−1p)Nnp−(p−1)Nn22p2),n≥0,Nn+1=((p−32+12p)I+(2−1p)Nnp−(p−1)Nn22p2)−pNn.$$\begin{array}{} \displaystyle \left\{ \begin{array}{*{20}{l}} {{X_0} = I,\quad {N_0} = A,}\\ {{X_{n + 1}} = {X_n}\left( {\frac{{\left( {p - \frac{3}{2} + \frac{1}{{2p}}} \right)I + \left( {2 - \frac{1}{p}} \right){N_n}}}{p} - \frac{{(p - 1)N_n^2}}{{2{p^2}}}} \right),\quad n \ge 0,}\\ {{N_{n + 1}} = {{\left( {\frac{{\left( {p - \frac{3}{2} + \frac{1}{{2p}}} \right)I + \left( {2 - \frac{1}{p}} \right){N_n}}}{p} - \frac{{(p - 1)N_n^2}}{{2{p^2}}}} \right)}^{ - p}}{N_n}.} \end{array}\right. \end{array}$$(18)

Note that (17) is useless, since it involves the matrix pth root of A.

Observe that we need to calculate Xn−1$\begin{array}{} \displaystyle X_{n}^{-1} \end{array}$ at each step of Chebychev’s method given by (16). We are interesting in obtaining other expression of Chebyshev’s method which avoids the computation of these inverses. So, if we consider the complex function f(x)=1xp−1a=0$\begin{array}{} \displaystyle f(x)=\frac{1}{x^{p}}-\frac{1}{a}=0 \end{array}$, where f:D⊆ℂ→ℂ$\begin{array}{} \displaystyle f: D\subseteq\mathbb{C}\rightarrow\mathbb{C} \end{array}$, then algorithm (17) is reduced to

x0∈D,xn+1=xn(1+1p(1−xnpa)+p+12p2(1−xnpa)2),n≥0.$$\begin{array}{} \displaystyle x_0\in D,\qquad x_{n+1} = x_n\left(1+\dfrac{1}{p}\left(1-\dfrac{x_n^p}{a}\right) + \dfrac{p+1}{2p^2}\left(1-\dfrac{x_n^p}{a}\right)^2\right),\quad n\geq 0. \end{array}$$(19)

As a consequence, the algorithm of Chebyshev’s method for solving ℱ(X)=0$\begin{array}{} \displaystyle \mathcal{F}(X)=0 \end{array}$, where ℱ:Θ→ℂr×r$\begin{array}{} \displaystyle \mathcal{F}:\Theta\rightarrow\mathbb{C}^{r \times r} \end{array}$, is then

X0∈Θ,Xn+1=Xn(I+1p(I−A−1Xnp)+p+12p2(I−A−1Xnp)2),n≥0.$$\begin{array}{} \displaystyle X_0\in\Theta,\quad X_{n+1} = X_n\left(I+\dfrac{1}{p}\left(I-A^{-1}X_n^p\right) + \dfrac{p+1}{2p^2}\left(I-A^{-1}X_n^p\right)^2\right),\quad n\geq 0. \end{array}$$(20)

Observe that we only need to calculate A⁻¹ and not Xn−1$\begin{array}{} \displaystyle X_{n}^{-1} \end{array}$ at each step, so that algorithm (20) is more efficient than algorithm (16).

For Halley’s method, a similar strategy in order to avoid the use of inverses is not possible. Remember that the algorithm of Halley’s method,

x0∈D,xn+1=xn−(11+12Lf(xn))f(xn)f'(xn),n≥0,$$\begin{array}{} \displaystyle {x_0} \in D,\qquad {x_{n + 1}} = {x_n} - \left( {\frac{1}{{1 + \frac{1}{2}{L_f}({x_n})}}} \right)\frac{{f({x_n})}}{{f'({x_n})}},\quad n \ge 0, \end{array}$$

always involves inverses (the quotient includes L_f (x_n)). This is the main advantage of Chebyshev’s method with respect to Halley’s method.

After that, we analyze the convergence of Chebyshev’s method. For this, we suppose AX₀ = X₀A and consider the following residual of algorithm (20):

R(Xn)=I−A−1Xnp,n≥0.$$\begin{array}{} \displaystyle R(X_n)=I-A^{-1}X_n^p,\quad n\geq 0. \end{array}$$

To prove the convergence of algorithm (20), we consider a submultiplicative matrix norm || ⋅ || defined in ℂr×r$\begin{array}{} \displaystyle \mathbb{C}^{r \times r} \end{array}$ and prove that {||R(X_n)|| is a scalar decreasing sequence convergent to zero. First of all, from AX₀ = X₀A, it follows AX_n = X_nA, AXnp=XnpA$\begin{array}{} \displaystyle AX_n^{p}=X_n^{p}A \end{array}$ and A−1Xnp=XnpA−1$\begin{array}{} \displaystyle A^{-1}X_n^{p}=X_n^{p}A^{-1} \end{array}$. As a consequence, we have

R(Xn)=I−(I−R(Xn−1))(I+1pR(Xn−1)+p+12p2R(Xn−1)2)p.$$\begin{array}{} \displaystyle R({X_n}) = I - (I - R({X_{n - 1}})){\left( {I + \frac{1}{p}R({X_{n - 1}}) + \frac{{p + 1}}{{2{p^2}}}R{{({X_{n - 1}})}^2}} \right)^p}. \end{array}$$

Next, taking into account Villareal’s formula,

P(y)p=(∑i=0maiyi)p=∑i=0mpPiyi,$$\begin{array}{} \displaystyle P(y)^p=\left(\sum_{i=0}^m a_{i}y^{i}\right)^p=\sum_{i=0}^{mp}P_{i}y^{i}, \end{array}$$

where P0=a0p$\begin{array}{} \displaystyle P_{0}=a_{0}^p \end{array}$, Pi=∑j=0i−1Pjai−ja0(i−j)(p+1)−ii$\begin{array}{} \displaystyle {P_{i}=\sum_{j=0}^{i-1}P_{j}\,\frac{a_{i-j}}{a_{0}}\,\frac{(i-j)(p+1)-i}{i}} \end{array}$, for i = 1,2,...,mp, and a_i = 0, for all i ≥ m + 1, it follows that

R(Xn)=I−(I−R(Xn−1))∑i=02pPiR(Xn−1)i,$$\begin{array}{} \displaystyle R(X_n) = I-(I-R(X_{n-1}))\sum_{i=0}^{2p} P_{i}R(X_{n-1})^{i}, \end{array}$$(21)

where a₀ = 1, a1=1p$\begin{array}{} \displaystyle a_{1}=\dfrac{1}{p} \end{array}$, a2=p+12p2$\begin{array}{} \displaystyle {a_2} = \frac{{p + 1}}{{2{p^2}}} \displaystyle \end{array}$, P₀ = 1 and Pi=∑j=0i−1Pjai−ja0(i−j)(p+1)−ii$\begin{array}{} \displaystyle \displaystyle{P_{i}=\sum_{j=0}^{i-1}P_{j}\,\frac{a_{i-j}}{a_{0}}\,\frac{(i-j)(p+1)-i}{i}} \end{array}$, for i = 1,2,...,2p, and a_i = 0, for all i ≥ 3. In addition, we have P₁ = P₂ = 1 and P₃ < 1.

Note that we can always obtain X_nA = AX_n. For this, it is enough to choose X₀ = I_r×r, as Guo does in [11]. Now, we establish the semilocal convergence of Chebyshev’s method in the following theorem.