Approximation of the p th Roots of a Matrix by Using Trapezoid Rule

The computation of the roots of positive definite matrices arises in nuclear magnetic resonance, control theory, lattice quantum chromo-dynamics QCD , and several other areas of applications. The Cauchy integral theorem which arises in complex analysis can be used for computing f(A), in particular the roots of A, where A is a square matrix. The Cauchy integral can be approximated by using the trapezoid rule. In this paper, we aim to give a brief overview of the computation of roots of positive definite matrices by employing integral representation. Some numerical experiments are given to illustrate the theoretical results.


Introduction
It is well known that contour integrals which form a component of the Cauchy integral theorem have an important role in complex analysis.The trapezoid rule is popular for the approximation of integrals due to its exponential accuracy if particular conditions are satisfied.It has been established in 1 that the trapezoid rule can be used to compute contour integrals to give a powerful algorithm for the computation of matrix functions.
Kellems 1 has studied the case of computing a matrix square root and a matrix exponential function by utilizing the trapezoid rule.In particular, he focused on the matrix exponential e A and its use in the heat equation.Only a few trapezoid rule points were required for very high accuracy.Davies and Higham 2 have investigated the computation of a matrixvector product f A b without explicitly computing f A .Their proposed methods were specific to the logarithm and fractional matrix powers which were based on quadrature and solution of an ordinary differential equation initial value problem, respectively.Hale et al. in 3 have presented new methods for the numerical computation of f A and f A b, where f A is a function such as A 1/2 or log A with singularities in −∞, 0 and A is a matrix with eigenvalues on or near 0, ∞ .The methods in 3 were based on combining contour integrals evaluated using the periodic trapezoid rule with conformal maps involving elliptic functions.
In this paper, we investigate computation of the pth roots of positive definite matrices by utilizing integral representation.Our approach is based on the work of Kellems 1 who compute the square root of positive definite matrix using trapezoid rule.The integral identity will be computed by employing trapezoid rule.We also study some matrix factorization procedure and their application in computing matrix roots using trapezoid rule.The outline of this paper is as follows.In Section 2 we introduce some basic definitions and also integral representation of function of matrices and trapezoid rule.In Section 3 we will obtain some formulas to compute the integral representation of the matrix pth root.Numerical experiments will be discussed in Section 4, and the conclusions will be presented in Section 5.

Approximation of the Matrix pth Roots
The Cauchy integral theorem states that the value of f a can be evaluated by an integral representation as follows: where Γ is a contour in C such that Γ enclose a and f z is an analytic and inside Γ 1 .The generalization of this formula in the matrix case can be presented as and can be defined element by element as follows: where the entries of zI −A −1 are analytic on Γ and also f z is analytic function in the neighborhood of the spectrum of A 4 .Cauchy's integral formula can be simplified by considering Γ to be a circle of radius r centered at some point z c , defined by z z c re iθ .This substitution gives us the following identity 1 : Writing z − z c re iθ and substituting into 2.4 gives A primary pth root of a square matrix A ∈ C n×n , with a p positive integer, is a solution of the matrix equation X p − A 0 that can be written as a polynomial of A. If A has distinct eigenvalues and none of which is zero, then A has exactly p primary pth roots.This result by 2.2 was obtained where f is any of the p analytic functions defined on the spectrum of A, such that f z p z and Γ is a closed contour which encloses the spectrum of A. If A has no nonpositive real eigenvalues, then there exists only one primary pth root whose eigenvalues lie in the sector S p {z ∈ C \ {0} : | arg z | < π/p} 5 .In this paper, we demonstrate the method for f A A 1/p .In order to calculate the integral 2.2 accurately, we first split the interval of integration a, b into N smaller uniform subintervals and then apply the trapezoidal rule on each of them.The composite trapezoidal rule is as follows 6 : where h b − a /N and x j a jh, j 1, . . ., N − 1.Since f a f b , 2.6 can be reformulated as For 2.5 , let the integrand be the function g θ .If we take N equally spaced points on Γ and consider that f 0 f 2π , then 2.7 can be written as which is the general formula for the computation of a matrix function when Γ is a circle 1 .
The function f z z 1/p is analytic i.e., f x −if y everywhere in C except at z 0. Consider a matrix A which has eigenvalues in the unit disk centered at z c .The contour is a disk of radius r, parameterized as z z c re iθ , and from 2.5 we have An important property of the trapezoidal rule approximation is that it has better accuracy than the standard matrix pth root algorithms 1 .Now we suppose the Random matrix and use trapezoid rule with z c 3 and r 2. Figure 1 shows us the convergence of this method for matrices of dimension 4, 8, 16, 32, and 64.In each case exponential accuracy was found.Since the eigenvalues are well clustered, a few points need to be used.For matrices with larger spectral radius or more scattered eigenvalues, the convergence will be slower.This is consistent with the finding in 1 for the case of square root p 2 .

Employing Matrix Decompositions
Matrix factorizations are utilized to compute the pth root of a matrix using trapezoid rule.Given a square matrix A, we are interested as to the simplest form of matrix B in C or R under unitary similarity transform A QBQ * or similarity transform A XBX −1 .Matrix B presents some information on A because many features and structure of matrices are invariant under similarity transform.In this part, three factorizations: Schur, Eigenvalue, and Hessenberg are investigated in relation to the use of the trapezoidal rule.

Schur Decomposition
One of the most applicable factorization of matrices is Schur decomposition which is presented in the following theorem 4 .

