Some Matrix Iterations for Computing Matrix Sign Function

Abstract

Some iterative methods are introduced and demonstrated for finding the matrix sign function. It is analytically shown that the new schemes are asymptotically stable. Convergence analysis along with the error bounds of the main proposed method is established. Different numerical experiments are employed to compare the behavior of the new schemes with the existing matrix iterations of the same type.

1. Introduction

Recently, the theory of matrix functions becomes an active topic of research in the field of advanced numerical linear algebra (see, e.g., [1–4]). In fact, the most common matrix function is the matrix inverse or the Moore-Penrose generalized inverse, routinely used in the scientific problems [5]. General matrix functions as well as the specific cases have been extensively discussed and developed in [6].

This paper is concerned with a special case known as matrix sign function, which is of clear importance in the theory of matrix functions [7]. Let us, as Higham considered in the fifth Chapter of [6], assume throughout this paper that the matrix 𝐴∈ℂ^𝑛×𝑛 has no eigenvalues on the imaginary axis. This assumption implies that the matrix sign function, 𝑆=sign(𝐴), (1) can be uniquely defined, whereas 𝐴 is a nonsingular square matrix. In order to define 𝑆, we remember the matrix sector function, for any positive integer 𝑝, can be defined by sect_𝑝(𝐴)=𝐴𝐴^{𝑝−1/𝑝}. (2)

Choosing 𝑝=2 in the matrix sector function will yield in the matrix sign function as 𝑆=𝐴(𝐴²)^−1/2. This also clearly puts on show the importance and the relevance of this matrix function to the other important matrix functions such as matrix square root.

Bini et al. in [8] proved that the principal 𝑝th root of the matrix 𝐴 is a multiple of the (2,1)-block of the matrix sign function, sign(𝐶), for the following block companion matrix 𝐶=0𝐼0𝐼⋱⋱⋱𝐼𝐴0∈ℂ^{𝑝𝑛×𝑝𝑛}. (3)

The matrix 𝑆 has the following properties.(1)𝑆²=𝐼 (𝑆 is involutory).(2)𝑆 is diagonalizable with eigenvalues ±1.(3)𝑆𝐴=𝐴𝑆.(4)If 𝐴 is real, then 𝑆 is real.(5)(𝐼+𝑆)/2 and (𝐼−𝑆)/2 are projectors onto the invariant subspaces associated with the eigenvalues in the right half-plane and left half-plane, respectively.

Although 𝑆 has eigenvalues of ±1, its norm can be arbitrarily large. Note that, for diagonalizable 𝐴, eigenvectors of 𝐴 are eigenvectors of 𝑆, with eigenvalues of −1 and 1, respectively. For more, see [9].

There are some other definitions for 𝑆 in the literature based on the Jordan canonical form and the integral representation. As indicated in [6], one of the most useful and widely applicable methods for computing 𝑆 is the matrix iteration of Newton given by 𝑋_𝑘+1=12𝑋_𝑘+𝑋_𝑘⁻¹. (4)

In 1991, a fundamental family of matrix iterations for finding the matrix sign function 𝑆 was introduced in [10] using Padé approximants to 𝑓(𝜉)=(1−𝜉)^−1/2 and the following characterization: (𝑧)𝑧sign=𝑠=𝑧^21/2=𝑧(1−𝜉)^1/2, (5) where 𝜉=1−𝑧². Let the (𝑚,𝑛)-Padé approximant to 𝑓(𝜉) be 𝑃_𝑚,𝑛(𝜉)/𝑄_𝑚,𝑛(𝜉), and 𝑚+𝑛≥1. The iteration 𝑧_𝑘+1=𝑧_𝑘𝑃_𝑚,𝑛1−𝑧²_𝑘𝑄_𝑚,𝑛1−𝑧²_𝑘∶=𝜑_{2𝑚+1,2𝑛} (6) has been proved to be convergent to 1 and −1 with order of convergence 𝑚+𝑛+1 form 𝑚≥𝑛−1. We remark that iterative methods of the type 𝑧_𝑘𝑃_𝑚,𝑛(1−𝑧²_𝑘)/𝑄_𝑚,𝑛(1−𝑧²_𝑘) are fixed-point type iterations, and if 𝑧_𝑘𝑃_𝑚,𝑛(1−𝑧²_𝑘)/𝑄_𝑚,𝑛(1−𝑧²_𝑘) does not converge then its reciprocal; that is, 𝑄_𝑚,𝑛(1−𝑧²_𝑘)/𝑧_𝑘𝑃_𝑚,𝑛(1−𝑧²_𝑘) converges to the sign matrix (e.g., forthcoming iterations (7) and (8)). Generally speaking, the iterations of Kenney and Laub (6), generated by the [ℎ/ℎ] and [(ℎ−1)/ℎ] Padé approximants, are globally convergent and their orders depend on ℎ. A discussion about such iterations was given in [11–13].

