Monday, April 18, 2011

Two diagonals

There is a matrix bidiagonalization algorithm in a good book Golub G.H., Van Loan C.F. Matrix Computation. 2nd edition. The John Hopkins Univ. Press, 1989, which relies on the Householder's reflections to create zeroes below the lower triangle part of any matrix and above the upper subdiagonal so that the resulting matrix would be something like:
This matrix is usually be used for the Singular Value Decomposition by elimination of values above the main diagonal (so that the resulting matrix performes a required scaling). There are two interesting moments:


------------------------------------
1. The algorithm 5.4.2, which performs bidiagonalization, in the book seems a bit wrong. Namely two lines can be omitted without harming the result. Those are: A(j+1:m,j) = v(j+1:m) and A(j,j+1:n) = v(j+2:m)' (in Matlab notation). These lines perform copying of the Householder's vector into the matrix, but this is not required because aff ected column and row already nulli fied by the Housholder transform and are not changed as the algorithm goes on. So these copyings simply do something wrong. Here is an example from the book: for input matrix

the following bidiagonalized matrix should be obtained:
But implementation of the algorithm yields the following result:

Note that values on the main diagonal and the upper sub-diagonal are almost the same, while there are a lot of "undesired" non-zeros as well (the di fferences are due to use of single-precision calculations). There are also signs misplaced, but I haven't found the source of this error (?) yet.

Now removing the copying of the Householder's vectors we obtain:


------------------------------------

2. I suppose that there is a good way for calculation of the matrix rank using this bidiagonalization procedure. There a lot of ways to calculate the rank, the simpliest one is to use Gauss eliminations, but it's not numerically stable, and a more reliable and convenient way is to use SVD and to calculate singular values larger than some threshold. One of the typical ways to calculate SVD is to use Golub-Reinsch algorithm, which firstly transforms the input matrix to the bidiagonal form (!) and then eliminates non-zeroes above the main diagonal. My idea is that the values on the main diagonal of bidiagonalized matrix can tell the number of singular values. If it's true then the matrix rank calculation can be performed via using bidiagonalization and then calculating non-zero values on the main diagonal. For the example above the initial matrix rank is 2 and the number of non-zeroes on the main diagonal of bidiagonalized form is also two.
------------------------------------


I may be wrong for these two notes, so any comments and critics are wellcome :).