On High-Order Iterative Schemes for the Matrix p th Root Avoiding the Use of Inverses

: This paper is devoted to the approximation of matrix p th roots. We present and analyze a family of algorithms free of inverses. The method is a combination of two families of iterative methods. The first one gives an approximation of the matrix inverse. The second family computes, using the first method, an approximation of the matrix p th root. We analyze the computational cost and the convergence of this family of methods. Finally, we introduce several numerical examples in order to check the performance of this combination of schemes. We conclude that the method without inverse emerges as a good alternative since a similar numerical behavior with smaller computational cost is obtained.


Introduction
The computation of operators on a matrix appears in many applications. The use of iterative methods has emerged as a useful technique to approximate it. This paper is devoted to the approximation of the matrix pth root. We recall that the principal pth root of a matrix A ∈ C r×r , where C is a set of complex numbers, without positive eigenvalues in R, as the unique solution X of with an argument of the eigenvalues in modulus lower than π p .
For this kind of problem, the well-known and famous Newton method takes the following form X 0 given, X n+1 = (p − 1)X n + AX 1−p n p , n ≥ 0.
This formula is a direct extension of the method applied to the scalar equation x p − a = 0. In this form, the method achieves quadratic convergence, but unfortunately, it is unstable [1]. However, as can be in seen in [2], a stable version can be found.
Similarly, we can derive the third-order Chebyshev method as where Θ = B ∈∈ C r×r s.t. B has no nonpositive real eigenvalues. In our paper [7], we proposed stable versions of this algorithm, we presented some numerical advantages of it with respect to Newton and other third-order methods such as Halley's method [8,9], and finally, we developed a general family of any order that includes both the Newton and Chebyshev methods.
In order to develop a method avoiding the use of inverses of the different iterates, we can consider the approximation of this less natural equation 1 x p − 1 a = 0. In this case, the method presented in [7] has the form      X 0 given, where This method has order m, and in particular, for m = 2 and m = 3, we recover the Newton and Chebyshev methods.
In the above Formula (2), we only need to compute the inverse of A. In the present paper, we propose to approximate the inverse of A by another iterative method. We can use our family introduced in [10] that has the form Our goal is to see that the new approach yields a similar numerical behavior but with the advantage that it is free of inverse operators, and in particular, has a smaller computational cost, which means that it requires fewer operations.
The computation of the pth root of a matrix appears, for example, in fractional differential equations, discrete representations of norms corresponding to finite element discretizations of fractional Sobolev spaces, and the computation of geodesic-midpoints in neural networks (see [11] and the references therein).

A General Method for Approximating the Matrix pth Root Free of Inverse Operators
As we mentioned in the introduction, for approximating A 1 p , we propose the use of a combination of two families. One for the approximation of A −1 and other, that use this approximation, for the approximation of the matrix pth root.
First step: for approximating

Second step: for approximating
where Y L denotes the final iteration computed by the above method (4) in the first step.

Convergence
To establish the convergence of the general method (5), we observe that Now, suppose that, from method (4), we consider n 0 ∈ N such that Y n 0 = Y L . Then, for the iterative process (4), using Theorem 3.2 in [10], we obtain On the other hand, for the iterative process (2), from Theorem 4.1 in [7], we have Finally, fixed Y L we take X 0 such that I − Y L X p 0 < 1. Then, from Theorem 4.1 in [7], the sequence {X n }, given by the general method (5), converges to (Y L ) 1 p ; therefore, {X n } is bounded. So, there exists M > 0 such that X n ≤ M for all n ∈ N.
then for all tolerance Tol > 0 given, there exists n * ∈ N such that I − A −1 X p n < Tol for n n * .
Proof. Given Tol > 0, from (8), there exists n 1 ∈ N such that I − Y L X p n < Tol 2 for all n n 1 .
On the other hand, as X n converges to (Y L ) 1 p , from Theorem 4.1 in [7], there exists Moreover, from (7), given Tol Tol 2 for all n n 3 . Therefore, from (6), there exists n * = max{n 1 , n 2 , n 3 } ∈ N such that I − A −1 X p n < Tol for all n n * .

Computational Cost and Efficiency
If A is a matrix q × q, and taking into account that the computational cost of the product of two matrices is q 3 and that adding or subtracting two matrices has a computational cost of q 2 , then the computational cost of • the computation of matrix Y L , given by means of using (4), is CC I (m, q) = mq 2 + (m + 1)q 3 and its computational efficiency is CE I (m, q) = (m + 1) the computation of method (2) has a computational cost of CC R (m, q) = mq 2 + (p + m)q 3 and its computational efficiency is CE R (m, q) = (m + 1) the computational cost of method (5) is obtained directly by CC(m, q) = 2mq 2 + (p + 2m + 1)q 3 and its computational efficiency is CE(m, q) = (m + 1) 1 2mq 2 +(p+2m+1)q 3 .
In Figures 1 and 2 the computational efficiency for fixed p = 2 and p = 3 respectively and different values of q and m are shown.  In most cases, when working with a family of methods, the best one is related to the problem considered (see the next figures) and the main way used to solve the problem. For example, the most efficient method (in terms of computation) for solving sparse problems, which appear when discretizing differential equations, is Chebyshev's method [7,10].
In Figures 3 and 4 the computational efficiency for fixed q = 15 and q = 20 respectively and different values of p and m are shown.

Applications Related to Differential Equations
In this section, we present a comparison of the original method (2) using A −1 and the new proposal (5), avoiding the use of any inverse operator. We refer to our paper [7] in order to see the advantages of the original method in comparison with other methods that appear in the literature.
We consider the approximation of the pth root of two matrices related to the discretization of differential equations. This type of matrix operations appears in the approximation of space-fractional diffusion problems [12].
We start with the following matrix where B = −1 − h 2 w, C = 2 + h 2 v and D = −1 + h 2 w. The matrix (9) can be seen as the result of applying a discretization process, using finite differences to the boundary value problem defined as Since the matrix A is sparse, the most efficient methods in both families are Chebyshevlike methods.
In Table 1, we compare our original method using A −1 (2) and the new combination using the approximation Y L of the inverse (5). We observe a similar numerical behavior of both methods. Thus, the approximation of the inverse seems a good alternative since has similar errors with smaller computational cost. Finally, we compute the pth matrix roots for the matrix This matrix is related to the discretization, using finite differences, of the laplacian operator that appear in many mathematical models. In Table 2, we observe again a similar numerical behavior of both schemes.

Conclusions
The function evaluation of a matrix appears in a large and growing number of applications. We have presented and studied a general family of iterative methods without using inverse operators, for the approximation of the matrix pth root. The family incorporates the approximation of A −1 , by an iterative method, into another iterative method for approximating A 1/p . The family includes methods of every order of convergence.
As it appears in [7,10], the most efficient method, in computational terms, for solving sparse problems, which appears when discretizing differential equations, is Chebyshev's method.
Some numerical examples related to the discretization of differential equations have been presented. We have concluded that the new approach (5) has a similar numerical behavior to that of (2) but with the advantage that it is free of inverse operators, and in particular, has a lower computational cost. Finally, we care to mention the existence of other strategies for computing fractional powers of a matrix, which do not use matrix iterations [13,14] or even take into account other techniques such as those appearing in [15].