rulinalg

rulinalg copied to clipboard

Somewhat related to #17 and #46.

There are two known bugs in the SVD implementation.

• [x] When computing the bidiagonal blocks we use the wrong criteria.

// This is wrong:
diag_abs_sum = T::min_positive_value() *
                                   (b_ii.abs() + *b.get_unchecked([i + 1, i + 1]));

// Instead:
diag_abs_sum = T::min_positive_value() *
                                   (b_ii.abs() + b.get_unchecked([i + 1, i + 1]).abs());

• [ ] We zero the off-diagonals incorrectly

We also zero the off-diagonals incorrectly before running the Golub-Kahan step. Currently we use a givens rotation to zero the off diagonal element in b. This is essentially equivalent to setting the value to zero directly... Instead we should use a (implicit) Given's rotation on the whole matrix and update the orthogonal matrices in turn.

• [ ] We do not need to transpose b if it was flipped before

Even if we transposed at the start, the matrix b is diagonal by definition and so we gain nothing by transposing it.

• [ ] Our stopping criterion is far too harsh

We frequently check that we are below T::min_positive_value(). This makes the algorithm take a very long time to converge. Instead we should use a larger (but still small) value.

In addition to these known issues we should explore the following:

• [x] Should the algorithm order the singular values automatically?
• [x] Why does the algorithm return negative singular values?

These last two may not be bugs but rather things we clean up after the algorithm.

I had a little time to look into this issue. It's really hard to find details of the implementation anywhere - but I checked out the netlib/eispack code. Their algorithm returns the singular values out of order but they are non-negative. Moreover it looks like they enforce the non-negative behaviour explicitly (you can see near the end of the file).

I had a little time to look into this issue. It's really hard to find details of the implementation anywhere - but I checked out the netlib/eispack code. Their algorithm returns the singular values out of order but they are non-negative. Moreover it looks like they enforce the non-negative behaviour explicitly (you can see near the end of the file).

Well, that certainly is no easy read! With all this goto jumping around, I can't tell from a quick glance what their criterion for flipping the sign is, but the fact that they do this suggests that you seem to be right in your assumption that the Golub-Reinsch algorithm seems to be able to return negative singular values. That's some good detective work right there!

By the way, I'll just leave this here as it's something we should consider.

Currently, the SVD implementation makes several tests against T::min_positive_value(). As you stated in #46, it's preferable to use the machine epsilon for the specific type instead, which is why you've made a PR to the num crate. However, I think we should probably multiply this factor by a small number (somewhere between 2-9, maybe?), because if we check against macheps directly, I think it's possible we will make a number of redundant iterations due to lack of precision. The "small multiple" is also very briefly mentioned in the paper you've previously mentioned.

Ticked off the problems tackled in #46 .

@AtheMathmo: The stopping criterion was also "fixed" in #46, in the sense that it uses 3 * eps where eps is machine epsilon instead of T::min_positive_value().

Thanks! Ticked that off too