# Approximating the Matrix Sign Function Using a Novel Iterative

Authors: F. Soleymani; P. S. Stanimirović; S. Shateyi; F. Khaksar Haghani
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This study presents a matrix iterative method for finding the sign of a square complex matrix. It is
shown that the sequence of iterates converges to the sign and has asymptotical stability, provided that the
initial matrix is appropriately chosen. Some illustrations are presented to support the theory.


---

## Body

## 1. Introduction

The generic matrix functionf(A) of a given matrix A∈Cn×n is defined formally by the integral representation
(1)f(A)=12πi∮γ‍f(ζ)(ζI-A)-1dζ,
where f:Ω→C is an analytic function, Ω⊆C, and γ is a closed curve which encircles all eigenvalues of A (it should be contained in the domain of analyticity of f). The integral representation (1) is known as the Cauchy integral formula [1]. The integral of a matrix M should be understood as the matrix whose entries are the integrals of the entries of M. However, this mathematically appealing formula for computing the matrix functions is complicated and needs complex analysis to be fully understandable. Hence, several important other strategies for computing the matrix functions have been proposed and investigated, such as the Jordan canonical form and iterative methods for applied numerical problems (see, e.g., [1–3]).In 1971, Roberts in [4] introduced the matrix sign function as a tool for model reduction and for solving Lyapunov and algebraic Riccati equations. He defined the sign function as a Cauchy integral and obtained the following integral representation of sign⁡(A):
(2)sign⁡(A)=S=2π∫0∞‍(t2I+A2)-1dt.The matrix sign function is widely exploited in numerical linear algebra, especially in the computation of invariant subspaces and solutions of Riccati equations [5–7]. Note that the application of this function enables a matrix to be decomposed into two components whose spectra lie on opposite sides of the imaginary axis. The matrix sign function is a valuable tool for the numerical solution of Sylvester and Lyapunov matrix equations (see, e.g., [8]). An application of a generalization of the Newton iteration for the matrix sign function to the solution of the generalized algebraic Bernoulli equations was considered in [9].Another application of this matrix function as a simple and direct method to derive some fundamental results in the theory of surface waves in anisotropic materials was presented in [10]. The authors of paper [11] investigated some practical iterations for matrix sector function which is a generalization of the matrix sign function.Due to the applicability of the matrix sign function along with the difficulty of representation (2), stable iterative methods have become some viable choices.The most important general family of matrix iterations for finding the matrix sign functionS was introduced in [12] using Padé approximants to f(ξ)=(1-ξ)-1/2 and the following characterization:
(3)sign⁡(z)=s=z(z2)1/2=z(1-ξ)1/2,
where ξ=1-z2 and ξ is less than 1 in magnitude. Let the (m,n)-Padé approximant to f(ξ) be Pm,n(ξ)/Qm,n(ξ) and m+n≥1. The iteration
(4)zk+1=zkPm,n(1-zk2)Qm,n(1-zk2)∶=φ2m+1,2n
has been proved to be convergent to 1 and -1 with the order of convergence m+n+1 for m≥n-1. Generally speaking, the iterations of Kenny and Laub (4), generated by the [m/m] and [(m-1)/m] Padé approximants, are globally convergent. A list of different iterations (in the scalar form) is given in Table 1.Table 1
Iteration functions for finding the sign (in the scalar form) from the Padé family (4).
n
=
0
n
=
1
n
=
2
n
=
3
n
=
4
m
=
0
x
2
x
1
+
x
2
8
x
3
+
6
x
2
-
x
4
16
x
5
+
15
x
2
-
5
x
4
+
x
6
128
x
35
+
140
x
2
-
70
x
4
+
28
x
6
-
5
x
8
m
=
1
x
2
(
3
-
x
2
)
x
(
3
+
x
2
)
1
+
3
x
2
4
x
(
1
+
x
2
)
1
+
6
x
2
+
x
4
-
8
x
(
3
+
5
x
2
)
-
5
-
45
x
2
-
15
x
4
+
x
6
64
x
(
3
+
7
x
2
)
35
+
420
x
2
+
210
x
4
-
28
x
6
+
3
x
8
m
=
2
x
8
(
15
-
10
x
2
+
3
x
4
)
x
4
15
+
10
x
2
-
x
4
1
+
5
x
2
x
(
5
+
10
x
2
+
x
4
)
1
+
10
x
2
+
5
x
4
6
x
+
20
x
3
+
6
x
5
1
+
15
x
2
+
15
x
4
+
x
6
-
16
x
(
3
+
7
x
2
(
2
+
x
2
)
)
-
7
-
140
x
2
-
210
x
4
-
28
x
6
+
x
8
m
=
3
-
x
(
-
35
+
35
x
2
-
21
x
4
+
5
x
6
)
16
x
(
35
+
35
x
2
-
7
x
4
+
x
6
)
8
+
56
x
2
35
x
+
105
x
3
+
21
x
5
-
x
7
6
+
84
x
2
+
70
x
4
x
(
7
+
35
x
2
+
21
x
4
+
x
6
)
1
+
7
x
2
(
3
+
5
x
2
+
x
4
)
8
(
x
+
7
x
3
+
7
x
5
+
x
7
)
1
+
28
x
2
+
70
x
4
+
28
x
6
+
x
8
m
=
4
x
(
315
-
420
x
2
+
378
x
4
-
180
x
6
+
35
x
8
)
128
315
x
+
420
x
3
-
126
x
5
+
36
x
7
-
5
x
9
64
+
576
x
2
105
x
+
420
x
3
+
126
x
5
-
12
x
7
+
x
9
16
+
288
x
2
+
336
x
4
63
x
+
420
x
3
+
378
x
5
+
36
x
7
-
x
9
8
(
1
+
3
x
2
(
9
+
7
x
2
(
3
+
x
2
)
)
)
x
(
3
+
x
2
)
(
3
+
27
x
2
+
33
x
4
+
x
6
)
1
+
36
x
2
+
126
x
4
+
84
x
6
+
9
x
8Note that Iannazzo in [13] pointed out that these iterations can be obtained from the general König family (which goes back to Schröder [14, 15]) applied to the equation z2-1=0. For a recent method in this area, one may refer to [16].A lot of known methods could be extracted from the Padé family (4). For example, the Newton’s iteration can be deduced as the reciprocal of the iteration corresponding to the case m=0 and n=1 in Table 1:
(5)Xk+1=12(Xk+Xk-1).Choosingm=n=2 yields the following fifth order method:
(6)Xk+1=[I+5Xk2(2I+Xk2)][Xk(5I+10Xk2+Xk4)]-1.Similarly, optionsm=1 and n=3, or m=3 and n=1, or m=0 and n=4 result in the following fifth order methods, respectively:
(7)Xk+1=-(-5I-45Xk2-15Xk4+Xk6)[8Xk(3I+5Xk2)]-1,(8)Xk+1=(8I+56Xk2)[Xk(35I+35Xk2-7Xk4+Xk6)]-1,(9)Xk+1=(35I+140Xk2-70Xk4+28Xk6-5Xk8)[128Xk]-1.The remaining sections of this work are organized as follows. Section2 presents the construction and the derivation of a new matrix iteration for finding S using the followingnew nonlinear equation solver:
(10)yk=xk-23f(xk)f′(xk),zk=xk-123f′(yk)+f′(xk)3f′(yk)-f′(xk)f(xk)f′(xk),xk+1=zk-f(zk)f[yk,zk],
where f[yk,zk]=(f(zk)-f(yk))/(zk-yk). (10) is a combination of Jarratt’s method [17] and a secant approach [18].Note that (10) is anovel three-step iterative method (for the scalar case). Section 2 also studies the stability of the method and verifies its asymptotical stability. In Section 3, we discuss some other aspects of the new method, applicable in the implementation. Therein, we derive a new inversion-free method as well as a scaled method for finding S. Section 4 is devoted to the numerical examples for illustrating the convergence behavior of the new method against the existing ones. We emphasize that our constructed solver possesses the same global convergence rate as (6), but numerical example will reveal that it has an equal or faster behavior for solving many randomly generated matrices. This would be a clear advantage of our solver in contrast to (6). Finally, concluding remarks will be drawn in Section 5.


