A parallel algorithm for computing the polar decomposition

The polar decomposition A = UH of a rectangular matrix A, where U is unitary and H is Hermitian positive semidefinite, is an important tool in various applications, including aerospace computations, factor analysis and signal processing. We consider a pth order iteration for computing U that involves p independent matrix inversions per step and which is hence very amenable to parallel computation. We show that scaling the iterates speeds convergence of the iteration but makes the iteration only conditionally stable, with the backward error typically κappa_2(A) times bigger than the unit roundoff. In our implementation of the iteration on the Kendall Square Research KSR1 virtaul shared memory MIMD computer we take p to be the number of processors (p ≤ 16 in our experiments). Our code is found to be significantly faster than two existing techniques for computing the polar decomposition: one a Newton iteration, the other based on the singular value decomposition