SVD produces wrong singular values

[cfunky]

The SVD with the default tolerance parameter eps can produce significantly wrong singular values. This behavior occurs in nalgebra version 0.30.1, potentially also in other versions (I did not try this). It looks like the problem occurs when, as in the example above, an singular value is close to zero and the tolerance parameter eps of Matrix::try_svd(...) or Matrix::try_svd_unordered(...) is set to a very low value, as is done internally, e.g., in Matrix::svd(...). Here's a minimal example reproducing this issue:

use nalgebra::{dmatrix, dvector};

fn main() {
    let x = dmatrix![-6.206610118536945f64, -3.67612186839874; -1.2755730783423473, 6.047238193479124];
    // let mut x_svd = x.try_svd(true, true, 1e-8, 0).unwrap();// This gives almost the same result as the line below.
    let mut x_svd = x.svd(true, true); // This gives almost the same result as the line above.
    println!("{}", x_svd.singular_values);
    x_svd.singular_values = dvector![1.0, 0.0];
    println!("{}", x_svd.singular_values);
    let y = x_svd.recompose().unwrap();
    let y_svd = y.svd(true, true); // This produces wrong singular values!
    // let y_svd = y.try_svd(true, true, 1e-8, 0).unwrap(); // This produces correct singular values! For small enough `eps` the wrong behavior occurs.
    println!("{}", y_svd.singular_values);
}

Link to playground

The above code prints

  ┌                   ┐
  │  7.81117060812432 │
  │ 5.405337300528416 │
  └                   ┘

  ┌   ┐
  │ 1 │
  │ 0 │
  └   ┘

  ┌                                    ┐
  │                 1.1503016006828526 │
  │ 0.00000000000000015415627596860353 │
  └                                    ┘

[Andlon]

Thanks for reporting! Ugh, that's pretty bad. Here's a modified version which also compiles for older versions of nalgebra:

use nalgebra::{DMatrix, DVector};

fn main() {
    let x = DMatrix::from_row_slice(2, 2, &[-6.206610118536945f64, -3.67612186839874, -1.2755730783423473, 6.047238193479124]);
    // let mut x_svd = x.try_svd(true, true, 1e-8, 0).unwrap();// This gives almost the same result as the line below.
    let mut x_svd = x.svd(true, true); // This gives almost the same result as the line above.
    println!("Original singular values: {}", x_svd.singular_values);
    x_svd.singular_values = DVector::from_row_slice(&[1.0, 0.0]);
    println!("Modified singular values: {}", x_svd.singular_values);
    let y = x_svd.recompose().unwrap();
    let y_svd = y.clone().svd(true, true); // This produces wrong singular values!
    println!("Recomposed singular values (default eps): {}", y_svd.singular_values);
    let y_svd = y.try_svd(true, true, 1e-8, 0).unwrap(); // This produces correct singular values! For small enough `eps` the wrong behavior occurs.
    println!("Recomposed singular values (eps = 1e-8): {}", y_svd.singular_values);
}

It seems to always have been broken, although the incorrect exact values produced seems to have changed a little between 0.29 and 0.30, probably due to the 2x2 and 3x3 special cases added in that release.

[Andlon]

Update: it seems setting eps = 1e-15 also fixes the issue in this case. I suspect the default epsilon is simply too close to unit round-off and something gets numerically completely truncated... Of course, I'm not sure how general this observation is. Surely it can be expected to depend on conditioning etc. I don't know the details of the algorithm well enough.

[Andlon]

I briefly looked at the SVD implementation. I believe it's probably based off of the SVD algorithm lined out in the book Matrix Computations by Golub & Van Loan. I couldn't find any comments though... @sebcrozet can hopefully clarify that. Here's an excerpt:

Note that it says that eps should be a small multiple of the unit round-off. The SVD impl uses T::RealField::default_epsilon(), which is exactly the unit round-off (~2.22e-16). So I guess we should multiply this at least by a factor of 2-4.

EDIT: Simply multiplying by a factor 2 seems to be sufficient in this case. I'd probably prefer something like 4 to have a bit more leeway, however. Is there a way we can test this problem more thoroughly?

[cfunky]

Increasing eps indeed seems to solve the problem in this case and also for matrices of different dimensions. Using the proptest crate and running the following code

use nalgebra::proptest::matrix;
use nalgebra::DVector;
use proptest::prelude::*;

proptest! {
    fn test_singular_values(x in matrix(-10f64..=10f64, 1..=5, 1..=5)) {
        let mut x_svd = x.svd(true, true);
        let mut true_singular_vals = DVector::zeros(x_svd.singular_values.len());
        true_singular_vals[0] = 1.0;
        x_svd.singular_values = true_singular_vals.clone();
        let y = x_svd.recompose().unwrap();
        let y_svd = y.try_svd(true, true, 4.0 * f64::EPSILON, 0).unwrap();
        prop_assert!(y_svd.singular_values.relative_eq(&true_singular_vals, 1e-12, 1e-12));
    }
}

fn main() {
    test_singular_values();
}

produces correct results.

Because you mentioned the SVD algorithm from the book by Golub & Van Loan, I took a look at the GNU scientific library (written in C) that uses this algorithm to compute the SVD. It seems that in their implementation the eps parameter is set to the machine epsilon and that, at least judging by my experiments, it works just fine. I think this might just be due to implementation differences (assuming nalgebra even uses this algorithm) leading to slightly different numerical values or it could point to a more serious flaw in the nalgebra implementation (possibly related to #983).

[sebcrozet]

nalgebra does implement the algorithm from Golub & Van Loan. So we should modify our epsilon to match the book’s recommendation.

[Andlon]

@cfunky: I had forgotten about or entirely missed #983. I would suspect there's more to the story here than just the eps parameter, especially in light of your observations from the GNU scientific library.

Given we rely on SVD for our work, this is personally very important to me to resolve, but unfortunately I also don't have the bandwidth to dig into this any time soon. Any assistance here would be much appreciated!

[cfunky]

Sorry for the delayed response. I'd be interested in investigating this further, as I also rely on nalgebra for work. I probably won't have a lot of time in the next couple of weeks due to upcoming deadlines, but would be happy to take a closer look afterwards.

[Andlon]

I'm somewhat cautiously closing this as I hope/believe it has been addressed by #1089. Feel free to re-open or make a new issue if similar issues crop up again, however!