## 2. A Novel Iterative Method

Throughout this work it is assumed thatA∈Cn×n has no eigenvalues on the imaginary axis, so that sign⁡(A) is defined. Note that this assumption implies that A is nonsingular.As discussed in the fundamental article of Kenny and Laub in 1995 [19], the construction of new matrix methods for the matrix sign S is related to the iterative methods for finding the solution of nonlinear scalar equations. Let us apply (10) to the following nonlinear matrix equation:
(11)X2=I,
in which I is the identity matrix. In fact, the sign S is a solution of the matrix equation (11). After application and simplification of its reciprocal, we obtain
(12)Xk+1=(7Xk+30Xk3+11Xk5)[I+20Xk2+25Xk4+2Xk6]-1,
where X0=A.Iannazzo in [20] mentioned that a matrix convergence is governed by the corresponding scalar convergence. Since the method (10) reads the following error equation:
(13)ek+1=13(c24-c22c3)ek5+O(ek6),
wherein cj=f(j)(α)/(j!f′(α)) and ek=xk-α, therefore its corresponding matrix method (12) converges with fifth order of convergence. But the question is whether this convergence is local or global. To answer this question, we draw the basins of attraction for the new scheme along with the existing methods of various orders.It is shown in Figures1–3 that methods (5), (6), and (12) are globally convergent, while methods (7), (8), and (9) are locally convergent (if one chooses a matrix A (in Figure 2 (left)) with one eigenvalue with negative real part, but in a green petal, then the matrix iteration will not converge to S).Figure 1
The basins of attraction for (5) (left) and (6) (right), for the polynomial x2-1=0 (shaded by the number of iterations to obtain the solution).Figure 2
The basins of attraction for (7) (left) and (8) (right), for the polynomial x2-1=0 (shaded by the number of iterations to obtain the solution).Figure 3
The basins of attraction for (9) (left) and (12) (right), for the polynomial x2-1=0 (shaded by the number of iterations to obtain the solution).The higher order convergence of (12) made its basins larger and lighter in contrast to (5). Hence, the new method could be of interest due to its global fifth order of convergence for finding S.Now, an important challenge, that must be proven, is to show the stability of the new method (12) for finding the matrix sign function. This will be done formally in what follows.Definition 1 (stability, see [1]).
Consider an iterationXk+1=g(Xk) with a fixed point X. Assume that g is Fréchet differentiable at X. The iteration is stable in a neighborhood of X if the Fréchet derivative Lg(X) has bounded powers; that is, there exists a constant c such that ∥Lgi(X)∥≤c for all i>0.We investigate the stability of (12) for finding S in a neighborhood of the solution of (11). In fact, we analyze how a small perturbation at the kth iterate is amplified or damped along the iterates. Stability concerns behavior close to convergence and so is an asymptotic property.Lemma 2.
The sequence{Xk}k=0k=∞ generated by (12) is asymptotically stable.Proof.
IfX0 is a function of A, then the iterates from (12) are all functions of A and hence commute with A. Commutativity properties are frequently used when deriving a matrix iteration for finding S.
In this study, we restrict the analysis to asymptotically small perturbations; that is, we use the differential error analysis.
LetΔXk be the numerical perturbation introduced at the kth iterate of (12). Next, one has
(14)X~k=Xk+ΔXk.
Here, we perform a first order error analysis; that is, we formally use approximations (ΔXk)i≈0, since (ΔXk)i, i≥2, is close to the zero (matrix). This formal manipulation is meaningful if ΔXk is sufficiently small. We have
(15)X~k+1=(7X~k+30X~k3+11X~k5)[I+20X~k2+25X~k4+2X~k6]-1=(7(Xk+ΔXk)+30(Xk+ΔXk)3+11(Xk+ΔXk)5)×[I+20(Xk+ΔXk)2+25(Xk+ΔXk)4+2(Xk+ΔXk)6]-1≈(S+10048ΔXk+52100SΔXkS)×(I-1912SΔXk-1912ΔXkS)≈S-2448ΔXk+2448SΔXkS,
where the following identities are used (for any nonsingular matrix B and the matrix C):
(16)(B+C)-1≈B-1-B-1CB-1.
Note that after some algebraic manipulations and usingΔXk+1=X~k+1-Xk+1, we have (assuming Xk≈sign⁡(A)=S for enough large k)
(17)ΔXk+1≈-12ΔXk+12SΔXkS.
We used the following facts on the matrix sign function S2=I, and S-1=S. We can now conclude that the perturbation at the iterate k+1 is bounded; that is,
(18)∥ΔXk+1∥≤12k+1∥SΔX0S-ΔX0∥.
Therefore, the sequence {Xk}k=0k=∞ generated by (12) is asymptotically stable. This ends the proof.Theorem 3.
LetA∈Cn×n have no pure imaginary eigenvalues. Then, the matrix sequence {Xk}k=0k=∞ defined by (12) converges to the matrix sign S.Proof .
LetA have a Jordan canonical form arranged as
(19)V-1AV=Λ=[C00N],
where V is a nonsingular matrix and C, N are square Jordan blocks corresponding to eigenvalues lying in C- (open left-half complex plane) and C+ (open right-half complex plane), respectively. Denote by λ1,…λp and λp+1,…λn values lying on the main diagonals of blocks C and N, respectively.
SinceV is invertible, we know that sign⁡(A) is diagonalizable and [1]
(20)sign⁡(A)=V[-Ip00In-p]V-1.
Therefore,
(21)sign⁡(Λ)=sign⁡(V-1AV)=V-1sign⁡(A)V=diag⁡(sign⁡(λ1),…,sign⁡(λp),1=diag⁡×sign⁡(λp+1),…,sign⁡(λn)).
On the other hand, if we defineDk=V-1XkV, then from the equation (12) we obtain
(22)Dk+1=(7Dk+30Dk3+11Dk5)[I+20Dk2+25Dk4+2Dk6]-1.
Notice that if D0 is a diagonal matrix, then all successive Dk are diagonal too. From (22), it is enough to prove that {Dk} converges to sign⁡(Λ), in order to ensure the convergence of the sequence generated by (12) to sign⁡(A).
We can write (22) as n uncoupled scalar iterations to solve g(x)=x2-1=0, given by
(23)dk+1i=7dki+30dki3+11dki51+20dki2+25dki4+2dki6,
where dki=(Dk)i,i and 1≤i≤n. From (22) and (23), it is enough to study the convergence of {dki} to sign⁡(λi), for all 1≤i≤n.
From (23) and since the eigenvalues of A are not pure imaginary, we have that sign⁡(λi)=si=±1. Thus, we attain
(24)dk+1i-1dk+1i+1=-(-1+dki)5(-1+2dki)(1+dki)5(1+2dki).
From the other point of view, since |d0i|=|λi|>0, we obtain
(25)lim⁡k→∞|dk+1i-1dk+1i+1|=0,
and lim⁡k→∞|dki|=1=|sign⁡(λi)|. This shows that {dki} is convergent. Now, one concludes that lim⁡k→∞Dk=sign⁡(Λ).
RecallingDk=V-1XkV, we have
(26)lim⁡k→∞Xk=V(lim⁡k→∞Dk)V-1=Vsign⁡(Λ)V-1=sign⁡(A),
and subsequently the convergence is established. The proof is complete.Theorem 4.
Consider the same conditions of Lemma2 and Theorem 3. Then the proposed method (12) has fifth order of convergence to the sign matrix S.Proof .
Clearly,Xk are rational functions of A and hence, like A, commute with S. On the other hand, we know that S2=I, S-1=S, S2j=I, and S2j+1=S, j≥1. Using the replacement Bk=I+20Xk2+25Xk4+2Xk6 (for the sake of simplicity), we have
(27)Xk+1-S=(7Xk+30Xk3+11Xk5)Bk-1-S=(7Xk+30Xk3+11Xk5-SBk)Bk-1=(-S+7Xk-20SXk2+30Xk3-25SXk4+11Xk5-2SXk6)Bk-1=(-S5+5S4Xk+2Xk-10S3Xk2-10SXk2+10S2Xk3+20Xk3-5SXk4-20SXk4-S5+Xk5+10Xk5-2SXk)Bk-1=((Xk-S)5+2Xk-10SXk2+20Xk3-20SXk4+10Xk5-2SXk)Bk-1=((Xk-S)5-2SXk×(-S+5Xk-10SXk2+10Xk3-5SXk4+Xk5))Bk-1=((Xk-S)5-2SXk(Xk-S)5)Bk-1=(Xk-S)5(I-2SXk)Bk-1.
Now, using any matrix norm from both sides of (27), we attain
(28)∥Xk+1-S∥≤(∥Bk-1∥∥I-2SXk∥)∥Xk-S∥5.
This reveals the fifth order of convergence for the new method (12). The proof is complete.


