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

Authors: F. Khaksar Haghani
Journal: Advances in Numerical Analysis
 (2014)

---

## Abstract

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.


---

## Body

## 1. Introduction

The polar decomposition ofA∈Cm×n factors A as the product
(1)A=UH,U*U=Ir,rank⁡(U)=r=rank⁡(A),
where U is unitary and H of order n is Hermitian positive semidefinite. The Hermitian factor H is always unique and can be expressed as
(2)H=(A*A)1/2,
and the unitary factor U∈Cm×n is unique if A 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 m≥n.This matrix decomposition has many applications in various fields. To give an example, general3×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 ofA 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]:
(3)U=2π∫0∞(t2I+A*A)-1dt.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 A2 by A*A 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:
(4)Uk+1=12(Uk+Uk-*),
while its following alternative for general rectangular cases was considered in [6] as
(5)Uk+1=12(Uk+Uk†*),
wherein U† stands for the Moore-Penrose generalized inverse. Note that, throughout this work, Uk-* stands for (Uk-1)*. 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 G∈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:
(6)Uk+1=[Uk(3I+Uk*Uk)][I+3Uk*Uk]-1.An initial matrixU0 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
(7)U0=1αA,
whereas α>0 is an estimate of ∥A∥2.The remaining sections of this paper are organized as follows. In Section2, we derive an iteration function for polar decomposition. The scheme is convergent to the unitary polar factor U, 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 ofA is based on applying a zero-finder to a particular map. That is, solving the following nonlinear (matrix) equation
(8)F(U)∶=U*U-I=0,
in which I 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:(9)xk+1=xk-20-9L(xk)20-19L(xk)f(xk)f′(xk),
with L(xk)=f′′(xk)f(xk)/f′(xk)2.Theorem 1.
Letα∈D be a simple zero of a sufficiently differentiable function f:D⊆C→C for an open interval D, which contains x0 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 x2-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.Attraction basins for (5) (a) and (11) (b), for the polynomial g(x)=x2-1.
(a)
(b)In terms of the fractal theory, it is necessary to find the global basins of attraction for a zeroz*:
(10)S(z*)∶={z∈C:Lk(z)⟶z*,ask⟶∞},
where Lk(x)=L(L(⋯(L(x)))) is the k-fold composition L∘⋯∘L of the iteration function L. Here, using (9), we have (in the reciprocal form)
(11)xk+1=L(xk)=38xk+42xk39+60xk2+11xk4,k=0,1,….To check the global convergence for the quadratic polynomialg(z)=z2-1 with the zeros ±1, (11) gives the following formula:
(12)B(z)=38z+42z39+60z2+11z4,
and we find
(13)B(z)+1B(z)-1=λ(U+1U-1)3,
wherein λ=-(9+11U)/(-9+11U).Let∂S denote the boundary of the set S. One of the basic notions in the fractal theory connected to iterative processes and convergence of an iterative function f is Julia set for the proposed operator L. Thus, when k→∞, we obtain(1)
if|(z+1)/(z-1)|<1, then |(Bk(z)+1)/(Bk(z)-1)|→0, and Bk(z)→-1;(2)
if|(z+1)/(z-1)|>1, then |(Bk(z)+1)/(Bk(z)-1)|→0, and Bk(z)→+1.Furthermore, we can conclude that the basins of attractionS(-1) and S(1) in the case of operator B are the half-planes on either side in relation to the line z=0 (the imaginary axis). Since ±1 are attractive fixed points of B, so the Julia set J(B) is the boundary of the basins of attraction S(-1) and S(1); that is,
(14)J(B)=∂S(-1)=∂S(1)={γi:γ∈R}.Actually, the Julia setJ(B) is just the line z=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:
(15)Uk+1=Uk[38I+42Yk][9I+60Yk+11Zk]-1,
where Yk=Uk*Uk, Zk=YkYk, and U0 is given by (7).Theorem 2.
Assume thatA∈Cm×n is an arbitrary matrix. Then, the matrix iterates {Uk}k=0k=∞ of (15) converge to U.Proof.
LetA have the following SVD form:
(16)A=PΣQ*,
where
(17)Σ=(Σr000),
and r=rank⁡(A). Zeros in Σ may not exist. We define the following sequence of matrices:
(18)Dk=P*UkQ.
On the other hand, using (15), one may obtain
(19)D0=Σ,Dk+1=[38Dk+42Dk3][9I+60Dk2+11Dk4]-1.
Since D0 is diagonal with positive diagonal and zero elements, it follows by induction that the sequence {Dk}k=0∞ is defined by
(20)Dk=(diag⁡(di(k))000),di(k)>0.
Accordingly, (19) represents r uncoupled scalar iterations as follows:
(21)di(0)=σi,1≤i≤r,di(k+1)=[38di(k)+42di(k)3][9+60di(k)2+11di(k)4]-1.
Simple manipulations yield the relation
(22)di(k+1)-1di(k+1)+1=-9-38di(k)+60di(k)2-42di(k)3+11di(k)49+38di(k)+60di(k)2+42di(k)3+11di(k)4.
Since σi is positive, (22) holds for each i. It follows that
(23)lim⁡k→∞|di(k+1)-1di(k+1)+1|⟶0.
That is to say,
(24)Dk⟶(Ir000).
Therefore, Uk→U and subsequently H=U*A. The proof is complete.Theorem 3.
LetA∈Cm×n be an arbitrary matrix. Then, new method (15) has third order to find the unitary polar factor of A.Proof.
Proposed scheme (15) transforms the singular values of Uk according to
(25)σi(k+1)=[38σi(k)+42σi(k)3][9+60σi(k)2+11σi(k)4]-1
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 k≥1 as follows:
(26)σi(k+1)-1σi(k+1)+1=-(-1+σi(k))3(-9+11σi(k))(1+σi(k))3(9+11σi(k)).
Now, we attain
(27)|σi(k+1)-1σi(k+1)+1|≤(|-9+11σi(k)9+11σi(k)|)|σi(k)-1σi(k)+1|3.
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 matrixUk before each cycle. An important scaling approach was derived in [13] in Frobenius norm as comes next:
(28)θk=(∥Uk†∥F∥Uk∥F)1/2.So, the new scheme can be expressed in the following accelerated form:(29)Computeθkby(28),k≥0,Mk=9I+60θk2Yk+11θk4Zk,Uk+1=θkUk[38I+42θk2Yk]Mk-1.


