Incorrect results from pseudo_inverse

[Ralith]

The matrix constructed by

na::DMatrix::from_column_slice(3, 8, &[10.0f32, 6.123234e-16, 20.0, 0.0, 10.0, -20.0, 0.0, 0.0, 0.0, 0.0, 10.0, 20.0, -10.0, 6.123234e-16, 20.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])

produces a pseudo_inverse of approximately:

  ┌                      ┐
  │  0.050 -0.000 -0.013 │
  │  0.000 -0.050  0.013 │
  │  0.000 -0.000 -0.000 │
  │  0.000 -0.050 -0.013 │
  │ -0.050 -0.000 -0.013 │
  │  0.000 -0.000  0.000 │
  │  0.000 -0.000  0.000 │
  │  0.000 -0.000  0.000 │
  └                      ┘

If f32 is replaced with f64, or if the entries having value 6.123234e-16 are zeroed out, pseudo_inverse instead yields approximately:

  ┌                      ┐
  │  0.050  0.000  0.013 │
  │  0.000  0.050 -0.013 │
  │  0.000  0.000  0.000 │
  │  0.000  0.050  0.013 │
  │ -0.050  0.000  0.013 │
  │  0.000  0.000  0.000 │
  │  0.000  0.000  0.000 │
  │  0.000  0.000  0.000 │
  └                      ┘

Note that the signs of the second and third columns are inverted. This result is the expected one. This is a massive change for a tiny perturbation of the input, and is breaking the algorithm I have built around this operation. Is this expected?

[rasmusgo]

The result is not a proper pseudo inverse for this matrix. Here is a test that fails for the given matrix:

#[test]
fn test_pseudo_inverse() {
    let s = 6.123234e-16;
    // let s = 0.0;
    let m = nalgebra::DMatrix::from_column_slice(
        3,
        8,
        &[
            10.0f32, s, 20.0, 0.0, 10.0, -20.0, 0.0, 0.0, 0.0, 0.0, 10.0, 20.0, -10.0, s, 20.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
        ],
    );
    let m_inv = m.clone().pseudo_inverse(1.0e-6).unwrap();
    println!("m: {:}", m);
    println!("m†: {:}", &m_inv);
    println!("mm†: {:}", &m * &m_inv);
    println!("mm†m: {:}", &m * &m_inv * &m);
    // Assert distance to identity matrix is small
    assert!(
        (&m * &m_inv).apply_metric_distance(&nalgebra::Matrix3::identity(), &nalgebra::UniformNorm)
            < 0.0001
    );
}

Also available as two tests in this rust playground.

mm† gets two -1 instead of 1 on the diagonal.

mm†: 
  ┌                                                                                                          ┐
  │                                  1                                  0                                  0 │
  │ -0.0000000000000000000000033087225                                 -1        -0.000000000000000007654043 │
  │                 -0.000000059604645                 -0.000000059604645                         -1.0000001 │

[rasmusgo]

Here is a test demonstrating that the problem comes from the SVD:

#[test]
// Exercises bug reported in issue #1313 of nalgebra (https://github.com/dimforge/nalgebra/issues/1313)
fn svd_regression_issue_1313() {
    let s = 6.123234e-16_f32;
    let m = nalgebra::dmatrix![
        10.0,   0.0, 0.0,  0.0, -10.0, 0.0, 0.0, 0.0;
           s,  10.0, 0.0, 10.0,     s, 0.0, 0.0, 0.0;
        20.0, -20.0, 0.0, 20.0,  20.0, 0.0, 0.0, 0.0;
    ];
    let svd = m.clone().svd(true, true);
    let m2 = svd.recompose().unwrap();
    assert_relative_eq!(&m, &m2, epsilon = 1e-5);
}

Printing the matrices (rounded to one decimal) reveals that the last row has been negated:

m: 
  ┌                                                 ┐
  │  10.0   0.0   0.0   0.0 -10.0   0.0   0.0   0.0 │
  │   0.0  10.0   0.0  10.0   0.0   0.0   0.0   0.0 │
  │  20.0 -20.0   0.0  20.0  20.0   0.0   0.0   0.0 │
  └                                                 ┘


m2: 
  ┌                                                 ┐
  │  10.0   0.0   0.0   0.0 -10.0   0.0   0.0   0.0 │
  │  -0.0 -10.0  -0.0 -10.0  -0.0  -0.0  -0.0  -0.0 │
  │ -20.0  20.0  -0.0 -20.0 -20.0   0.0   0.0   0.0 │
  └                                                 ┘

[rasmusgo]

The problem seems to be in the bidiagonalization:

#[test]
fn bidiagonal_regression_issue_1313() {
    let s = 6.123234e-16_f32;
    let mut m = nalgebra::dmatrix![
        10.0,   0.0, 0.0,  0.0, -10.0, 0.0, 0.0, 0.0;
        s,     10.0, 0.0, 10.0,     s, 0.0, 0.0, 0.0;
        20.0, -20.0, 0.0, 20.0,  20.0, 0.0, 0.0, 0.0;
    ];
    m.unscale_mut(m.camax());
    let bidiagonal = m.clone().bidiagonalize();
    let (u, d, v_t) = bidiagonal.unpack();
    let m2 = &u * d * &v_t;
    println!("m = {:.1}", m);
    println!("m2 = {:.1}", m2);
    assert_relative_eq!(m, m2, epsilon = 1e-6);
}

m = 
  ┌                                         ┐
  │  0.5  0.0  0.0  0.0 -0.5  0.0  0.0  0.0 │
  │  0.0  0.5  0.0  0.5  0.0  0.0  0.0  0.0 │
  │  1.0 -1.0  0.0  1.0  1.0  0.0  0.0  0.0 │
  └                                         ┘


m2 = 
  ┌                                         ┐
  │  0.5  0.0  0.0  0.0 -0.5  0.0  0.0  0.0 │
  │ -0.0 -0.5 -0.0 -0.5 -0.0 -0.0 -0.0 -0.0 │
  │ -1.0  1.0 -0.0 -1.0 -1.0  0.0  0.0  0.0 │
  └                                         ┘

[Ralith]

Fixed by #1314.