Two-sided Jacobi algorithm for Singular Value Decomposition (SVD)

Introduction

The SVD of a (real square, for simplicity) matrix A is a representation of the matrix in the form

A = U D V^T ,

where matrix D is diagonal with non-negative elements and matrices U and V are orghogonal. The diagonal elements of matrix D can always be chosen non-negative by multiplying the relevant columns of matrix U with (-1).

SVD can be used to solve a number of problems in linear algebra.

Problem

Implement the two-sided Jacobi SVD algorithm.

Algorithm (as described in the "Eigenvalues" chapter of the book)

In the cyclic Jacobi eigenvalue algorithm for symmetric matrices we applied the elementary Jacobi transformation

A → J^T A J

where J ≐ J(θ,p,q) is the Jacobi matrix (where the angle θ is chosen to eliminate the off-diagonal elements A_pq=A_qp) to all upper off-diagonal matrix elements in cyclic sweeps until the matrix becomes diagonal.

The two-sided Jacobi SVD algorithm for general real square matrices is the same as the Jacobi eigenvalue algorithm, except that the elementary transformation is slightly different,

A → J^T G^T A J ,

where

• G ≐ G(θ,p,q) is the Givens (Jacobi) matrix where the angle is chosen such that the matrix A' ≐ G^T(θ,p,q) A has identical off-diagonal elements, A'_pq=A'_qp:

  tan(θ) = (A_pq - A_qp)/(A_qq + A_pp)

  (check that this is correct).
• J ≐ J(θ,p,q) is the Jacobi matrix where the angle is chosen to eliminate the off-diagonal elements A'_pq=A'_qp.

The matrices U and V are updated after each transformation as

U ← U G J , V ← V J .

Of course you should not actually multiply Jacobi matrices—which costs O(n³) operations—but rather perform updates which only cost O(n) operations.