It is attempted to present an iteration method for finding polar decomposition. The approach is categorized in the scope of Newton-type methods. Error analysis and rate of convergence are studied. Some illustrations are also given to disclose the numerical behavior of the proposed method.

1. Introduction

The polar decomposition of factors as the product 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 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 linear or 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 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 . Similar notations are used throughout.

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 , 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 must be employed in matrix fixed-point type methods such as (4)–(6). An initial approximation for the unitary factor of any matrices can be expressed as whereas is an estimate of .

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.

2. 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 root-finding 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 for an open interval , which contains 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 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.

In terms of the fractal theory, it is necessary to find the global basins of attraction for a zero : where is the -fold composition of the iteration function . Here, using (9), we have (in the reciprocal form)

To check the global convergence for the quadratic polynomial with the zeros , (11) gives the following formula: and we find wherein .

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(1)if , then , and ;(2)if , then , and .

Furthermore, we can conclude that the basins of attraction and in the case of operator are the half-planes on either side in relation to the line (the imaginary axis). Since are attractive fixed points of , so the Julia set is the boundary of the basins of attraction and ; that is,

Actually, the Julia set is just the line 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 is given by (7).

Theorem 2. Assume that is an arbitrary matrix. Then, the matrix iterates of (15) converge to .

Proof. Let have the following SVD form: where and . Zeros in may not exist. We define the following sequence of matrices: On the other hand, using (15), one may obtain Since is diagonal with positive diagonal and zero elements, it follows by induction that the sequence is defined by Accordingly, (19) represents uncoupled scalar iterations as follows: Simple manipulations yield the relation Since is positive, (22) holds for each . It follows that That is to say, Therefore, and subsequently . The proof is complete.

Theorem 3. Let 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 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:

3. 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 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 matrices: m = 310; n = 300; number = 6; SeedRandom;Table[A[l] = RandomComplex[{-10-10 I,10+10 I}, {m, n}];, {l, number}];The results of comparison are carried out in Tables 1 and 2 applying the tolerance with . 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.

4. 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 , , 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.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.


The author sincerely thanks Islamic Azad University, Shahrekord Branch, Shahrekord, Iran, for financial support.