# Approximation of thepth Roots of a Matrix by Using Trapezoid Rule

Authors: Amir Sadeghi; Ahmad Izani Md. Ismail
Journal: International Journal of Mathematics and Mathematical Sciences
 (2012)

---

## Abstract

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


---

## Body

## 1. 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 eA 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 matrix-vector 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 A1/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 thepth 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.


## 2. Approximation of the Matrixpth Roots

The Cauchy integral theorem states that the value off(a) can be evaluated by an integral representation as follows:
(2.1)f(a)=12πi∮Γf(z)(z-a)dz,
where Γ is a contour in ℂ 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
(2.2)f(A)=12πi∮Γf(z)(zI-A)-1dz
and can be defined element by element as follows:
(2.3)f(A)=fij⟹fkj=12πi∮Γf(z)ekT(zI-A)-1ejdz,
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 zc, defined by z=zc+reiθ. This substitution gives us the following identity [1]:
(2.4)f(A)=12π∫02πf(z)(zI-A)-1reiθdθ.
Writing z-zc=reiθ and substituting into (2.4) gives
(2.5)f(A)=12π∫02π(z-zc)f(z)(zI-A)-1dθ.
A primary pth root of a square matrix A∈ℂn×n, with a p positive integer, is a solution of the matrix equation Xp-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 Sp={z∈ℂ∖{0}:|arg(z)|<π/p} [5]. In this paper, we demonstrate the method for f(A)=A1/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]:
(2.6)∫abf(x)dx≈h2(f(x0)+2∑j=1N-1f(xj)+f(xN)),
where h=(b-a)/N and xj=a+jh,j=1,…,N-1. Since f(a)=f(b), (2.6) can be reformulated as
(2.7)∫abf(x)dx≈b-aN∑j=0N-1f(xj).
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
(2.8)f(A)≈1N∑j=0N-1g(θj),
which is the general formula for the computation of a matrix function when Γ is a circle [1].The functionf(z)=z1/p is analytic (i.e., fx=-ify) everywhere in ℂ except at z=0. Consider a matrix A which has eigenvalues in the unit disk centered at zc. The contour is a disk of radius r, parameterized as z=zc+reiθ, and from (2.5) we have
(2.9)A1/p≈1N∑j=0N-1(zj-zc)zj1/p(zjI-A)-1.
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 withzc=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).Figure 1
Error in trapezoid rule for computingA1/p as matrix size increases.


## 3. Employing Matrix Decompositions

Matrix factorizations are utilized to compute thepth 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.


### 3.1. 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).
LetA∈ℂn×n; then there exists a unitary Q∈ℂn×n such that
(3.1)Q*AQ=T=D+N,
where D=diag(λ1,…,λn) and N is strictly upper triangular. Further, Q=diag(q1,…,qn) is a column partitioning of the unitary matrix Q where qi is referred to as Schur vectors, and from AQ=QT Schur vectors satisfy
(3.2)Aqk=λkqk+∑i=1k-1nikqi,k=1,…,n.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 factorizationA=QTQT.(ii)
MatrixT is upper triangular and so we then set Rj,j=Tj,j1/p.(iii)
Else operate column-by-column onT to produce the upper triangular pth root matrix R.This algorithm uses(28+(p-1)/3)n3 flops in total. The matrix pth root is given as B=QRQT. It can be verified that
(3.3)T=Rp⟹QTQT=QRpQT⟹A=(QRQT)…(QRQT)⟹A=Bp.
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
(3.4)zjI-A=zjI-QTQ*=Q(zjI-T)Q*.
Using this in (2.9) gives
(3.5)A1/p≈1NQ(∑j=0N-1(zj-zc)zj1/p(zjI-T)-1)Q*.


### 3.2. Eigenvalue Decomposition

This factorization is also called spectral decomposition and is presented as follows [4].Theorem 3.2 (Eigenvalue decomposition theorem).
LetA∈ℂn×n; there exists a nonsingular X∈ℂn×n which can diagonalize A(3.6)X-1AX=diag⁡(λ1,…,λn)
if and only if the geometric multiplicities of all eigenvalue λi are equal to their algebraic multiplicities.
Utilizing the property ofX, one has
(3.7)zjI-A=zjI-XDX-1=X(zjI-D)X-1.
Replacing this into (2.9) yields
(3.8)A1/p≈1NX(∑j=0N-1(zj-zc)zj1/p(zjI-D)-1)X-1.


### 3.3. Hessenberg Decomposition

This factorization is also called spectral decomposition and is presented as follows [4].Theorem 3.3 (Hessenberg decomposition theorem).
LetA∈ℝn×n; then there exists a orthogonal matrix Q∈ℝn×n such that
(3.9)QTAQ=H,
where H is a Hessenberg matrix which means that the elements below the subdiagonal are zero.
Applying the Hessenberg factorization and the orthogonality ofQ, one can write
(3.10)zjI-A=zjI-QHQT=Q(zjI-H)QT.
Substituting this into (2.9) will give
(3.11)A1/p≈1NQ(∑j=0N-1(zj-zc)zj1/p(zjI-H)-1)QT.


