A Third-Order Newton-Type Method for Finding Polar Decomposition

and the unitary factorU ∈ Cm×n is unique ifA is nonsingular [1]. The exponent 1/2 denotes the principal square root, that is, the one whose eigenvalues lie in the right half-plane. Here, we assume thatm ≥ n. This matrix decomposition has many applications in various fields. To give an example, general 3 × 3 linear or 4 × 4 homogenous matrices can be formed by composing primitive matrices for translation, rotation, scale, and so on. Current 3D computer graphics systems manipulate and interpolate parametric forms of these primitives to generate scenes and motion [2]. Hence, decomposing a composite matrix is necessary. This paper follows one of such ways, known as the polar decomposition (1). Practical interest in the polar decomposition stems mainly from the fact that the unitary polar factor of A is the nearest unitary matrix to A in any unitarily invariant norm [3]. Apart from (1), the polar decomposition can be defined by the following integral formula [4]:


Introduction
The polar decomposition of  ∈ C × factors  as the product  = ,  *  =   , rank () =  = rank () , where  is unitary and  of order  is Hermitian positive semidefinite.The Hermitian factor  is always unique and can be expressed as and the unitary factor  ∈ C × is unique if  is nonsingular [1].The exponent 1/2 denotes the principal square root, that is, the one whose eigenvalues lie in the right half-plane.Here, we assume that  ≥ .This matrix decomposition has many applications in various fields.To give an example, general 3 × 3 linear or 4 × 4 homogenous matrices can be formed by composing primitive matrices for translation, rotation, scale, and so on.Current 3D computer graphics systems manipulate and interpolate parametric forms of these primitives to generate scenes and motion [2].Hence, decomposing a composite matrix is necessary.This paper follows one of such ways, known as the polar decomposition (1).
Practical interest in the polar decomposition stems mainly from the fact that the unitary polar factor of  is the nearest unitary matrix to  in any unitarily invariant norm [3].
Apart from (1), the polar decomposition can be defined by the following integral formula [4]: Formula (3) illustrates a guiding principle that any property or iteration involving the matrix sign function can be converted into one for the polar decomposition using the replacement  2 by  *  and vice versa.
Here, we concentrate on the iterative expressions for finding (1), since integral representation (3) has some complicated structure and requires complex analysis.
Newton's method for square nonsingular cases introduced in [5] is as follows: while its following alternative for general rectangular cases was considered in [6] as wherein  † stands for the Moore-Penrose generalized inverse.Note that, throughout this work,  − *  stands for ( −1  ) * .Similar notations are used throughout.

Advances in Numerical Analysis
Authors in [7] derived important results for (4).They discussed that although Newton's method for the polar decomposition immediately destroys the underlying group structure, when  ∈ U, it forces equality between the adjoint and the conjugate transpose of each iterate.This implies that the Newton iterates approach the group at the same rate that they approach unitarity.
The cubically convergent method of Halley has been developed for polar decomposition in [8] as follows: An initial matrix  0 must be employed in matrix fixedpoint type methods such as (4)-( 6).An initial approximation for the unitary factor of any matrices can be expressed as whereas  > 0 is an estimate of ‖‖ 2 .
The remaining sections of this paper are organized as follows.In Section 2, we derive an iteration function for polar decomposition.The scheme is convergent to the unitary polar factor , and the rate of convergence is three since the proposed formulation transforms the singular values of the matrices produced per cycle with a cubical rate to one.Some illustrations are also provided to support the theoretical aspects of the paper in Section 3. Finally, conclusions are drawn in Section 4.

