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.

M. Sharifi 1 and S. Karimi Vanani 1 and F. Khaksar Haghani 1 and M. Arab 1 and S. Shateyi 2

Academic Editor:Predrag S. Stanimirovic

1, Department of Mathematics, Islamic Azad University, Shahrekord Branch, Shahrekord, Iran
2, Department of Mathematics and Applied Mathematics, University of Venda, Thohoyandou 0950, South Africa

Received 20 July 2014; Revised 29 August 2014; Accepted 30 August 2014; 29 March 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

It is known that the function of sign in the scalar case is defined for any [figure omitted; refer to PDF] not on the imaginary axis by [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] is a matrix with no eigenvalues on the imaginary axis. To define this matrix function formally, let [figure omitted; refer to PDF] be a Jordan canonical form arranged so that [figure omitted; refer to PDF] , where the eigenvalues of [figure omitted; refer to PDF] lie in the open left half-plane and those of [figure omitted; refer to PDF] lie in the open right half-plane; then [figure omitted; refer to PDF] where [figure omitted; refer to PDF] . A simplified definition of the matrix sign function for Hermitian case (eigenvalues are all real) is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is a diagonalization of [figure omitted; refer to PDF] .

The importance of computing [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] is a square root of the identity matrix, it is not equal to [figure omitted; refer to PDF] or [figure omitted; refer to PDF] unless the spectrum of [figure omitted; refer to PDF] lies entirely in the open right half-plane or open left half-plane, respectively. Hence, in general, [figure omitted; refer to PDF] is a nonprimary square root of [figure omitted; refer to PDF] .

In this paper, we focus on iterative methods for finding [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF]

It should be remarked that iterative methods, such as (6), and the Newton-Schultz iteration [figure omitted; refer to PDF] or the cubically convergent Halley method [figure omitted; refer to PDF] are all special cases of the Pade family proposed originally in [4]. The Pade 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 [figure omitted; refer to PDF] and when of course [figure omitted; refer to PDF] is one solution of this equation (see for more [5]).

Here, we consider the following root-solver: [figure omitted; refer to PDF] with [figure omitted; refer to PDF] . In what follows, we observe that (10) possesses third order of convergence.

Theorem 1.

Let [figure omitted; refer to PDF] be a simple zero of a sufficiently differentiable function [figure omitted; refer to PDF] , which contains [figure omitted; refer to PDF] as an initial approximation. Then the iterative expression (10) satisfies [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] .

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): [figure omitted; refer to PDF] where [figure omitted; refer to PDF] . This is named PM1 from now on.

The proposed scheme (12) is not a member of Pade family [4]. Furthermore, applying (10) on the scalar equation [figure omitted; refer to PDF] 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.