## 3.1. 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).
LetA∈ℂn×n; then there exists a unitary Q∈ℂn×n such that
(3.1)Q*AQ=T=D+N,
where D=diag(λ1,…,λn) and N is strictly upper triangular. Further, Q=diag(q1,…,qn) is a column partitioning of the unitary matrix Q where qi is referred to as Schur vectors, and from AQ=QT Schur vectors satisfy
(3.2)Aqk=λkqk+∑i=1k-1nikqi,k=1,…,n.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 factorizationA=QTQT.(ii)
MatrixT is upper triangular and so we then set Rj,j=Tj,j1/p.(iii)
Else operate column-by-column onT to produce the upper triangular pth root matrix R.This algorithm uses(28+(p-1)/3)n3 flops in total. The matrix pth root is given as B=QRQT. It can be verified that
(3.3)T=Rp⟹QTQT=QRpQT⟹A=(QRQT)…(QRQT)⟹A=Bp.
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
(3.4)zjI-A=zjI-QTQ*=Q(zjI-T)Q*.
Using this in (2.9) gives
(3.5)A1/p≈1NQ(∑j=0N-1(zj-zc)zj1/p(zjI-T)-1)Q*.


## 3.2. Eigenvalue Decomposition

This factorization is also called spectral decomposition and is presented as follows [4].Theorem 3.2 (Eigenvalue decomposition theorem).
LetA∈ℂn×n; there exists a nonsingular X∈ℂn×n which can diagonalize A(3.6)X-1AX=diag⁡(λ1,…,λn)
if and only if the geometric multiplicities of all eigenvalue λi are equal to their algebraic multiplicities.
Utilizing the property ofX, one has
(3.7)zjI-A=zjI-XDX-1=X(zjI-D)X-1.
Replacing this into (2.9) yields
(3.8)A1/p≈1NX(∑j=0N-1(zj-zc)zj1/p(zjI-D)-1)X-1.


## 3.3. Hessenberg Decomposition

This factorization is also called spectral decomposition and is presented as follows [4].Theorem 3.3 (Hessenberg decomposition theorem).
LetA∈ℝn×n; then there exists a orthogonal matrix Q∈ℝn×n such that
(3.9)QTAQ=H,
where H is a Hessenberg matrix which means that the elements below the subdiagonal are zero.
Applying the Hessenberg factorization and the orthogonality ofQ, one can write
(3.10)zjI-A=zjI-QHQT=Q(zjI-H)QT.
Substituting this into (2.9) will give
(3.11)A1/p≈1NQ(∑j=0N-1(zj-zc)zj1/p(zjI-H)-1)QT.


## 4. 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 ifX̃ is approximated value of A1/p by using different methods, then the absolute error and residual errors can be considered as follows:
(4.1)e(X̃)=‖X̃p-A‖F,Res(X̃)=‖X̃p-A‖F‖A‖F,
respectively, where ∥·∥F is Frobenius norm.Test 1
For the first experiment, consider20×20 Random matrix in the form
(4.2)A=randn(N)/sqrt(N)+3*eye(N),
which has positive eigenvalues. We estimate the pth root of A for different values of p using trapezoid rule. Furthermore, the Schur, eigenvalue, and Hessenberg factorizations are utilized for speeding up the computations. Moreover, the residual error Res(X̃) and used time 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.Table 1
Comparison residual error for computing thepth root of Random matrix.
pSchur decompositionEigenvalue decompositionHessenberg decomposition24.660633743656308e-0155.030320198311745e-0151.554377109452457e-015162.311964104013483e-0143.524734357427803e-0147.749850409855744e-015528.706212825269891e-0141.121352227197161e-0132.904341268646203e-0141282.358332191022634e-0132.524139574310145e-0137.105511057327197e-01420123.448827113054843e-0123.472268532385117e-0121.170919968312386e-012Table 2
Comparison time in seconds for computing thepth root of Random matrix.
pSchur decompositionEigenvalue decompositionHessenberg decomposition20.0548550.0284010.017409160.1617440.2265240.141455520.4047490.7613850.6638781280.9689701.9141471.80349020120.0133310.0250150.018379Test 2
In this example consider10×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 number 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.Table 3
Comparison absolute error for different numbers of points for different matrices.
No. of pointsSchur decompositionEigenvalue decompositionHessenberg decomposition24.2354768914775933.6176967113303099.01827121836307047.6428225550780702.26854880846257385.56976811035167486.50973162179347216.0595306707345210.951463218971532160.0359106106439740.3429295484835090.001435059307873320.0000027708870730.3429295484835090.000000210760699640.0000000000003750.0034292954848350.0000000000003491280.0000000000000020.0034292954848350.0000000000000022560.0000000000000020.0034292954848350.0000000000000025120.0000000000000040.0034292954848350.00000000000000310240.0000000000000010.0003429295484840.000000000000000Test 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 ofA and absolute errors and also time in seconds are demonstrated in these figures. Figures 2–7 show the comparison of the residual errors and timings for these methods.
According to Figure2 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 220 points in trapezoid rule can achieve error of 10-6 in the last experiment.Residual error for computingA1/p for Random matrix.
(a)(b)(c)(d)Time in seconds for computingA1/p for Random matrix.
(a)(b)(c)(d)Residual error for computingA1/p for Hilbert matrix.
(a)(b)(c)(d)Time in seconds for computingA1/p for Hilbert matrix.
(a)(b)(c)(d)Residual error for computingA1/p for Lehmer matrix.
(a)(b)(c)(d)Time in seconds for computingA1/p for Lehmer matrix.
(a)(b)(c)(d)


## 5. Conclusion

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



---

*Source: 634698-2012-01-24.xml*