## 3. Multiplication-Rich and Scaled Variants of the New Method

In general, reduction of a problem in numerical linear algebra to the matrix inversion problem is not an advisable technique. The iterations (5)–(9) as well as (12) require explicit computation of a matrix inverse in each iterative step. Since explicit usage of the inverse matrix is relatively rare in numerical analysis, there is a normal aspiration to approximate the inverse. Such discussions and variants for (12) will be given in the following subsections.


### 3.1. Solve a Matrix Equation Instead of the Matrix Inverse

One of the ways to avoid explicit usage of the inverse matrix is to solve corresponding system of linear matrix equations instead of the matrix inverse. This is done in Algorithm1.Algorithm 1:The new method for computing the matrix sign.
(1)Given a suitable X0∈Cn×n
(2)for k=0,1,… until convergence do
(3)Bk=I+20Xk2+25Xk4+2Xk6
(4)Solve the linear system BkXk+1=7Xk+30Xk3+11Xk5 for Xk+1
(5)end for


### 3.2. Use Approximation of the Matrix Inverse

Another tendency is to give matrix multiplication-rich iterations which retain the convergence rate of the method. For example, the inverseXk-1 in (5) can be replaced by one step of Schulz method for the matrix inverse, which has the form Xk(2I-Xk2). This replacement produces the Newton-Schulz iteration
(29)Xk+1=12Xk(3I-Xk2),X0=A,
which is multiplication-rich and retains the quadratic convergence of Newton’s method. However, it is only locally convergent, with convergence guaranteed for ∥I-A2∥<1 (see [1]). We apply similar idea to (12). For the sake of simplicity, we use the notation Bk=I+20Xk2+25Xk4+2Xk6. Replacing Bk-1 with Bk(2I-Bk2) in (12), we get the following iterative rule for computing the matrix sign function:
(30)Xk+1=Xk(7I+30Xk2+11Xk4)×(I+20Xk2+25Xk4+2Xk6)×(I-40Xk2-450Xk4-1004Xk6-705Xk8-100Xk10-4Xk12).Formula (30) is multiplication-rich with convergence guaranteed for ∥I-A2∥<1.Note that inverse-free algorithms are suitable for the implementation on vector and parallel computers. The iterative scheme (30) includes 9 matrix multiplications, while the complexity of (12) contains 6 matrix multiplications and one matrix inversion to achieve fifth convergence order.