A lot of known methods could be extracted from the Padé family (6). For example, the well-known Halley’s matrix iteration of order three can be deduced as follows: 𝑋_𝑘+1=𝐼+3𝑋²_𝑘𝑋_𝑘3𝐼+𝑋²_𝑘⁻¹. (7) Another fourth-order method could be attained as follows: 𝑋_𝑘+1=𝐼+6𝑋²_𝑘+𝑋⁴_𝑘4𝑋_𝑘𝐼+𝑋²_𝑘⁻¹. (8)

Note that, for lower order methods such as (4), the convergence is slow; that is, initially convergence can be slow if |𝑥_𝑘|≫1. Hence, a scaling approach (a.k.a. norm scaling) to accelerate the beginning of this phase is necessary and can be done in what follows [6]: 𝑋₀𝜇=𝐴,_𝑘=𝑋_𝑘⁻¹𝑋_𝑘,𝑋_𝑘+1=12𝜇_𝑘𝑋_𝑘+𝜇_𝑘⁻¹𝑋_𝑘⁻¹. (9)

Such an approach could be done to refine the initial matrix and to provide a much more robust initial matrix to arrive at the convergence phase rapidly.

The rest of this paper has been organized as follows. Section 2 gives the basic idea of obtaining other higher order solvers for computing 𝑆, while a fourth-order family of methods has been introduced. Section 3 is devoted to find the best method of this family in terms of the lowest computational cost. An analysis will be given to show that the new matrix iteration is asymptotically stable with local quartic convergence. To find a method with fourth order and global convergence, we also give another matrix iteration therein. Convergence analysis along with the error bounds of the proposed method is established. Numerical studies will be included in Section 4 to compare the efficiency and the stability of the schemes for finding 𝑆 using different tests. Finally, a short conclusion of the study will be drawn in Section 5.

2. Basic Idea

The connection between the matrix iteration of Newton and Newton’s root-finding method may not be clear at the first sight. Generally speaking, in the theory of matrix functions, many of the matrix functions could effectively be calculated by the existing iterative methods for finding the solution of nonlinear equations [14].

To illustrate further, apply Newton’s method on the following matrix equation: 𝑋²=𝐼, (10) in which 𝐼 is the identity matrix of the appropriate size; it would yield in the matrix Newton’s iteration (4). Note that 𝑆 is one solution of (10). Note that, in the last decade, many efficient higher-order iterative methods have been developed for solving nonlinear equations [15], and some of them have been extended to solve nonlinear matrix equation; see, for example, [16, 17]. But our work is the first to discuss the application of high order root solvers for matrix sign function.

Let us consider the following fourth-order family of iterative methods [18]: 𝑦_𝑘=𝑥_𝑘−𝑓^′𝑥_𝑘⁻¹𝑓𝑥_𝑘,𝑥_𝑘+1=𝑦_𝑘−𝑓^′𝑥_𝑘⁻¹𝑓𝑦_𝑘×1+2𝑓𝑥_𝑘⁻¹𝑓𝑦_𝑘+𝛽𝑓𝑥_𝑘⁻¹𝑓𝑦_𝑘², (11) wherein 𝛽 is a free parameter in ℝ. Applying the uniparametric family (11) on matrix equation (10) results in the following novel family of matrix iterations: 𝑋_𝑘+1=128𝑋⁷_𝑘⁻¹×−𝛽−𝐼+𝑋²_𝑘⁴+8𝑋²_𝑘−5𝑋⁴_𝑘+15𝑋⁶_𝑘+5𝑋⁸_𝑘. (12)

The free parameter 𝛽 plays an important role in the next section to derive the best possible matrix iteration out of (12).

3. Main Results

