Computing the Singular Value Decomposition of 3x3 matrices with minimal branching and elementary floating point operations

Abstract

A numerical method for the computation of the Singular Value Decomposition of 3x3 matrices is presented.The proposed methodology robustly handles rank-deficient matrices and guarantees orthonormality of the computed rotational factors. The algorithm is tailored to the characteristics of SIMD or vector processors. In particular, it does not require any explicit branching beyond simple conditional assignments (as in the C++ ternary operator "?:", or the SSE4.1 instruction VBLENDPS), enabling trivial data-level parallelism for any number of operations. Furthermore, no trigonometric or other expensive operations are required; the only floating point operations utilized are addition, multiplication, and an inexact (yet fast) reciprocal square root which is broadly available on current SIMD/vector architectures. The performance observed approaches the limit of making the 3x3 SVD a memory-bound (as opposed to CPU-bound) operation on current SMP platforms.