HESTENES' METHOD

To compute the SVD of A, an orthogonal matrix V is generated such that the transformed matrix AV = W has orthogonal columns. Scaling each column of W so that its 2-norm is unity we get :

$$W = \tilde{U} \tilde{\Sigma}$$

where $\tilde{U}$ is an $m \times n$ matrix and $\tilde{\Sigma}$ is an $n \times n$ nonnegative diagonal matrix. The corresponding SVD of A is given by :

 

$$A = \tilde{U} \tilde{\Sigma} V^T$$ (2)

Comparing equation (2) with equation (1), we notice the different sizes of U and $\tilde{U}$ as well as $\Sigma$ and $\tilde{\Sigma}$. However, since the null columns of U and $\tilde{U}$ are associated with zero diagonal elements of $\Sigma$ and $\tilde{\Sigma}$, respectively, the non-null columns of U and $\tilde{U}$ form orthogonal bases for range(A) and there is no essential difference between equations (1) and (2). Successive Jacobi rotations are applied to A in order to obtain W, giving us the following iteration :

A_k+1 = A_k J_k

where $J_{k} = J(p,q,\theta)$ and $\theta$ is determined so that the columns of p and q of A_k+1 are orthogonal. Denoting column p of A_k as a^(k)_p, the requirements on $\theta$ can be stated as:

$$\left( \begin{array} {cc} a^{(k+1)}_{p}&a^{(k+1)}_{q}\\ \end... ...ight) \left( \begin{array} {cc} c&s\\ -s&c\\ \end{array}\right)$$

where c and s are the cosine and sine of $\theta$, such that a^(k+1)T_p a^(k+1)_q = 0, leading to the following determination of c and s :

a^(k+1)_p = c a^(k)_p - s a^(k)_q

a^(k+1)_q = s a^(k)_p + c a^(k)_q

a^(k+1)T_p a^(k+1)_q = c s a^(k)T_p a^(k)_p + c² a^(k)T_p a^(k)_q - s² a^(k)T_q a^(k)_p - c s a^(k)T_q a^(k)_q = 0

Denoting $\alpha = a^{(k)T}_{p} a^{(k)}_{p}, \beta = a^{(k)T}_{q} a^{(k)}_{q}, \gamma = a^{(k)T}_{p} a^{(k)}_{q} = a^{(k)T}_{q} a^{(k)}_{p}$ we get :

$$c s \alpha + (c^{2} - s^{2}) \gamma - c s \beta = 0$$

Dividing by c s :

$$\alpha + \left( \frac{c}{s}-\frac{s}{c} \right) \gamma - \beta = 0$$

Denoting $t = \frac{s}{c}$ and $\xi = \frac{\beta - \alpha}{2 \gamma}$ we get :

$$t^{2} + 2 \xi t - 1 = 0$$

giving us $t = \frac{sign(\xi)}{\vert\xi\vert+\sqrt{1+\xi^{2}}}, c = \frac{1}{\sqrt{1+t^{2}}}$ and s = t c.

For every iteration we choose different column pairs (p,q) and obtain a matrix J_k that will orthogonalize those columns. As we are computing W, we can also compute V beginning with the $n \times n$ identity matrix and performing the following with each iteration :

V_k+1 = V_k J_k

There are several ways of choosing the pairs of columns (p,q) that are orthogonalized at every step. The approach used here is one that is suitable for parallel implementation on the CM.