In order to reduce the computational complexity of (12), the parameter 𝛽 must be chosen as if the number of matrix-matrix products gets down along the number of matrix inversions.

Choosing 𝛽=0 will simplify the whole family into the following method with reasonable computational cost in contrast to its convergence order: 𝑋_𝑘+1=16𝑋⁵_𝑘⁻¹𝐼−5𝑋²_𝑘+15𝑋⁴_𝑘+5𝑋⁶_𝑘. (13)

We now rewrite obtained iteration (13) as efficiently as possible to reduce the number of matrix-matrix multiplications in what follows: 𝑋_𝑘+1=1𝑌16^5𝑘−5𝑌³_𝑘+15𝑌_𝑘+5𝑋_𝑘, (14) where 𝑌_𝑘=𝑋_𝑘⁻¹ and 𝑋₀=𝐴.

Definition 1 (stability [6]). Consider an iteration 𝑋_𝑘+1=𝑔(𝑋_𝑘) with a fixed point 𝑋. Assume that 𝑔 is Fréchet differentiable at 𝑋. The iteration is stable in a neighborhood of 𝑋 if the Fréchet derivative 𝐿_𝑔(𝑋) has bounded powers; that is, there exists a constant 𝑐 such that ‖𝐿^𝑖_𝑔(𝑋)‖≤𝑐 for all 𝑖>0.

Now, we first investigate the stability of (14) for the matrix sign function in a neighborhood of the solution of (10). In fact, we analyze how a small perturbation at the 𝑘th iterate is amplified or damped along the iterates. Note that a general way for assessing the stability of some matrix iterations has been studied by Iannazzo in [19]. The forthcoming approach follows these results of [19].

Lemma 2. The sequence {𝑋_𝑘}^𝑘=∞_𝑘=0 generated by (14) is asymptotically stable.

Proof. Stability concerns behavior close to convergence and so is an asymptotic property. Let Δ𝑋_𝑘 be the numerical perturbation introduced at the 𝑘th iterate of (14). Next, one has 𝑋_𝑘=𝑋_𝑘+Δ𝑋_𝑘. (15) Here, we perform a first-order error analysis; that is, we formally neglect quadratic terms such as (Δ𝑋_𝑘)². This formal manipulation is meaningful if Δ𝑋_𝑘 is sufficiently small. We have 𝑋_𝑘+1=1𝑌16^5𝑘−5𝑌³_𝑘+15𝑌_𝑘+5𝑋_𝑘=1𝑋16^𝑘−5−5𝑋_𝑘⁻³+15𝑋_𝑘⁻¹+5𝑋_𝑘=1𝑋16^𝑘+Δ𝑋_𝑘⁻⁵−5𝑋_𝑘+Δ𝑋_𝑘⁻³+15𝑋_𝑘+Δ𝑋_𝑘⁻¹+5𝑋_𝑘+Δ𝑋_𝑘≈1𝑋16^𝑘−1−𝑋_𝑘⁻¹Δ𝑋_𝑘⋅𝑋_𝑘⁻¹⁵−5𝑋_𝑘⁻¹−𝑋_𝑘⁻¹Δ𝑋_𝑘⋅𝑋_𝑘⁻¹³+15𝑋_𝑘⁻¹−𝑋_𝑘⁻¹Δ𝑋_𝑘⋅𝑋_𝑘⁻¹+5𝑋_𝑘+Δ𝑋_𝑘, (16) where the following identity has been used (for any nonsingular matrix 𝐵 and the matrix 𝐶): (𝐵+𝐶)⁻¹≈𝐵⁻¹−𝐵⁻¹𝐶𝐵⁻¹. (17)
Note that the commutativity between 𝑋_𝑘 and Δ𝑋_𝑘 is not used throughout this paper because it does not hold. Further simplifying yields to 𝑋_𝑘+11≈𝑆+2Δ𝑋_𝑘−12𝑆Δ𝑋_𝑘𝑆, (18) where 𝑆²=𝐼, 𝑆⁻¹=𝑆, and for large enough 𝑘 we have 𝑋_𝑘≈sign(𝐴)=𝑆, and (Δ𝑋_𝑘)^𝑖, 𝑖≥2 is close to zero (matrix) and can be neglected by choosing (Δ𝑋_𝑘)^𝑖≈0. After some algebraic manipulation and using Δ𝑋_𝑘+1=𝑋_𝑘+1−𝑋_𝑘+1, we conclude that Δ𝑋_𝑘+1≈12Δ𝑋_𝑘−12𝑆Δ𝑋_𝑘𝑆. (19) Applying (19) recursively, and after some algebraic manipulations, we have Δ𝑋_𝑘+1≤12^𝑘+1Δ𝑋₀−𝑆Δ𝑋₀𝑆. (20) From (20), we can conclude that the perturbation at the iterate 𝑘+1 is bounded. Therefore, the sequence {𝑋_𝑘} generated by (14) is asymptotically stable.

