We introduce a backward stable algorithm for computing the CS decomposition
of a partitioned 2nΓn matrix with orthonormal columns, or a
rank-deficient partial isometry. The algorithm computes two nΓn polar
decompositions (which can be carried out in parallel) followed by an
eigendecomposition of a judiciously crafted nΓ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