I have an implementation of the Golub-Kahan algorithm which sort of works but I recently discovered a bug. I've provided some sources at the end which I've used for my implementation. I've also linked my source in case it is of interest to any other readers.
Part of the SVD algorithm is to ensure that the matrix is in irreducible bidiagonal form. This means that if we have a zero in the diagonal with a non-zero super-diagonal we must zero the row and split the matrix into two bidiagonal parts. This is described on page 454 of the Matrix Computations book (linked below).
For example:
$$ \begin{bmatrix} b11 & b12 & 0 \\ 0 & 0 & b23 \\ 0 & 0 & b33 \\ \end{bmatrix} $$
Here $b22$ is $0$ and so we must zero the second row. In my current implementation I am essentially just setting $b23 = 0$, which I realise now is foolish. Instead I should be using Given's rotations to zero this row. In the case above I understand how to do this - but I am unsure how to do this in more complex cases. For example:
$$ \begin{bmatrix} b11 & b12 & 0 & 0 & 0 \\ 0 & 0 & b23 & 0 & 0 \\ 0 & 0 & b33 & b34 & 0 \\ 0 & 0 & 0 & b44 & b45 \\ 0 & 0 & 0 & 0 & b55 \\ \end{bmatrix} $$
Once again I need to zero the second row. With one Given's rotation acting on rows 2 and 3 I can get:
$$ \begin{bmatrix} b11 & b12 & 0 & 0 & 0 \\ 0 & 0 & 0 & x & 0 \\ 0 & 0 & y & z & 0 \\ 0 & 0 & 0 & b44 & b45 \\ 0 & 0 & 0 & 0 & b55 \\ \end{bmatrix} $$
But after this how can I guarantee that any rotation will keep $b23 = 0$? In the Matrix Computations book they describe the next step as:
$$ \begin{bmatrix} b11 & b12 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & x' \\ 0 & 0 & y' & z' & 0 \\ 0 & 0 & 0 & b44 & b45 \\ 0 & 0 & 0 & 0 & b55 \\ \end{bmatrix} $$
But I can't figure out a Given's rotation that would get me here.
As an aside, but still related to this algorithm, should I expect the algorithm to produce non-negative, ordered singular values? This is not the case currently but perhaps fixing this step will fix this?
[1] - Computation of the Singular Value Decomposition, algorithm 1b
So I managed to figure this one out. My confusion stemmed from the fact that I thought the Given's rotation had to be acting on adjacent rows. To go from the second matrix in the question body to the third - we use a Given's rotation acting on the 2nd and 4th rows:
$$ G = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & c(\theta) & 0 & -s(\theta) & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & s(\theta) & 0 & c(\theta) & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} $$
So that the 1st, 3rd and 5th rows remain as they are. And we choose $\theta$ s.t the 4th element in row 2 is zeroed.
$$ G * B = \begin{bmatrix} b11 & b12 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & x' \\ 0 & 0 & y & z & 0 \\ 0 & 0 & 0 & u & v \\ 0 & 0 & 0 & 0 & b55 \\ \end{bmatrix} $$
Finally to zero $x'$ we repeat the above procedure with the 2nd and 5th rows.
I haven't implemented this yet but I'm pretty confident it works!