Remark 3. If 𝑋₀ is a function of 𝐴, then the iterates from (14) or (23) are all functions of 𝐴 and hence commute with 𝐴.

The proposed matrix iteration requires one matrix inverse per iteration along four matrix-matrix products to achieve the convergence order four, while Halley’s method requires one matrix inversion and three matrix-matrix products to reach order three.

The basins of attraction for the two iterative methods of (4) and (14) to solve 𝑥²−1=0 in the complex plane have been drawn in Figure 1. It shows that the new scheme has larger basins due to its higher convergence order. Note that the roots have been identified with two white points.

(a)

(b)

Figure 1 

The basins of attraction for (4) (a) and fourth-order method (14) (b) for the polynomial 𝑥2−1=0 (shaded by the number of iterations to obtain the solution).

Iannazzo in [20] discussed that the matrix convergence is governed by scalar convergence. That is to say, the fourth-order convergence of new method (14) might not be global unlike Newton’s iteration (4). To illustrate further, if one chooses a matrix 𝐴 (in Figure 1(b)) with one eigenvalue with negative real part, but in a yellow petal, then the matrix iteration will not converge to 𝑆. Therefore, it could be mentioned that scheme (14) converges with local fourth order. This restriction has encouraged us to propose another fourth-order method with global convergence in what follows and leave behind the previous iteration with interest in terms of theoretical analysis only.

Let us apply the following fourth-order nonlinear solver [16] on matrix equation (10): 𝑦_𝑘=𝑥_𝑘−2⁻¹𝑓^′𝑥_𝑘⁻¹𝑓𝑥_𝑘,𝑧_𝑘=𝑥_𝑘−𝑓^′𝑦_𝑘⁻¹𝑓𝑥_𝑘,𝑥_𝑘+1=𝑧_𝑘−𝑧_𝑘−𝑦_𝑘⁻¹𝑓𝑧_𝑘−𝑓𝑦_𝑘⁻¹𝑓𝑧_𝑘, (21) which reads the error equation 𝑒_𝑘+1=𝑐³₂2−𝑐₂𝑐₃8𝑒⁴_𝑘+𝑂𝑒⁵_𝑘, (22) wherein 𝑐_𝑗=𝑓^(𝑗)(𝛼)/𝑗!𝑓^′(𝛼) and  𝑒_𝑘=𝑥_𝑘−𝛼, and then to obtain its corresponding matrix iteration, as follows 𝑋_𝑘+1=𝐼+18𝑌_𝑘+13𝑍_𝑘𝑋_𝑘7𝐼+𝑌_𝑘𝐼+3𝑌_𝑘⁻¹, (23) where 𝑌_𝑘=𝑋_𝑘𝑋_𝑘,  𝑍_𝑘=𝑌_𝑘𝑌_𝑘, and 𝑋₀=𝐴. Figure 2 shows the basins of attraction for new method (23) and scheme (7), while both reveal global convergence. The higher order of convergence for (23) made its basins larger and lighter.

(a)

(b)

Figure 2 

The basins of attraction for (7) (a) and method (23) (b) for the polynomial 𝑥2−1=0 (shaded by the number of iterations to obtain the solution).

Remark 4. In this paper, we restrict the analyses to asymptotically small perturbations; that is, we use the differential error analysis.

Lemma 5. The sequence {𝑋_𝑘}^𝑘=∞_𝑘=0 generated by (23) is asymptotically stable.

