On a Cubically Convergent Iterative Method for Matrix Sign

Abstract

We propose an iterative method for finding matrix sign function. It is shown that the scheme has global behavior with cubical rate of convergence. Examples are included to show the applicability and efficiency of the proposed scheme and its reciprocal.

1. Introduction

It is known that the function of sign in the scalar case is defined for any not on the imaginary axis by An extension of (1) for the matrix case was given firstly by Roberts in [1]. This extended matrix function is of clear importance in several applications (see, e.g., [2] and the references therein).

Assume that is a matrix with no eigenvalues on the imaginary axis. To define this matrix function formally, let be a Jordan canonical form arranged so that , where the eigenvalues of lie in the open left half-plane and those of lie in the open right half-plane; then where . A simplified definition of the matrix sign function for Hermitian case (eigenvalues are all real) is where is a diagonalization of .

The importance of computing is also due to the fact that the sign function plays a fundamental role in iterative methods for matrix roots and the polar decomposition [3].

Note that although is a square root of the identity matrix, it is not equal to or unless the spectrum of lies entirely in the open right half-plane or open left half-plane, respectively. Hence, in general, is a nonprimary square root of .

In this paper, we focus on iterative methods for finding . In fact, such methods are Newton-type schemes which are in essence fixed-point-type methods by producing a convergent sequence of matrices via applying a suitable initial matrix.

The most famous method of this class is the quadratic Newton method defined by

It should be remarked that iterative methods, such as (6), and the Newton-Schultz iteration or the cubically convergent Halley method are all special cases of the Padé family proposed originally in [4]. The Padé approximation belongs to a broader category of rational approximations. Coincidentally, the best uniform approximation of the sign function on a pair of symmetric but disjoint intervals can be expressed as a rational function.

Note that although (7) does not possess a global convergence behavior, on state-of-the-art parallel computer architectures, matrix inversions scale less satisfactorily than matrix multiplications do, and subsequently (7) is useful in some problems. However, due to local convergence behavior, it is excluded from our numerical examples in this work.

The rest of this paper is organized as follows. In Section 2, we discuss how to construct a new iterative method for finding (3). It is also shown that the constructed method is convergent with cubical rate. It is noted that its reciprocal iteration obtained from our main method is also convergent. Numerical examples are furnished to show the higher numerical accuracy for the constructed solvers in Section 3. The paper ends in Section 4 with some concluding comments.

2. A New Method

The connection of matrix iteration methods with the sign function is not immediately obvious, but in fact such methods can be derived by applying a suitable root-finding method to the nonlinear matrix equation and when of course is one solution of this equation (see for more [5]).

Here, we consider the following root-solver: with . In what follows, we observe that (10) possesses third order of convergence.

Theorem 1. Let be a simple zero of a sufficiently differentiable function , which contains as an initial approximation. Then the iterative expression (10) satisfies where , .

Proof. The proof would be similar to the proofs given in [6].

Applying (10) on the matrix equation (9) will result in the following new matrix fixed-point-type iteration for finding (3): where . This is named PM1 from now on.

The proposed scheme (12) is not a member of Padé family [4]. Furthermore, applying (10) on the scalar equation provides a global convergence in the complex plane (except the points lying on the imaginary axis). This global behavior, which is kept for matrix case, has been illustrated in Figure 1 by drawing the basins of attraction for (6) and (8). The attraction basins for (7) (local convergence) and (12) (global convergence) are also portrayed in Figure 2.

(a)

(b)

Figure 1 

Attraction basins for (6) (a) and (8) (b) for the polynomial .

(a)

(b)

Figure 2 

Attraction basins of (7) (a) and (12) (b) for the polynomial .

Theorem 2. Let have no pure imaginary eigenvalues. Then, the matrix sequence defined by (12) converges to , choosing .

Proof. We remark that all matrices, whether they are diagonalizable or not, have a Jordan normal form , where the matrix consists of Jordan blocks. For this reason, let have a Jordan canonical form arranged as where is a nonsingular matrix and are square Jordan blocks corresponding to eigenvalues lying in and , respectively. We have If we define , then, from the method (12), we obtain Note that if is a diagonal matrix, then, based on an inductive proof, all successive are diagonal too. From (15), it is enough to show that converges to . We remark that the case at which is not diagonal will be discussed later in the proof.
In the meantime, we can write (15) as uncoupled scalar iterations to solve , given by where and . From (15) and (16), it is enough to study the convergence of to .
It is known that . Thus, we attain Since , we have and . This shows that is convergent.
In the convergence proof, may not be diagonal. Since the Jordan canonical form of some matrices may not be diagonal, thus, one cannot write (15) as uncoupled scalar iterations (16). We comment that in this case our method is also convergent. To this goal, we must pursue the scalar relationship among the eigenvalues of the iterates for the studied rational matrix iteration.
In this case, the eigenvalues of are mapped from the iterate to the iterate by the following relation: So, (19) clearly shows that the eigenvalues in the general case are convergent to ; that is to say, Consequently, we have The proof is ended.

