A Backward Stable Algorithm for Computing the CS Decomposition via the Polar Decomposition

TL;DR

This paper presents a backward stable, parallelizable algorithm for computing the CS decomposition of matrices using polar decompositions and eigendecomposition, demonstrating high numerical stability through numerical experiments.

Contribution

The paper introduces a novel backward stable algorithm for CS decomposition that leverages parallelizable polar and eigendecompositions, enhancing stability and efficiency.

Findings

Algorithm is backward stable when polar and eigendecompositions are stable.

The method is highly parallelizable due to the nature of the subcomputations.

Numerical experiments confirm excellent stability and performance.

Abstract

We introduce a backward stable algorithm for computing the CS decomposition of a partitioned $2n \times n$ matrix with orthonormal columns, or a rank-deficient partial isometry. The algorithm computes two $n \times n$ polar decompositions (which can be carried out in parallel) followed by an eigendecomposition of a judiciously crafted $n \times n$ Hermitian matrix. We prove that the algorithm is backward stable whenever the aforementioned decompositions are computed in a backward stable way. Since the polar decomposition and the symmetric eigendecomposition are highly amenable to parallelization, the algorithm inherits this feature. We illustrate this fact by invoking recently developed algorithms for the polar decomposition and symmetric eigendecomposition that leverage Zolotarev's best rational approximations of the sign function. Numerical examples demonstrate that the resulting algorithm for computing the CS decomposition enjoys excellent numerical stability.