Proof. Using the same assumptions as in the Proof of Lemma 2 (perturbations are restricted to a neighborhood of the solution), we can write the following: 𝑋_𝑘+1=𝐼+18𝑌_𝑘+13𝑍_𝑘𝑋_𝑘7𝐼+𝑌_𝑘𝐼+3𝑌_𝑘⁻¹=𝐼+18𝑋²_𝑘+13𝑋⁴_𝑘×𝐼+3𝑋_𝑘⁻¹7𝐼+𝑋_𝑘⁻¹𝑋_𝑘⁻¹≈𝐼+18𝑆+Δ𝑋_𝑘𝑆²+13𝑆+Δ𝑋_𝑘𝑆⁴×𝐼+3𝑆+Δ𝑋_𝑘𝑆⁻¹7𝐼+𝑆+Δ𝑋_𝑘𝑆⁻¹×𝑆+Δ𝑋_𝑘𝑆⁻¹. (24) Further simplifying and by considering the terminology of Lemma 2, we attain Δ𝑋_𝑘+1≈𝐼+118𝑆Δ𝑋_𝑘𝑆−158𝑆Δ𝑋_𝑘𝑆−158Δ𝑋_𝑘1≈−2Δ𝑋_𝑘−12𝑆Δ𝑋_𝑘𝑆. (25) It shows that the main proposed iterative scheme (23) is asymptotically stable, since recursively one has ‖Δ𝑋_𝑘+1‖≤(1/2^𝑘+1)‖Δ𝑋₀+𝑆Δ𝑋₀𝑆‖. The proof is complete.

Remark 6. The analysis done does not take into account the influence of the rounding errors on the convergence. Sometimes these errors may lead to slower convergence or even to the failure of the method. We remark concerning this potential danger for problems solving in single precision arithmetic or in lower precision.

Theorem 7. Let 𝐴∈ℂ^𝑛×𝑛 have no pure imaginary eigenvalues. Then, for the proposed iterates {𝑋_𝑘}^𝑘=∞_𝑘=0 in (23), (𝑋_𝑘(7𝐼+𝑌_𝑘)(𝐼+3𝑌_𝑘))⁻¹ is defined per stage.

Proof. We must show that 𝑋_𝑘(7𝐼+𝑌_𝑘)(𝐼+3𝑌_𝑘) which is obtained at each iteration is nonsingular, since the inverse of the matrix 𝑋_𝑘(7𝐼+𝑌_𝑘)(𝐼+3𝑌_𝑘) must be computed per computing step. Toward this goal, it is enough to show that the eigenvalues of the computed matrix at the end of each iteration are in the open half-plane.
Using the initial matrix 𝑋₀=𝐴 and based on the fact that 𝐴 has no eigenvalues on the imaginary axis, the eigenvalues of the initial matrix are in the open half-plane. Let 𝜆 be the eigenvalue of the matrix 𝑋_𝑘 in the 𝑘th iterate. We have 𝜆=𝑟exp𝑖𝜃, where 𝑖=√−1. Hence, 𝜆+𝜆⁻¹=𝑟+𝑟⁻¹cos(𝜃)+𝑖𝑟−𝑟⁻¹(𝜃),𝜆sin₁7+𝜆²1+3𝜆²⁻¹=(cos(𝜃)−𝑖sin(𝜃))³𝑟22𝑟²+7+3𝑟⁴cos(2𝜃)+𝑖−7+3𝑟⁴()=1sin2𝜃7𝑒^𝑖𝜃𝑟+22𝑒^3𝑖𝜃𝑟³+3𝑒^5𝑖𝜃𝑟⁵, (26) and consequently 1+18𝜆²+13𝜆⁴𝜆¹7+𝜆²1+3𝜆²⁻¹=((𝜃)(𝜃×cos−𝑖sin))18𝑟²+1+13𝑟⁴cos(2𝜃)+𝑖−1+13𝑟⁴()×𝑟sin2𝜃7+𝑟²cos(𝜃)+𝑖−7+𝑟²(𝜃)×sin1+3𝑟²cos(𝜃)+𝑖−1+3𝑟²sin(𝜃)⁻¹. (27) Therefore, the eigenvalues of 𝑋_𝑘(7𝐼+𝑌_𝑘)(𝐼+3𝑌_𝑘) remain in their open half-plane under mapping (23). And 𝑋_𝑘(7𝐼+𝑌_𝑘)(𝐼+3𝑌_𝑘) is defined and is nonsingular for all 𝑘.

