Approximation of the pth Roots of a Matrix by Using Trapezoid Rule

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

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 𝑒^𝐴 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 𝑓(𝐴)𝑏 without explicitly computing 𝑓(𝐴). 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 𝑓(𝐴) and 𝑓(𝐴)𝑏, where 𝑓(𝐴) is a function such as 𝐴^1/2 or log(𝐴) with singularities in (−∞,0] and 𝐴 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 𝑝th 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 𝑝th root. Numerical experiments will be discussed in Section 4, and the conclusions will be presented in Section 5.

2. Approximation of the Matrix pth Roots

The Cauchy integral theorem states that the value of 𝑓(𝑎) can be evaluated by an integral representation as follows: 𝑓(𝑎)=12𝜋𝑖_Γ𝑓(𝑧)(𝑧−𝑎)𝑑𝑧, (2.1) where Γ is a contour in ℂ such that Γ enclose 𝑎 and 𝑓(𝑧) is an analytic and inside Γ [1]. The generalization of this formula in the matrix case can be presented as 𝑓(𝐴)=12𝜋𝑖_Γ𝑓(𝑧)(𝑧𝐼−𝐴)⁻¹𝑑𝑧 (2.2) and can be defined element by element as follows: 𝑓(𝐴)=𝑓_𝑖𝑗⟹𝑓_𝑘𝑗=12𝜋𝑖_Γ𝑓(𝑧)𝑒^𝑇_𝑘(𝑧𝐼−𝐴)⁻¹𝑒_𝑗𝑑𝑧, (2.3) where the entries of (z𝐼−𝐴)⁻¹ are analytic on Γ and also 𝑓(𝑧) is analytic function in the neighborhood of the spectrum of 𝐴 [4]. Cauchy's integral formula can be simplified by considering Γ to be a circle of radius 𝑟 centered at some point 𝑧_𝑐, defined by 𝑧=𝑧_𝑐+𝑟𝑒^𝑖𝜃. This substitution gives us the following identity [1]: 𝑓(𝐴)=12𝜋₀^2𝜋𝑓(𝑧)(𝑧𝐼−𝐴)⁻¹𝑟𝑒^𝑖𝜃𝑑𝜃. (2.4) Writing 𝑧−𝑧_𝑐=𝑟𝑒^𝑖𝜃 and substituting into (2.4) gives 𝑓(𝐴)=12𝜋₀^2𝜋𝑧−𝑧_𝑐𝑓(𝑧)(𝑧𝐼−𝐴)⁻¹𝑑𝜃. (2.5) A primary 𝑝th root of a square matrix 𝐴∈ℂ^𝑛×𝑛, with a 𝑝 positive integer, is a solution of the matrix equation 𝑋^𝑝−𝐴=0 that can be written as a polynomial of 𝐴. If 𝐴 has ℓ distinct eigenvalues and none of which is zero, then 𝐴 has exactly 𝑝^ℓ primary 𝑝th roots. This result by (2.2) was obtained where 𝑓 is any of the 𝑝^ℓ analytic functions defined on the spectrum of 𝐴, such that 𝑓(𝑧)^𝑝=𝑧 and Γ is a closed contour which encloses the spectrum of 𝐴. If 𝐴 has no nonpositive real eigenvalues, then there exists only one primary 𝑝th root whose eigenvalues lie in the sector 𝑆_𝑝={𝑧∈ℂ⧵{0}∶|arg(𝑧)|<𝜋/𝑝} [5]. In this paper, we demonstrate the method for 𝑓(𝐴)=𝐴^1/𝑝.

In order to calculate the integral (2.2) accurately, we first split the interval of integration [𝑎,𝑏] into 𝑁 smaller uniform subintervals and then apply the trapezoidal rule on each of them. The composite trapezoidal rule is as follows [6]: ^𝑏_𝑎𝑓(𝑥)ℎ𝑑𝑥≈2𝑓𝑥₀+2^𝑁−1_𝑗=1𝑓𝑥_𝑗+𝑓𝑥_𝑁, (2.6) where ℎ=(𝑏−𝑎)/𝑁 and 𝑥_𝑗=𝑎+𝑗ℎ,𝑗=1,…,𝑁−1. Since 𝑓(𝑎)=𝑓(𝑏), (2.6) can be reformulated as ^𝑏_𝑎𝑓(𝑥)𝑑𝑥≈𝑏−𝑎𝑁^𝑁−1_𝑗=0𝑓𝑥_𝑗. (2.7) For (2.5), let the integrand be the function 𝑔(𝜃). If we take 𝑁 equally spaced points on Γ and consider that 𝑓(0)=𝑓(2𝜋), then (2.7) can be written as 𝑓(𝐴)≈1𝑁^𝑁−1_𝑗=0𝑔𝜃_𝑗, (2.8) which is the general formula for the computation of a matrix function when Γ is a circle [1].

