This post reminds quickly how to diagonalize a dense matrix A using the Singular Value Decomposition (SVD).

Let $A$ a $m_Xn$ real matrix, it always exist the decomposition: $A = U\Sigma V^T$, with $U$ a $m_Xm$ orthogonal matrix, $V$ ($n_Xn$) another orthogonal matrix and $\Sigma$ a $m_Xn$ diagonal matrix on the top $m-n$ rows.

To solve that there are two main approaches:

- The
**Jacobi based SVD**method which is the simplest, - The
**QR Decomposition based SVD**method which is computationally more efficient.

# Jacobi based SVD

$A$ is a symmetric matrix and we apply plane rotations $J^{(0)}$… $J^{(n)}$ in order to diagonalize iteratively the $A$ matrix :

$$

A^{(k+1)} = J^{(k)T}A^{(k)} J^{(k)}

$$

with $A^{(0)} = A$, and $A^{(\infty)} = \Sigma$

The $J$ matrix is chosen orthogonal to ensure that $J^{(0)}.J^{(1)}$… $J^{(n)}$ stays orthogonal, thus has the form:

$$

J(i,j,\phi) =

\begin{pmatrix}

1 & & & & & & & \\

& \ddots & & & & & & \\

& & c & & s & & & \\

& & & \ddots & & & & \\

& & -s & & c & & & \\

& & & & & \ddots & & \\

& & & & & & 1 & \\

\end{pmatrix}

$$

intersection of rows $i$, $j$ and columns $i$, $j$ are filled with a plane rotation matrix.

The goal is to find the rotation matrix J such that non-diagonal elements of $A$ become 0. The parameter to compute is $\phi$ and can be found by solving:

\[

\begin{pmatrix}

\hat{x} & 0 \\

0 & \hat{y}

\end{pmatrix} =

\begin{pmatrix}

c & -s \\

s & c

\end{pmatrix}

\begin{pmatrix}

x & w \\

w & y

\end{pmatrix}

\begin{pmatrix}

c & s \\

-s & c

\end{pmatrix}

\]

that leads to an equation: $w(c^2-s^2) + cs(x-y) = 0$ that gives the $\phi$ angle (for more details, the reader should refer to “A Stream Algorithm for the SVD“ Volker Strumpen, Henry Hormann, and Anant Agarwal).

Note that this is also possible to decompose and diagonalize non symmetric matrices by adding a step to make the $i$, $j$ rows and columns to become symmetric first. Then one can apply previous computations.

**QR Decomposition based SVD**

This principle is more efficient than the previous one. It is based on the QR Decomposition method to first obtain a bi-diagonal matrix B. Then, in a second step, this B matrix is transformed using the Golub-Kahan-Reinsch algorithm. This method is now a standard when solving the SVD problems. However it seems that it does not apply very well to parallel programming.