Theorem 8. Let 𝐴∈ℂ^𝑛×𝑛 has no pure imaginary eigenvalues. Then, the proposed iterates {𝑋_𝑘}^𝑘=∞_𝑘=0 of (23) converge to the sign matrix 𝑆.

Proof. For our analysis, we assume that 𝐴 is diagonalizable; that is, there exists a nonsingular matrix 𝑉 such that 𝑉⁻¹𝐴𝑉=Λ=diag𝜆₁,𝜆₂,…,𝜆_𝑛, (28) where 𝜆₁,𝜆₂,…,𝜆_𝑛 are the eigenvalues of 𝐴. Note that we know that [10] sign(Λ)=sign𝑉⁻¹𝐴𝑉=𝑉⁻¹sign(𝐴)𝑉, (29) for any nonsingular matrix 𝑉. On the other hand, if we define 𝐷_𝑘=𝑉⁻¹𝑋_𝑘𝑉, then we have from (23) that 𝐷_𝑘+1=𝐼+18𝐷²_𝑘+13𝐷⁴_𝑘𝐷_𝑘7𝐼+𝐷²_𝑘𝐼+3𝐷²_𝑘⁻¹. (30) Notice that if 𝐷₀ is a diagonal matrix, then all successive 𝐷_𝑘 are diagonal too. From (30), it is enough to prove that {𝐷_𝑘} converges to the sign of Λ  and then ensure the convergence of the sequence generated by (23) to sign(𝐴).
Therefore, we can write (30) as 𝑛 uncoupled scalar iterations to solve 𝑔(𝑥)=0, with 𝑔(𝑥)=𝑥²−1, given by 𝑑^𝑖_𝑘+1=1+18𝑑^𝑖_𝑘²+13𝑑^𝑖_𝑘⁴𝑑^𝑖_𝑘7𝐼+𝑑^𝑖_𝑘²𝐼+3𝑑^𝑖_𝑘², (31) where 𝑑^𝑖_𝑘=(𝐷_𝑘)_𝑖,𝑖 and 1≤𝑖≤𝑛. On the other hand, sign(𝐷_𝑘)=sign(Λ), for all 𝑘≥0. From (30) and (31), it is enough to study the convergence of {𝑑^𝑖_𝑘} to the sign of 𝜆_𝑖, for all 1≤𝑖≤𝑛.
From (31) and since the eigenvalues of 𝐴 are not pure imaginary, we have that sign(𝜆_𝑖)=𝑠_𝑖=±1. Thus, we attain 𝑑^𝑖_𝑘+1−1𝑑^𝑖_𝑘+1=+1−6+18𝑑^𝑖_𝑘−22𝑑^𝑖_𝑘²+13𝑑^𝑖_𝑘³−3𝑑^𝑖_𝑘⁴8+18𝑑^𝑖_𝑘+22𝑑^𝑖_𝑘²+13𝑑^𝑖_𝑘³+3𝑑^𝑖_𝑘⁴. (32) Since |𝑑^𝑖₀|=|𝜆_𝑖|>0, we obtain lim_𝑘→∞𝑑^𝑖_𝑘+1−1𝑑^𝑖_𝑘+1+1=1, (33) and lim_𝑘→∞|𝑑^𝑖_𝑘|=1=|sign(𝜆_𝑖)|. This shows that {𝑑^𝑖_𝑘} is convergent. Now, it could be easy to conclude that lim_𝑘→∞𝐷_𝑘=sign(Λ). Recalling 𝐷_𝑘=𝑉⁻¹𝑋_𝑘𝑉, we have lim_𝑘→∞𝑋_𝑘=𝑉lim_𝑘→∞𝐷_𝑘𝑉⁻¹=𝑉sign(Λ)𝑉⁻¹, (34) and subsequently the convergence is established. The proof is complete.

Remark 9. It is clear that convergence will be slow if either 𝜌(𝐴)≫1 or 𝐴 has eigenvalues close to the imaginary axis. Hence, it is better to first construct a robust seed by scaled method (9).

Theorem 10. Let 𝐴∈ℂ^𝑛×𝑛 have no pure imaginary eigenvalues. Then, new method (23) has fourth order to find the sign matrix 𝑆.

