One-sided Jacobi algorithm for Singular Value Decomposition

Introduction

The singular value decomposition (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 one-sided Jacobi SVD algorithm.

Algorithm

In this method the elementary iteration is given as

A → A J(θ,p,q)

where the indices (p,q) are swept cyclicly (p=1..n, q=p+1..n) and where the angle θ is chosen such that the columns number p and q of the matrix AJ(θ,p,q) are orthogonal. One can show that the angle should be taken from the following equation (you should use atan2 function),

tan(2θ)=2a_p^Ta_q /(a_q^Ta_q - a_p^Ta_p)

where a_i is the i-th column of matrix A (check this).

After the iterations converge and the matrix A'=AJ (where J is the accumulation of the individual rotations) has orthogonal columns, the SVD is simply given as

A=UDV^T

where

V=J, D_ii=||a'_i||, u_i=a'_i/||a'_i||,

where a'_i is the i-th column of matrix A' and u_i us the i-th column of matrix U.