Figure 1: Attraction basins for (6) (a) and (8) (b) for the polynomial [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Figure 2: Attraction basins of (7) (a) and (12) (b) for the polynomial [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Theorem 2.

Let [figure omitted; refer to PDF] have no pure imaginary eigenvalues. Then, the matrix sequence [figure omitted; refer to PDF] defined by (12) converges to [figure omitted; refer to PDF] , choosing [figure omitted; refer to PDF] .

Proof.

We remark that all matrices, whether they are diagonalizable or not, have a Jordan normal form [figure omitted; refer to PDF] , where the matrix [figure omitted; refer to PDF] consists of Jordan blocks. For this reason, let [figure omitted; refer to PDF] have a Jordan canonical form arranged as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is a nonsingular matrix and [figure omitted; refer to PDF] are square Jordan blocks corresponding to eigenvalues lying in [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , respectively. We have [figure omitted; refer to PDF] If we define [figure omitted; refer to PDF] , then, from the method (12), we obtain [figure omitted; refer to PDF] Note that if [figure omitted; refer to PDF] is a diagonal matrix, then, based on an inductive proof, all successive [figure omitted; refer to PDF] are diagonal too. From (15), it is enough to show that [figure omitted; refer to PDF] converges to [figure omitted; refer to PDF] . We remark that the case at which [figure omitted; refer to PDF] is not diagonal will be discussed later in the proof.

In the meantime, we can write (15) as [figure omitted; refer to PDF] uncoupled scalar iterations to solve [figure omitted; refer to PDF] , given by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . From (15) and (16), it is enough to study the convergence of [figure omitted; refer to PDF] to [figure omitted; refer to PDF] .

It is known that [figure omitted; refer to PDF] . Thus, we attain [figure omitted; refer to PDF] Since [figure omitted; refer to PDF] , we have [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . This shows that [figure omitted; refer to PDF] is convergent.

In the convergence proof, [figure omitted; refer to PDF] may not be diagonal. Since the Jordan canonical form of some matrices may not be diagonal, thus, one cannot write (15) as [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] are mapped from the iterate [figure omitted; refer to PDF] to the iterate [figure omitted; refer to PDF] by the following relation: [figure omitted; refer to PDF] So, (19) clearly shows that the eigenvalues in the general case are convergent to [figure omitted; refer to PDF] ; that is to say, [figure omitted; refer to PDF] Consequently, we have [figure omitted; refer to PDF] The proof is ended.

Theorem 3.

Let [figure omitted; refer to PDF] have no pure imaginary eigenvalues. Then the proposed method (12) converges cubically to the sign matrix [figure omitted; refer to PDF] .

Proof.

Clearly, [figure omitted; refer to PDF] are rational functions of [figure omitted; refer to PDF] and, hence, like [figure omitted; refer to PDF] , commute with [figure omitted; refer to PDF] . On the other hand, we know that [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . Using the replacement [figure omitted; refer to PDF] , we have [figure omitted; refer to PDF] Now, using any matrix norm from both sides of (22), we attain [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] . 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, [figure omitted; refer to PDF] is replaced by [figure omitted; refer to PDF] . Subsequently, we can present the accelerated forms of our proposed methods as follows: [figure omitted; refer to PDF] or [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . The different scaling factors for [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] .

Methods	NM	HM	PM1	PM2
IT	14	9	8	8
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	1.99077	3	3	3

Table 2: Results of comparisons for Example 6 using [figure omitted; refer to PDF] .

Methods	NM	HM	PM1	PM2
IT	10	7	6	6
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	2.00228	3.00001	3.00015	3

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

Example 4.

In this example, we compare the methods for the following [figure omitted; refer to PDF] complex matrix:

n = 500; SeedRandom [figure omitted; refer to PDF] ;

A = RandomComplex[{-100 - I, 100 + I},{n,n}];

We apply here double precision arithmetic with the stop termination [figure omitted; refer to PDF] . Results are given in Figure 3.

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

[figure omitted; refer to PDF]

Example 5 (academic test).

We compute the matrix sign for the following complex test problem: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] We apply here 600-digit fixed point arithmetic in our calculations with the stop termination [figure omitted; refer to PDF] . The results for this example are illustrated in Table 1. We report the COCs in [figure omitted; refer to PDF] .

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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] . 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: [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] . 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.

Acknowledgment

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

Conflict of Interests

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

References

[1] J. D. Roberts, "Linear model reduction and solution of the algebraic Riccati equation by use of the sign function," International Journal of Control , vol. 32, no. 4, pp. 677-687, 1980.

[2] C. S. Kenney, A. J. Laub, "The matrix sign function," IEEE Transactions on Automatic Control , vol. 40, no. 8, pp. 1330-1348, 1995.

[3] N. J. Higham Functions of Matrices: Theory and Computation , Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 2008.

[4] C. Kenney, A. J. Laub, "Rational iterative methods for the matrix sign function," SIAM Journal on Matrix Analysis and Applications , vol. 12, no. 2, pp. 273-291, 1991.

[5] F. Soleymani, P. S. Stanimirovic, S. Shateyi, F. K. Haghani, "Approximating the matrix sign function using a novel iterative method," Abstract and Applied Analysis , vol. 2014, 2014.

[6] F. Soleymani, "Some high-order iterative methods for finding all the real zeros," Thai Journal of Mathematics , vol. 12, no. 2, pp. 313-327, 2014.

[7] C. Kenney, A. J. Laub, "On scaling Newton's method for polar decomposition and the matrix sign function," SIAM Journal on Matrix Analysis and Applications , vol. 13, no. 3, pp. 698-706, 1992.

[8] S. Wagon Mathematica in Action , Springer, New York, NY, USA, 2010., 3rd.

[9] F. Soleymani, E. Tohidi, S. Shateyi, F. Haghani, "Some matrix iterations for computing matrix sign function," Journal of Applied Mathematics , vol. 2014, 2014.