The function 𝑓(𝑧)=𝑧^1/𝑝 is analytic (i.e., 𝑓_𝑥=−𝑖𝑓_𝑦) everywhere in ℂ except at 𝑧=0. Consider a matrix 𝐴 which has eigenvalues in the unit disk centered at 𝑧_𝑐. The contour is a disk of radius 𝑟, parameterized as 𝑧=𝑧_𝑐+𝑟𝑒^𝑖𝜃, and from (2.5) we have 𝐴^1/𝑝≈1𝑁^𝑁−1_𝑗=0𝑧_𝑗−𝑧_𝑐𝑧_𝑗^1/𝑝𝑧_𝑗𝐼−𝐴⁻¹. (2.9) An important property of the trapezoidal rule approximation is that it has better accuracy than the standard matrix 𝑝th root algorithms [1].

Now we suppose the Random matrix and use trapezoid rule with 𝑧_𝑐=3 and 𝑟=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 (𝑝=2).

Figure 1 

Error in trapezoid rule for computing 𝐴1/𝑝 as matrix size increases.

3. Employing Matrix Decompositions

Matrix factorizations are utilized to compute the 𝑝th root of a matrix using trapezoid rule. Given a square matrix 𝐴, we are interested as to the simplest form of matrix 𝐵 in 𝐶 or 𝑅 under unitary similarity transform 𝐴=𝑄𝐵𝑄^∗ or similarity transform 𝐴=𝑋𝐵𝑋⁻¹. Matrix 𝐵 presents some information on 𝐴 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). Let 𝐴∈ℂ^𝑛×𝑛; then there exists a unitary 𝑄∈ℂ^𝑛×𝑛 such that 𝑄^∗𝐴𝑄=𝑇=𝐷+𝑁, (3.1) where 𝐷=diag(𝜆₁,…,𝜆_𝑛) and 𝑁 is strictly upper triangular. Further, 𝑄=diag(𝑞₁,…,𝑞_𝑛) is a column partitioning of the unitary matrix 𝑄 where 𝑞_𝑖 is referred to as Schur vectors, and from 𝐴𝑄=𝑄𝑇 Schur vectors satisfy 𝐴𝑞_𝑘=𝜆_𝑘𝑞_𝑘+^𝑘−1_𝑖=1𝑛_𝑖𝑘𝑞_𝑖,𝑘=1,…,𝑛. (3.2)

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 𝐴=𝑄𝑇𝑄^𝑇.(ii)Matrix 𝑇 is upper triangular and so we then set 𝑅_𝑗,𝑗=𝑇^1/𝑝_𝑗,𝑗.(iii)Else operate column-by-column on 𝑇 to produce the upper triangular 𝑝th root matrix 𝑅.

This algorithm uses (28+(𝑝−1)/3)𝑛³ flops in total. The matrix 𝑝th root is given as 𝐵=𝑄𝑅𝑄^𝑇. It can be verified that 𝑇=𝑅^𝑝⟹𝑄𝑇𝑄^𝑇=𝑄𝑅^𝑝𝑄^𝑇⟹𝐴=𝑄𝑅𝑄^𝑇…𝑄𝑅𝑄^𝑇⟹𝐴=𝐵^𝑝. (3.3) This can be used to speed up the trapezoid rule method: implement a preliminary factorization of 𝐴, 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 𝑄, we then have 𝑧_𝑗𝐼−𝐴=𝑧_𝑗𝐼−𝑄𝑇𝑄^∗=𝑄𝑧_𝑗𝐼−𝑇𝑄^∗. (3.4) Using this in (2.9) gives 𝐴^1/𝑝≈1𝑁𝑄^𝑁−1_𝑗=0𝑧_𝑗−𝑧_𝑐𝑧_𝑗^1/𝑝𝑧_𝑗𝐼−𝑇⁻¹𝑄^∗. (3.5)

3.2. Eigenvalue Decomposition

This factorization is also called spectral decomposition and is presented as follows [4].

Theorem 3.2 (Eigenvalue decomposition theorem). Let 𝐴∈ℂ^𝑛×𝑛; there exists a nonsingular 𝑋∈ℂ^𝑛×𝑛 which can diagonalize 𝐴𝑋⁻¹𝐴𝑋=diag𝜆₁,…,𝜆_𝑛 (3.6) if and only if the geometric multiplicities of all eigenvalue 𝜆_𝑖 are equal to their algebraic multiplicities.
Utilizing the property of 𝑋, one has 𝑧_𝑗𝐼−𝐴=𝑧_𝑗𝐼−𝑋𝐷𝑋⁻¹=𝑋𝑧_𝑗𝐼−𝐷𝑋⁻¹. (3.7) Replacing this into (2.9) yields 𝐴^1/𝑝≈1𝑁𝑋^𝑁−1_𝑗=0𝑧_𝑗−𝑧_𝑐𝑧_𝑗^1/𝑝𝑧_𝑗𝐼−𝐷⁻¹𝑋⁻¹. (3.8)