Proof. Clearly, the 𝑋_𝑘 are rational functions of 𝐴 and hence, like 𝐴, commute with 𝑆. On the other hand, we know that 𝑆²=𝐼, 𝑆⁻¹=𝑆, 𝑆^2𝑗=𝐼, and 𝑆^2𝑗+1=𝑆, 𝑗≥1. Choosing 𝐵_𝑘=𝑋_𝑘(7𝐼+𝑌_𝑘)(𝐼+3𝑌_𝑘) (for the sake of simplicity), we have 𝑋_𝑘+1−𝑆=𝐼+18𝑌_𝑘+13𝑍_𝑘𝐵_𝑘⁻¹=−𝑆𝐼+18𝑋²_𝑘+13𝑋⁴_𝑘−𝑆𝐵_𝑘𝐵_𝑘⁻¹=𝐼+18𝑋²_𝑘+13𝑋⁴_𝑘−7𝑋_𝑘𝑆−22𝑋³_𝑘𝑆−3𝑋⁵_𝑘𝑆𝐵_𝑘⁻¹=𝑋_𝑘−𝑆⁴−3𝑋⁵_𝑘𝑆⁵+12𝑋⁴_𝑘𝑆⁴−18𝑋³_𝑘𝑆+12𝑋²_𝑘𝑆²−𝑋_𝑘𝑆³𝐵_𝑘⁻¹=𝑋_𝑘−𝑆⁴−3𝑋_𝑘𝑆×𝑋⁴−4𝑋³_𝑘𝑆+6𝑋²_𝑘𝑆²−4𝑋_𝑘𝑆³+𝐼𝐵_𝑘⁻¹=𝑋_𝑘−𝑆⁴−3𝑋_𝑘𝑆𝑋_𝑘−𝑆⁴𝐵_𝑘⁻¹=𝑋_𝑘−𝑆⁴𝐼−3𝑋_𝑘𝑆𝐵_𝑘⁻¹. (35) Now, using a matrix operator norm from both sides of (35), we attain 𝑋_𝑘+1−𝑆≤𝐵_𝑘⁻¹𝐼−3𝑋_𝑘𝑆𝑋_𝑘−𝑆⁴. (36) This reveals the fourth order of convergence for new method (23). The proof is ended.

4. Numerical Results

This section addresses issues related to the numerical precision of the computation of matrix sign function using Mathematica 8 built-in precision [21, 22]. The value of machine precision that produced the results included here is 15.96 digits, which corresponds to a 53-digit binary double precision number with a mantissa [23].

For numerical comparisons in this section, we have used methods (4) denoted by “Newton,” (7) denoted by “Halley,” (8) denoted by “M4,” (14) denoted by “PM1,” and (23) denoted by “PM2.”

It must be noted that the Newton-Schulz iteration, which replaces the inverse of the matrix 𝑋_𝑘 in each iteration of (4) by the Schulz inverse-finder [6] and can be written as 𝑋_𝑘+1=12𝑋_𝑘3𝐼−𝑋²_𝑘, (37) will not be considered in the numerical comparisons. Because due to the use of Schulz iteration instead of the matrix inversion, though the quadratic convergence of (4) will remain unchanged, it fully demands a good initial matrix and might be more risky to diverge in contrast to scheme (4). In fact, its convergence is guaranteed only if ‖𝐼−𝐴²‖<1; see Figure 3(a).

(a)

(b)

Figure 3 

The basins of attraction for (37) and (a) for the complex polynomial 𝑥2−1=0 (shaded by the number of iterations to obtain the solution) and the distribution of the eigenvalues of 𝐴[1] for Example 1.

We report the running time using the command AbsoluteTiming[] for the elapsed CPU time (in seconds) in the experiments. The computer specifications are Microsoft Windows XP Intel(R), Pentium(R) 4, and CPU 3.20 GHz, with 4 GB of RAM.

Example 1. The aim of this example is to compare different methods for finding the matrix sign function of a randomly generated dense 600×600 matrix as follows:  n = 600; SeedRandom[22]; A[1] = RandomReal[{-100, 100}, {n, n}].In this test example, the prescribed tolerance is ‖𝑋²_𝑘−𝐼‖_𝐹≤10⁻⁸ and the maximum number of iterations is set to 100.

The results of comparisons in terms of the number of iterations and the computational time have been reported in Table 1, for various matrix iterations in finding the matrix sign function numerically. Note that whatever the eigenvalues of a matrix are closer to the imaginary axis, the speed of convergence for the different methods becomes slower and more risky to face with singular matrices 𝑋_𝑘, whose inverse could not be computed; see Figure 3(b).

