Numerical behaviour of Higham's scaled method for polar decomposition (Q1397923)

As is well known, a nonsingular complex matrix \(A\) has a unique polar decomposition \( A=U_AH_A, \) where \(U_A\) is a unitary matrix and \(H_A\) is Hermitian positive-definite. The scaled method of \textit{N. J. Higham} [SIAM J. Sci. Stat. Comput. 7, 1160--1174 (1986; Zbl 0607.65014)] is an iterative method for computing the unitary factor \(U_A\). It consists in computing a sequence of matrices \((X_k)_{k=0}^{\infty}\), by \[ X_{k+1}= \frac{1}{2} \biggl(\gamma_kX_k+ \frac{1}{\gamma_k}X_k^{-H}\biggr),\quad X_0=A, \;\gamma_k>0. \] In exact arithmetic, this iteration always converges to \(U_A\) if \(\gamma=1\). Some other choices for \(\gamma\) accelerate the convergence. However, if \(cond(A)\) is large, then the numerical computation of the inverse of \(X_k\) in the initial stages could cause a loss of accuracy. Extensive numerical tests have shown that the loss of accuracy is in reality rather small, except in extremely ill-conditioned situations. The authors present a detailed error analysis to explain this phenomenon and at the same time improve the algorithm by studying the influence of the choice of the parameter \(\gamma\) on the convergence behaviour. In essence, the authors formulate assumptions on the quality of the numerical matrix inversion used in the iteration, on the scaling, the machine precision, the dimension of the matrix and its condition number. The conclusion is that the computed \(X_k\) is an acceptable approximation of \(U_A\) after a few steps if the conditions are satisfied. The error analysis is tedious but confirmed by extensive numerical tests.