Discrete Dynamics in Nature and Society

Volume 2015 (2015), Article ID 649423, 11 pages

http://dx.doi.org/10.1155/2015/649423

## A Novel Iterative Method for Polar Decomposition and Matrix Sign Function

^{1}Department of Mathematics, Islamic Azad University, Zahedan Branch, Zahedan, Iran^{2}Faculty of Sciences and Mathematics, University of Niš, Visegradska 33, 18000 Niš, Serbia^{3}Faculty of Computer Science, Goce Delčev University, Goce Delčev 89, 2000 Štip, Macedonia

Received 3 April 2015; Accepted 21 June 2015

Academic Editor: Gian I. Bischi

Copyright © 2015 F. Soleymani et al. 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.

#### Abstract

We define and investigate a globally convergent iterative method possessing sixth order of convergence which is intended to calculate the polar decomposition and the matrix sign function. Some analysis of stability and computational complexity are brought forward. The behaviors of the proposed algorithms are illustrated by numerical experiments.

#### 1. Introduction

Recently, a large number of papers intended for calculating simple and multiple zeros of the nonlinear equation have been published. There are several attempts to overcome the difficulties in designing new iterative methods for solving or in the application of such iterative methods in computing matrix functions (see [1, 2]). Our current work discusses an application of a high order iterative method for solving nonlinear scalar equations in calculating the polar decomposition and the matrix sign function.

We now restate some preliminaries about the polar decomposition and the matrix sign.

Let denote the linear space of all complex matrices. The polar decomposition of a complex matrix is defined aswherein is a Hermitian positive semidefinite matrix of the order and is a subunitary matrix (partial isometry). Here, the inequality is assumed. A matrix is subunitary if for any , where and denote the linear space spanned by columns of the matrix (range of ) and the null space of , respectively. Note that if , then , and is an orthonormal Stiefel matrix.

It is assumed that the singular value decomposition (SVD) of has the following form:wherein and are unitary matrices andThe Hermitian factor is always unique and it can be expressed as , and the unitary factor is unique if is nonsingular. Note that, here, the exponent denotes the principal square root: the one whose eigenvalues lie in the right half-plane.

The polar decomposition may be easily constructed from an SVD of the given matrix . However, the SVD is a substantial calculation that displays much more about the structure of than does the polar decomposition. Constructing the polar decomposition from the SVD destroys this extra information and wastes the arithmetic work used to compute it. It is intuitively more appealing to use the polar decomposition as a preliminary step in the computation of the SVD.

On the other hand, a related issue to the polar decomposition is the matrix sign decomposition, which is defined for a matrix having no pure imaginary eigenvalues. The most concise definition of the matrix sign decomposition is given in [3]Here, is the matrix sign function, introduced by Roberts [4]. We herewith remark that, in this work, whenever we write about the computation of a matrix sign function, we mean a square matrix with no eigenvalues on the imaginary axis.

Now we briefly review some of the most important iterative methods for computing the matrix polar decomposition.

Among many iterations available for computing the polar decomposition, the most practically useful is the scaled Newton iteration [5] as well as the recently proposed dynamically weighted Halley iteration [6]. The method of Newton introduced for computing the polar decomposition (via unitary polar factor) is defined in [5] by the iterative schemefor the square nonsingular cases and by the following alternative for general rectangular cases [7]:Note that stands for the Moore-Penrose generalized inverse, , and .

The cubically convergent method of Halleyhas been developed in [8] for computing the unitary polar factor.

The particular formula (7) is also applicable to singular or rectangular matrices. This scheme has further been developed adaptively in [6]. In fact the dynamically weighted Halley method (DWH) for computing unitary polar factor has been introduced as follows:The parameters , , and are dynamically chosen to accelerate the convergence. They are computed bywhereIn (9), is a lower bound for the smallest singular value of . Fortunately, once is obtained, then effective and sharp bounds can be attained at no cost from the following recurrence:

An initial matrix must be employed in a matrix fixed-point type method so as to arrive at the convergence phase. Such a sharp initial approximation for the unitary factor can be expressed aswhereas is an estimate of (a safe choice of which is ).

The rest of this paper unfolds the contents in what follows. Section 2 derives a new iterative scheme for solving nonlinear equations. Additionally, by applying the illustration via basins of attraction, we find a scheme with a global convergence. The global convergence of the constructed solver is verified analytically. In Section 3, we generalize the proposed nonlinear equation solver into the iterative method for finding the polar decomposition. Furthermore, we prove that the new scheme is convergent. The computational complexity is discussed in Section 4. In Section 5, we define an extension of the proposed method in numerical computation of the matrix sign function and show its asymptotic stability. Section 6 is devoted to the application of the contributed methods in solving two numerical examples (one in double precision arithmetic and the other in a high precision computing environment). Finally, Section 7 draws a conclusion of this paper.

#### 2. Chebyshev-Halley Type Iteration and Its Extension

Many of the iterative methods for computing matrix functions can be deduced by applying a nonlinear equation solver to a special mapping. For example, applying Newton’s method to the mappingin which is the identity matrix of the appropriate size, could result in iterates (5) (note that the equivalent form of (13) for the matrix sign is ). This reveals a close relation between matrix functions and iterative root-finding methods [3].

Gutiérrez and Hernández in [9] developed a Chebyshev-Halley type scheme in Banach space for finding simple zeros of , which can be written in what follows:wherein , , and the convergence order is cubic.

If we decide to apply (14) for solving (13), we will obtain a family of at least cubically convergent schemes in finding the polar decomposition. But an important barrier occurred and it would be the nonglobal convergence of some of the schemes.

On the other hand, Iannazzo in [10] showed that the matrix convergence is governed by the scalar convergence* for the so-called pure matrix iterations*. We remark that this is true for the matrix sign function when the scalars are the eigenvalues, while for polar decomposition it is true when scalars are singular values. This is also well illustrated in the textbook [3].

Hence, we should find a member from (14), so that the derived method is new and also possesses the global convergence. Toward this goal, we employ the theory of basins of attraction for (14) so as to solve the quadratic polynomial in a square of the complex plane whereas the maximal number of iterations is fixed to 100 and the stopping criterion is .

This is done in Figures 1–3 for different values of . In Figures 1(a) and 2(a), the convergence in the process of solving is local, and therefore convergence of the matrix iterations could happen only for very sharp initial matrices. Moreover, the method in Figure 3(b) behaves chaotically and also includes divergence points. The divergence points are included in the black areas and, accordingly, should be outside of interest.