The Singular Value Decomposition: Anatomy of Optimizing an Algorithm for Extreme Scale

Abstract

The computation of the singular value decomposition, or SVD, has a long history with many improvements over the years, both in its implementations and algorithmically. Here, we survey the evolution of SVD algorithms for dense matrices, discussing the motivation and performance impacts of changes. There are two main branches of dense SVD methods: bidiagonalization and Jacobi. Bidiagonalization methods started with the implementation by Golub and Reinsch in Algol60, which was subsequently ported to Fortran in the EIS-PACK library, and was later more efficiently implemented in the LINPACK library, targeting contemporary vector machines. To address cache-based memory hierarchies, the SVD algorithm was reformulated to use Level 3 BLAS in the LAPACK library. To address new architectures, ScaLAPACK was introduced to take advantage of distributed computing, and MAGMA was developed for accelerators such as GPUs. Algorithmically, the divide and conquer and MRRR algorithms were developed to reduce the number of operations. Still, these methods remained memory bound, so two-stage algorithms were developed to reduce memory operations and increase the computational intensity, with efficient implementations in PLASMA, DPLASMA, and MAGMA. Jacobi methods started with the two-sided method of Kogbetliantz and the one-sided method of Hestenes. They have likewise had many developments, including parallel and block versions and preconditioning to improve convergence. In this paper, we investigate the impact of these changes by testing various historical and current implementations on a common, modern multicore machine and a distributed computing platform. We show that algorithmic and implementation improvements have increased the speed of the SVD by several orders of magnitude, while using up to 40 times less energy.

Backward error and condition number of a generalized Sylvester equation, with application to the stochastic Galerkin method

“…From a computational point of view evaluating (5.1) is not expensive because B 1 is a diagonal matrix, and m n, that is Y is a tall, skinny matrix, and therefore σ m (Y ) can be computed very efficiently [17,Sec 5.4]. Similarly the expression for µ B in (3.20) can also be simplified, which gives…”

Reducing the amount of out‐of‐core data access for GPU‐accelerated randomized SVD

“…EISPACK was designed to run on a single‐core CPU and was replaced by LINPACK, which first implemented the SVD algorithm with basic linear algebra subprogram (BLAS) interface. The performance of LINPACK was limited by the BLAS1 implementation and benefited little from multicore architectures . LAPACK redesigned the SVD algorithm to use BLAS3 routines wherever possible to improve the performance on the multicore CPUs.…”

“…This low efficiency is due to the computation of tall‐skinny GEMM, which is closer to GEMV (BLAS2 routines) than GEMM (BLAS3 routines). The BLAS2 routines are less efficient than the BLAS3 routines due to vector accesses that degrade cache hit rate on both the multicore CPU and GPU; BLAS3 routines are 20 to 40 times more efficient than BLAS2 routines . Regarding the in‐core performance of tall‐skinny GEMM, Chen et al achieved 1.1 to 3.0× speedups over cuBLAS for tall‐skinny matrices with up to 16 columns.…”

“…Both underlying methods can be implemented using general matrix‐matrix multiplication (GEMM) routines, whereas general matrix‐vector multiplication (GEMV) is required for deterministic SVD computation. GEMM is more suitable for modern parallel computers, where it achieves 20 to 40 times higher flop/s than GEMV …”

“…First, large matrix data may not fit into the GPU memory due to its limited capacity. Second, communication across distinct memory hierarchies or networks often constitutes a performance bottleneck due to the increasing gap between arithmetic and communication performance …”