A Third-Order Method
The procedure of constructing a new iterative method for the unitary factor of  is based on applying a zero-finder to a particular map.That is, solving the following nonlinear (matrix) equation in which  is the identity matrix, by an appropriate rootfinding method could yield new iteration functions (see, e.g., [9,10]).Therefore, we first introduce the following iterative expression for finding the simple zeros of nonlinear equations: with Theorem 1.Let  ∈  be a simple zero of a sufficiently differentiable function  :  ⊆ C → C for an open interval , which contains  0 as an initial approximation of .Then, iterative expression (9) has third order of convergence.
Proof.The proof is similar to the proofs given in [11].So, it is skipped over.
Drawing the attraction basins of ( 9) for finding the solution of the polynomial equation  2 −1 = 0 in the complex plane reveals that the application of ( 9) for finding matrix sign function and consequently the polar decomposition has global convergence (see Figure 1).However, it is necessary to show this global behavior analytically.
Let  denote the boundary of the set .One of the basic notions in the fractal theory connected to iterative processes and convergence of an iterative function  is Julia set for the proposed operator .Thus, when  → ∞, we obtain Furthermore, we can conclude that the basins of attraction (−1) and (1) in the case of operator  are the halfplanes on either side in relation to the line  = 0 (the imaginary axis).Since ±1 are attractive fixed points of , so the Julia set () is the boundary of the basins of attraction (−1) and (1); that is, Actually, the Julia set () is just the line  = 0 for (12), and thus the new third-order method ( 11) is globally convergent.
By taking into account this global behavior, we extend (11) as follows: where   =  *    ,   =     , and  0 is given by (7).
Proof.Let  have the following SVD form: where and  = rank().Zeros in Σ may not exist.We define the following sequence of matrices: On the other hand, using (15), one may obtain Since  0 is diagonal with positive diagonal and zero elements, it follows by induction that the sequence {  } ∞ =0 is defined by Accordingly, (19) represents  uncoupled scalar iterations as follows: Simple manipulations yield the relation That is to say, Therefore,   →  and subsequently  =  * .The proof is complete.
Theorem 3. Let  ∈ C × be an arbitrary matrix.Then, new method (15) has third order to find the unitary polar factor of .
Proof.Proposed scheme (15) transforms the singular values of   according to and leaves the singular vectors invariant.From (25), it is enough to show that convergence of the singular values to unity has third order for  ≥ 1 as follows: Now, we attain This reveals the third order of convergence for new method (15).The proof is ended.
The proposed method is not a member of Padé family of iterations given in [12], with global convergence.So, it is interesting from both theoretical and computational points of view.
The speed of convergence can be slow at the beginning of the process; so, it is necessary to scale the matrix   before each cycle.An important scaling approach was derived in [13] in Frobenius norm as comes next: So, the new scheme can be expressed in the following accelerated form: Compute   by (28) ,  ≥ 0, (29)

Numerical Examples
We have tested contributed method (15) denoted by PMP, using the programming package Mathematica 8 in double precision [14].Apart from this scheme, several iterative methods, such as (5) denoted by NMP and (6) denoted by HMP, and accelerated Newton's method given by Compute   by (28) ,  ≥ 0, have been tested and compared.The stopping termination in this work is wherein  is the tolerance.
Example 1.In this experiment, we compare the behavior of different methods.We used the following six complex randomly generated rectangular 310 × 300 matrices: m = 310; n = 300; number = 6; SeedRandom [345]; The results of comparison are carried out in Tables 1 and 2 applying the tolerance  = 10 −10 with  0 = .It could easily be observed that there is a clear reduction in the number of iterations using PMP.
These theoretical results of Section 2 have been confirmed by numerical examples here.So, we demonstrate the convergence behavior of proposed iteration (15).Note that the superiority of PMP in contrast to HMP is understandable from the fact that PMP has larger attraction basins, and subsequently it could cluster the singular values to unity much faster than HMP.

Discussion
Matrix decomposition is well established as an important part of computer graphics.Just as every nonzero complex number  =   admits a unique polar representation with  ∈ R + , (−, +], every matrix  can be decomposed into a product of the unitary polar factor  and a positive semidefinite matrix .The polar decomposition is of interest in many applications, for example, whenever it is required to orthogonalize a matrix.In this paper, we have developed a new method for finding the unitary polar factor .It has been shown that the convergence is global and its rate is three.Scaling form of our proposed method has also been given.From numerical results, we observe that accuracy in successive approximations increases, showing stable nature of the method.Also, like the existing methods, the presented method shows consistent convergence behavior.Further improvement of convergence rate can be considered for future studies.

Table 1 :
Results of comparisons for Example 1 in terms of number of iterations.

Table 2 :
Results of comparisons for Example 1 in terms of elapsed time (s).