Theorem 3. Let have no pure imaginary eigenvalues. Then the proposed method (12) converges cubically to the sign matrix .

Proof. Clearly, are rational functions of and, hence, like , commute with . On the other hand, we know that , , , and , . Using the replacement , we have Now, using any matrix norm from both sides of (22), we attain This reveals the cubical rate of convergence for the new method (12). The proof is complete.

It should be remarked that the reciprocal iteration obtained from (12) is also convergent to the sign matrix (3) as follows: where . This is named PM2. Similar convergence results as the ones given in Theorems 2-3 hold for (24).

A scaling approach to accelerate the beginning phase of convergence is normally necessary since the convergence rate cannot be seen in the initial iterates. Such an idea was discussed fully in [7] for Newton’s method. An effective way to enhance the initial speed of convergence is to scale the iterates prior to each iteration; that is, is replaced by . Subsequently, we can present the accelerated forms of our proposed methods as follows: or where and . The different scaling factors for in (27) are borrowed from Newton’s method. For this reason it is important to show the behavior of the accelerator methods (25)-(26) and this will be done in the next section.

3. Numerical Examples

In this section, the results of comparisons in terms of number of iterations and the residual norms have been reported for various matrix iterations. We compare PM1 and PM2 with (6) denoted by NM and (8) denoted by HM. The programming package Mathematica [8] is used throughout this section. In Tables 1 and 2, IT stands for the number of iterates.

Table 1 

Results of comparisons for Example 5 using .

Table 2 

Results of comparisons for Example 6 using .

Note that the computational order of convergence for matrix iterations in finding can be estimated by [9] where , , and are the last three approximations.

Example 4. In this example, we compare the methods for the following complex matrix:
n = 500; SeedRandom;
A = RandomComplex[{-100 - I, 100 + I},{n,n}];
We apply here double precision arithmetic with the stop termination . Results are given in Figure 3.

Figure 3 

Convergence history versus number of iterations for different methods in Example 4.

Example 5 (academic test). We compute the matrix sign for the following complex test problem: where We apply here 600-digit fixed point arithmetic in our calculations with the stop termination . The results for this example are illustrated in Table 1. We report the COCs in .

Iterative schemes PM1 and PM2 are evidently believed to be more favorable than the other compared methods due to their fewer number of iterations and acceptable accuracy. Hence, the proposed methods with properly chosen initial matrix can be helpful in finding the sign of a nonsingular complex matrix.

Example 6. Here we rerun Example 5 using the scaling approaches (27) with the stop termination . The results for this example are illustrated in Table 2. We used the determinantal scaling for all compared methods. The numerical results uphold the theoretical discussions of Section 2.

A price paid for the high order convergence is the increased amount of matrix multiplications and inversions. This is a typical consequence. However the most important advantage of the presented methods in contrast to the methods of the same orders, such as (8), is their larger attraction basins. This superiority basically allows the new methods to converge to a required tolerance in one lower iteration than their same order methods. Hence, studying the thorough computational efficiency index of the proposed methods may not be an easy task and it must be pursued experimentally. In an experimental manner, if the costs of one matrix-matrix product and one matrix inversion are unity and 1.5 of unity, respectively, then we have the following efficiency indices for different methods: , , and . Note that for Newton’s method we have one matrix-matrix product per cycle due to the computation of stopping criterion. Other similar computations for efficiency indices for different examples show similar behaviors to the above mentioned one.

4. Summary

Matrix functions are used in many areas of linear algebra and arise in numerous applications in science and engineering. The function of a matrix can be defined in several ways, of which the following three are generally the most useful: Jordan canonical form, polynomial interpolation, and finally Cauchy integral.

In this paper, we have focus on iterative methods for this purpose. Hence, a third order nonlinear equation solver has been employed for constructing a new method for . It was shown that the convergence is global via attraction basins in the complex plane and the rate of convergence is cubic. Furthermore, PM2 as the reciprocal of the method PM1 with the same convergence properties was proposed. The acceleration of PM1 and PM2 via scaling was also illustrated simply.

Finally some numerical examples in both double and multiple precisions were performed to show the efficiency of PM1 and PM2. Further researches must be forced to extend the obtained iterations for computing polar decompositions in future studies.

Conflict of Interests

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

Acknowledgment

The authors would like to thank the referees for their helpful corrections and suggestions.