Theorem 3.1 Schur decomposition theorem . Let
where D diag λ 1 , . . ., λ n and N is strictly upper triangular.Further, Q diag q 1 , . . ., q n is a column partitioning of the unitary matrix Q where q i is referred to as Schur vectors, and from AQ QT Schur vectors satisfy One of the most famous algorithms to compute matrix roots is Smith's algorithm proposed in 7 .Generally, this algorithm is presented as follows.
i Compute the Schur factorization A QT Q T .
ii Matrix T is upper triangular and so we then set R j,j T 1/p j,j .
iii Else operate column-by-column on T to produce the upper triangular pth root matrix R.
This algorithm uses 28 p − 1 /3 n 3 flops in total.The matrix pth root is given as B QRQ T .It can be verified that

3.3
This can be used to speed up the trapezoid rule method: implement a preliminary factorization of A, operate on the factored matrix, and then combine the factors at the end of the computation 1 .Using the Schur factorization and the unitary of Q, we then have Using this in 2.9 gives

Eigenvalue Decomposition
This factorization is also called spectral decomposition and is presented as follows 4 .
Theorem 3.2 Eigenvalue decomposition theorem .Let A ∈ C n×n ; there exists a nonsingular if and only if the geometric multiplicities of all eigenvalue λ i are equal to their algebraic multiplicities.
Utilizing the property of X, one has

International Journal of Mathematics and Mathematical Sciences
Replacing this into 2.9 yields 3.8

Hessenberg Decomposition
This factorization is also called spectral decomposition and is presented as follows 4 .
Theorem 3.3 Hessenberg decomposition theorem .Let A ∈ R n×n ; then there exists a orthogonal matrix Q ∈ R n×n such that where H is a Hessenberg matrix which means that the elements below the subdiagonal are zero.
Applying the Hessenberg factorization and the orthogonality of Q, one can write

3.10
Substituting this into 2.9 will give 3.11

Numerical Experiments
In this section we present some numerical experiment to illustrate the theory which is developed.All the computations have been carried out using MATLAB 7.10 Ra .We assume positive definite matrices with positive nonzero eigenvalues.These matrices, which are given in MATLAB gallery, are used to compute roots of matrices.Recall that if X is approximated value of A 1/p by using different methods, then the absolute error and residual errors can be considered as follows: in proposed method for the computation of matrix pth root, for p 2, 16, 52, 128, and 2012, are compared.The results are observed in Tables 1 and 2. It should be mentioned that the number of the point in trapezoid rule is considered as N 128.As can be shown in the result, the Schur factorization has almost exactly the same error as MATLAB's algorithm.The Hessenberg factorization gives the best accuracy among the three methods.In fact, the Hessenberg factorization is quicker than some algorithms in MATLAB.The trapezoid rule for the computation of the matrix pth root may be more effective than several other algorithms but this is dependent on the spectrum of A. This is consistent with the finding in 1 for the case of square root.

Test 2
In this example consider 10 × 10 Random, Hilbert, and Lehmer matrices.We first fix the value of p and then use trapezoid rule to compute matrix pth root.The relation between the num-   ber of points in trapezoid rule and obtained absolute errors e X is investigated.In the implementation, p 100 is supposed and the number of points is increased as N 2 , for 1, . . ., 10.Comparison between errors for various matrices is illustrated in Table 3.It can be seen that errors by increasing the number of points in trapezoid rule will be reduced.

Test 3
In the last test Random matrix, Lehmer matrix, and Hilbert matrix are supposed.We have computed the 5th, 17th, 64th, and 128th root of these matrices using trapezoidal rule and Smith's algorithm.Furthermore, in this example the absolute errors are estimated.As shown in the figures, the accuracy is measured in either the Frobenius norm or the 2-norm for different matrices.The relations between dimension of A and absolute errors and also time in     2 which has four parts, the absolute error in all cases p 5, 17, 64, and 128 in trapezoid rule are smaller than Smith's method.In addition, in Figure 3 it is illustrated that except the case of p 5, the time of computation of the pth root using Smith's method is longer than trapezoidal rule.Also, Figures 4 and 5 are given difference between error and also time in trapezoid rule and Smith's method for Hilbert matrix.The residual error for trapezoid rule is large while for Smith's algorithm is small.For Hilbert matrix except case p 5 , in all cases Smith's method is more time consuming than trapezoid rule.Finally, Figures 6 and 7 show the computation of the pth root of Lehmer matrix which in the most cases reveal that trapezoid rule is more expensive in time and also error than Smith's method by using N 128 points.It must be mentioned that by increasing the number of points, considerable accurate solution can be obtained.For example, using 2 20 points in trapezoid rule can achieve error of 10 −6 in the last experiment.

Conclusion
In this paper, we have studied the use of trapezoidal rule in conjunction with the Cauchy integral theorem to compute the pth roots of matrices.It was found that the technique is feasible and accurate.

N
= number of points in trapezoid rule

Figure 1 :
Figure 1: Error in trapezoid rule for computing A 1/p as matrix size increases.

Figure 4 :
Figure 4: Residual error for computing A 1/p for Hilbert matrix.

Figure 5 :
Figure 5: Time in seconds for computing A 1/p for Hilbert matrix.

Figure 6 :
Figure 6: Residual error for computing A 1/p for Lehmer matrix.

Figure 7 :
Figure 7: Time in seconds for computing A 1/p for Lehmer matrix.

Table 1 :
Comparison residual error for computing the pth root of Random matrix.

Table 2 :
Comparison time in seconds for computing the pth root of Random matrix.

Table 3 :
Comparison absolute error for different numbers of points for different matrices.aredemonstrated in these figures.Figures2-7show the comparison of the residual errors and timings for these methods.According to Figure seconds