## 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
(30)Computeθkby(28),k≥0,Uk+1=12[θkUk+θk-1Uk†*],
have been tested and compared. The stopping termination in this work is
(31)Rk+1=∥Uk+1-Uk∥∞∥Uk∥∞≤ϵ,
wherein ϵ is the tolerance.Example 1.
In this experiment, we compare the behavior of different methods. We used the following six complex randomly generated rectangular310×300 matrices: m = 310; n = 300; number = 6; SeedRandom[345];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 ϵ=10-10 with U0=A. It could easily be observed that there is a clear reduction in the number of iterations using PMP.Table 1
Results of comparisons for Example1 in terms of number of iterations.
Methods
NMP
HMP
ANMP
PMP
#1
11
8
9
6
#2
11
8
9
6
#3
11
8
9
6
#4
11
8
9
6
#5
11
8
9
6
#6
12
8
9
7Table 2
Results of comparisons for Example1 in terms of elapsed time (s).
Methods
NMP
HMP
ANMP
PMP
#1
6.0312500
3.0156250
4.7656250
2.7812500
#2
5.2812500
3.0781250
4.7500000
2.8437500
#3
5.2500000
2.9375000
4.7656250
2.7343750
#4
5.2656250
3.0312500
4.7500000
2.8125000
#5
5.2656250
2.9375000
4.7031250
2.7500000
#6
5.7343750
3.0312500
4.7656250
3.3125000These theoretical results of Section2 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 numberz=reiθ admits a unique polar representation with r∈R+, φ(-π,+π], every matrix A can be decomposed into a product of the unitary polar factor U and a positive semidefinite matrix H. 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 factorU. 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.



---

*Source: 576325-2014-09-30.xml*