3.3. Hessenberg Decomposition

This factorization is also called spectral decomposition and is presented as follows [4].

Theorem 3.3 (Hessenberg decomposition theorem). Let 𝐴∈ℝ^𝑛×𝑛; then there exists a orthogonal matrix 𝑄∈ℝ^𝑛×𝑛 such that 𝑄^𝑇𝐴𝑄=𝐻, (3.9) where 𝐻 is a Hessenberg matrix which means that the elements below the subdiagonal are zero.
Applying the Hessenberg factorization and the orthogonality of 𝑄, one can write 𝑧_𝑗𝐼−𝐴=𝑧_𝑗𝐼−𝑄𝐻𝑄^𝑇=𝑄𝑧_𝑗𝐼−𝐻𝑄^𝑇. (3.10) Substituting this into (2.9) will give 𝐴^1/𝑝≈1𝑁𝑄^𝑁−1_𝑗=0𝑧_𝑗−𝑧_𝑐𝑧^1/𝑝_j𝑧_𝑗𝐼−𝐻⁻¹𝑄^𝑇. (3.11)

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 if 𝑋 is approximated value of 𝐴^1/𝑝 by using different methods, then the absolute error and residual errors can be considered as follows: 𝑒𝑋=𝑋^𝑝−𝐴_𝐹,𝑋=𝑋Res^𝑝−𝐴_𝐹‖𝐴‖_𝐹, (4.1) respectively, where ‖⋅‖_𝐹 is Frobenius norm.

Test 1
For the first experiment, consider 20×20 Random matrix in the form 𝐴=randn(𝑁)/sqrt(𝑁)+3∗eye(𝑁), (4.2) which has positive eigenvalues. We estimate the 𝑝th root of 𝐴 for different values of 𝑝 using trapezoid rule. Furthermore, the Schur, eigenvalue, and Hessenberg factorizations are utilized for speeding up the computations. Moreover, the residual error Res(𝑋) and used time in proposed method for the computation of matrix 𝑝th root, for 𝑝=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 𝑁=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 𝑝th root may be more effective than several other algorithms but this is dependent on the spectrum of 𝐴. This is consistent with the finding in [1] for the case of square root.

Table 1 

Comparison residual error for computing the 𝑝th root of Random matrix.

Table 2 

Comparison time in seconds for computing the 𝑝th root of Random matrix.

Test 2
In this example consider 10×10 Random, Hilbert, and Lehmer matrices. We first fix the value of 𝑝 and then use trapezoid rule to compute matrix 𝑝th root. The relation between the number of points in trapezoid rule and obtained absolute errors 𝑒(𝑋) is investigated. In the implementation, 𝑝=100 is supposed and the number of points is increased as 𝑁=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.

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 𝐴 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 Figure 2 which has four parts, the absolute error in all cases (𝑝=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 𝑝=5, the time of computation of the 𝑝th 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 (𝑝=5), in all cases Smith's method is more time consuming than trapezoid rule. Finally, Figures 6 and 7 show the computation of the 𝑝th 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 𝑁=128 points. It must be mentioned that by increasing the number of points, considerable accurate solution can be obtained. For example, using 2²⁰ points in trapezoid rule can achieve error of 10⁻⁶ in the last experiment.

(a)

(b)

(c)

(d)

Figure 2 

Residual error for computing 𝐴1/𝑝 for Random matrix.

(a)

(b)

(c)

(d)

Figure 3 

Time in seconds for computing 𝐴1/𝑝 for Random matrix.

(a)

(b)

(c)

(d)

Figure 4 

Residual error for computing 𝐴1/𝑝 for Hilbert matrix.

(a)

(b)

(c)

(d)

Figure 5 

Time in seconds for computing 𝐴1/𝑝 for Hilbert matrix.

(a)

(b)

(c)

(d)

Figure 6 

Residual error for computing 𝐴1/𝑝 for Lehmer matrix.

(a)

(b)

(c)

(d)

Figure 7 

Time in seconds for computing 𝐴1/𝑝 for Lehmer matrix.

5. Conclusion

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

Acknowledgments

The authors would like to thank Dr. Anthony Kellems and the anonymous referees for their helpful comments which improved the presentation.