Table 1 

Results of comparisons for Example 1.

Note that method (13) is not competitive in terms of computational cost and local convergence and hence we will remove it from further consideration. Unlike it, new method (23) has global convergence with asymptotical stability and could be considered as an alternative over the existing iterative methods for finding 𝑆.

Example 2. In this example, 15 random dense 120×120 matrices are considered to compare the behavior of different methods in what follows:  n = 120; number = 15; SeedRandom[1]; Table[A[l] = RandomReal[{-5, 5},{n, n}];, {l, number}].For this test, the prescribed tolerance is ‖𝑋_𝑘+1−𝑋_𝑘‖₂≤10⁻¹² and the maximum number of iterations is set to 100.

The results of comparisons are reported in Figure 4. New method (23) beats its competitors in terms of the number of iterations, while both (23) and (8) are the best iterations in terms of the computational time. Note that, in the last two examples, we have used double precision arithmetic in our calculations.

(a)

(b)

Figure 4 

The comparisons of different matrix iterations in terms of the number of iterations and the computational time for 15 different test matrices in Example 2.

Although this example showed the robustness of new method (23), there is an approach to observe the order of convergence of different iteration methods numerically. To be more precise, the computational order of convergence for matrix iterations in finding the matrix sign function can be estimated by 𝑋COC=log²_𝑘+1−𝐼/𝑋²_𝑘−𝐼log𝑋²_𝑘−𝐼/𝑋²_𝑘−1,−𝐼 (38) wherein 𝑋_𝑘−1,𝑋_𝑘,𝑋_𝑘+1 are the last three approximations for finding 𝑆 in the convergence phase.

Example 3 (an academical test). Let us find the computational order of convergence for different methods when finding the matrix sign for the well-known Wilson matrix as follows: 𝐴=.1078775658610975910 (39) In order to find the COC, we herein apply 64-digit fixed point arithmetic in our calculations.

The convergence history along the COCs (in the infinity norm) using formula (38) for different methods is illustrated in Figure 5 and Table 2, applying the stopping termination ‖𝑋²_𝑘−𝐼‖_∞≤10⁻¹⁶. Results show that new method (23) is quite fast and its computational order of convergence for academical tests in high precision computing environment is around 4.

Table 2 

Results of comparisons for Example 3.

Figure 5 

Convergence history based on the logarithm of the residuals ‖𝑋2𝑘−𝐼‖∞ in Example 3.

4.1. Scaling

Main proposed iteration (23) is quite fast and reliable due to the discussions in Sections 3 and 4. However, a way is open for speeding up its initial phase of convergence.

An effective way to enhance the initial speed of convergence is to scale the iterates prior to each iteration; that is, 𝑋_𝑘 is replaced by 𝜇_𝑘𝑋_𝑘. Such an approach can simply be done in what follows: 𝑋₀𝜇=𝐴,_𝑘=𝑋_𝑘⁻¹𝑋_𝑘,𝑋_𝑘+1=𝐼+18𝜇²_𝑘𝑌_𝑘+13𝜇⁴_𝑘𝑍_𝑘×𝜇_𝑘𝑋_𝑘7𝐼+𝜇²_𝑘𝑌_𝑘𝐼+3𝜇²_𝑘𝑌_𝑘⁻¹, (40) where lim_𝑘→∞𝜇_𝑘=1 and lim_𝑘→∞𝑋_𝑘=𝑆.

5. Discussion

A function of a matrix can be defined and computed in several ways, such as Cauchy integral, polynomial interpolation, and Jordan canonical form. However, another approach is to use iteration methods for such computations. Several matrix functions can be computed by iteration 𝑋_𝑘+1=𝑔(𝑋_𝑘), with an appropriate initial matrix 𝑋₀ where, for reasons of computational cost, 𝑔 is usually a polynomial or a rational function.

Under this motivation, in this paper we have introduced and demonstrated some fourth-order matrix methods for finding the matrix sign function. The proposed methods consist of one matrix inversion per cycle and are asymptotically stable. The consistency and efficiency of the contributed methods have also been tested numerically for finding the matrix sign functions to support the theoretical parts.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors are grateful to the anonymous referees for their valuable comments and suggestions for improving the paper.