### 3.3. A Hybrid Method

An efficient algorithm for finding the sign by avoiding the computation of matrix inverse is to use (12) or (33) in initial iterations, until Xk2 is close enough to I, and in a subsequent stage apply (30). Such a switching approach was proposed originally in [19]. A similar idea is exploited in Algorithm 2.Algorithm 2:The new (inversion-free) method for computing the matrix sign.
(
1)Given a suitable X0∈Cn×n
(2)use (12) until ∥Xk2-I∥F<ε=1
(3)set X0=Xk
(4)for k=0,1,… until convergence do
(5)Ck=(I-40Xk2-450Xk4-1004Xk6-705Xk8-100Xk10-4Xk12)
(6)Xk+1=Xk(7I+30Xk2+11Xk4)(I+20Xk2+25Xk4+2Xk6)Ck
(7)end for


### 3.4. Scaling Method

For lower order methods, such as (5), the convergence is slow at the beginning. In fact, once the error is sufficiently small (in practice, less than, say, 0.5), successive errors decrease rapidly, each being approximately the square of the previous one. However, initial convergence can be slow if the iteration Xk has a large eigenvalue, that is, in the case ∥Xk∥≫1. Hence, a scaling approach to accelerate the beginning of this phase is necessary and can be done in what follows [7] for the Newton’s method:
(31)X0=A,μk=isthescalingparametercomputedby(32),Xk+1=12(μkXk+μk-1Xk-1),
wherein
(32)μk={∥Xk-1∥∥Xk∥,(norm⁡scaling),ρ(Xk-1)ρ(Xk),(spectralscaling),|det⁡(Xk)|-1/n,(determinantalscaling).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 new iteration (12) is quite fast and reliable due to the discussions in Sections 2 and 3. However, a way is open for speeding up its initial phase of convergence via the concept of scaling.An effective way to enhance the initial speed of convergence is to scale the iterates prior to each iteration, that is,Xk is replaced by μkXk. Such an idea can simply be done in what follows:
(33)X0=A,μk=isthescalingparametercomputedby(32),Xk+1=(7μkXk+30μk3Xk3+11μk5Xk5)×[I+20μk2Xk2+25μk4Xk4+2μk6Xk6]-1,
where lim⁡k→∞μk=1 and lim⁡k→∞Xk=S. Algorithm 3 illustrates an efficient method based on (33) for finding S.Algorithm 3:The new scaling method for computing the matrix sign.
(1)Given a suitable X0∈Cn×n
(2)use (33) until ∥Xk2-I∥F<ε=1
(3)set X0=Xk
(4)for k=0,1,… until convergence do
(5)Ck=(I-40Xk2-450Xk4-1004Xk6-705Xk8-100Xk10-4Xk12)
(6)Xk+1=Xk(7I+30Xk2+11Xk4)(I+20Xk2+25Xk4+2Xk6)Ck
(7)end for


## 3.1. Solve a Matrix Equation Instead of the Matrix Inverse

One of the ways to avoid explicit usage of the inverse matrix is to solve corresponding system of linear matrix equations instead of the matrix inverse. This is done in Algorithm1.Algorithm 1:The new method for computing the matrix sign.
(1)Given a suitable X0∈Cn×n
(2)for k=0,1,… until convergence do
(3)Bk=I+20Xk2+25Xk4+2Xk6
(4)Solve the linear system BkXk+1=7Xk+30Xk3+11Xk5 for Xk+1
(5)end for


## 3.2. Use Approximation of the Matrix Inverse

Another tendency is to give matrix multiplication-rich iterations which retain the convergence rate of the method. For example, the inverseXk-1 in (5) can be replaced by one step of Schulz method for the matrix inverse, which has the form Xk(2I-Xk2). This replacement produces the Newton-Schulz iteration
(29)Xk+1=12Xk(3I-Xk2),X0=A,
which is multiplication-rich and retains the quadratic convergence of Newton’s method. However, it is only locally convergent, with convergence guaranteed for ∥I-A2∥<1 (see [1]). We apply similar idea to (12). For the sake of simplicity, we use the notation Bk=I+20Xk2+25Xk4+2Xk6. Replacing Bk-1 with Bk(2I-Bk2) in (12), we get the following iterative rule for computing the matrix sign function:
(30)Xk+1=Xk(7I+30Xk2+11Xk4)×(I+20Xk2+25Xk4+2Xk6)×(I-40Xk2-450Xk4-1004Xk6-705Xk8-100Xk10-4Xk12).Formula (30) is multiplication-rich with convergence guaranteed for ∥I-A2∥<1.Note that inverse-free algorithms are suitable for the implementation on vector and parallel computers. The iterative scheme (30) includes 9 matrix multiplications, while the complexity of (12) contains 6 matrix multiplications and one matrix inversion to achieve fifth convergence order.


## 3.3. A Hybrid Method

An efficient algorithm for finding the sign by avoiding the computation of matrix inverse is to use (12) or (33) in initial iterations, until Xk2 is close enough to I, and in a subsequent stage apply (30). Such a switching approach was proposed originally in [19]. A similar idea is exploited in Algorithm 2.Algorithm 2:The new (inversion-free) method for computing the matrix sign.
(
1)Given a suitable X0∈Cn×n
(2)use (12) until ∥Xk2-I∥F<ε=1
(3)set X0=Xk
(4)for k=0,1,… until convergence do
(5)Ck=(I-40Xk2-450Xk4-1004Xk6-705Xk8-100Xk10-4Xk12)
(6)Xk+1=Xk(7I+30Xk2+11Xk4)(I+20Xk2+25Xk4+2Xk6)Ck
(7)end for


## 3.4. Scaling Method

For lower order methods, such as (5), the convergence is slow at the beginning. In fact, once the error is sufficiently small (in practice, less than, say, 0.5), successive errors decrease rapidly, each being approximately the square of the previous one. However, initial convergence can be slow if the iteration Xk has a large eigenvalue, that is, in the case ∥Xk∥≫1. Hence, a scaling approach to accelerate the beginning of this phase is necessary and can be done in what follows [7] for the Newton’s method:
(31)X0=A,μk=isthescalingparametercomputedby(32),Xk+1=12(μkXk+μk-1Xk-1),
wherein
(32)μk={∥Xk-1∥∥Xk∥,(norm⁡scaling),ρ(Xk-1)ρ(Xk),(spectralscaling),|det⁡(Xk)|-1/n,(determinantalscaling).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 new iteration (12) is quite fast and reliable due to the discussions in Sections 2 and 3. However, a way is open for speeding up its initial phase of convergence via the concept of scaling.An effective way to enhance the initial speed of convergence is to scale the iterates prior to each iteration, that is,Xk is replaced by μkXk. Such an idea can simply be done in what follows:
(33)X0=A,μk=isthescalingparametercomputedby(32),Xk+1=(7μkXk+30μk3Xk3+11μk5Xk5)×[I+20μk2Xk2+25μk4Xk4+2μk6Xk6]-1,
where lim⁡k→∞μk=1 and lim⁡k→∞Xk=S. Algorithm 3 illustrates an efficient method based on (33) for finding S.Algorithm 3:The new scaling method for computing the matrix sign.
(1)Given a suitable X0∈Cn×n
(2)use (33) until ∥Xk2-I∥F<ε=1
(3)set X0=Xk
(4)for k=0,1,… until convergence do
(5)Ck=(I-40Xk2-450Xk4-1004Xk6-705Xk8-100Xk10-4Xk12)
(6)Xk+1=Xk(7I+30Xk2+11Xk4)(I+20Xk2+25Xk4+2Xk6)Ck
(7)end for


## 4. Experimental Results

There are basically two general ways to terminate a matrix iteration for findingS, that is, to stop when Xk has relative error below a tolerance or to stop when some suitable residual (such as Xk2-I) is below the tolerance. The relevant aim will in general be problem dependent. However, the stop termination is really important in matrix methods.The considered stopping termination in this section would be the safe strategy introduced in [1] as follows:
(34)∥Xk2-I∥*∥Xk∥*2≤ϵ.For comparisons, we have used the matrix globally convergent methods (5), (6), and (12) using Mathematica 8 built-in precision, [21]. We used Mathematica function Inverse[] to compute the required matrix inverse. Implementation details of the function Inverse[] are based on efficient row reduction (Gaussian elimination) based on numerical approximation.Example 5.
In this test, we examine the behavior of different iterative methods for finding the matrix sign function of 10 randomly generated complex square70×70 matrices as follows: n  =  70;  number  =  10;  SeedRandom[12345];Table[A[l]  = RandomComplex[{-5-I, 5+I},  {n, n}];,  {l, number}].The results of comparisons in terms of the number of iterations have been reported in Table2. We remark that whatever the eigenvalues of a matrix are closer to the imaginary axis, the speed of convergence for different methods becomes slower and more risky to face with singular matrices Xk, whose inverses could not be computed.Table 2
Results of comparisons for Example5.
Matrices
#
1
#
2
#
3
#
4
#
5
#
6
#
7
#
8
#
9
#
10
Number of iterations for (5)
15
16
14
14
15
14
13
13
13
13
Number of iterations for (6)
7
7
6
6
7
6
6
6
6
6
Number of iterations for (12)
6
7
6
6
6
6
5
6
6
5In this example, we have used the stopping criterion (34) with ϵ=10-10, matrix norm one, and X0=A as the initial matrix.Example 6.
In order to confirm the acceleration via scaling, we have reran Example5 with the norm scaling. The results are summarized in Table 3. The numerics reverify the effectiveness of the new method (12) in finding the matrix sign.Table 3
Results of comparisons for Example6 with norm scaling.
Matrices
#
1
#
2
#
3
#
4
#
5
#
6
#
7
#
8
#
9
#
10
Number of iterations for the scaled (5)
12
11
11
12
13
12
11
10
14
10
Number of iterations for the scaled (6)
5
6
5
5
5
5
5
5
5
5
Number of iterations for the scaled (12)
5
5
5
5
5
5
4
5
5
4


## 5. Summary

Interest in the sign function grew steadily starting from the 1970s and 1980s, initially among engineers and later among numerical analysts. Following such a trend, in this study we proposed a fifth order new iterative method for finding the matrix signS.Applying the basins of attraction in the complex plane, based on the results from [20], we concluded that the introduced method (12) has global convergence. We then theoretically found that it is asymptotically stable. Some numerical experiments have been studied in order to show the faster convergence on the basis of a smaller number of iterations.Several modifications of the introduced method have been established, such as a new inversion-free method, a composite method, and a scaled version.



---

*Source: 105301-2014-07-17.xml*