nalgebra icon indicating copy to clipboard operation
nalgebra copied to clipboard

SVD fails to order singular values for 3x3 matrix

Open rasmusgo opened this issue 2 years ago • 3 comments

I had a bug in my physics simulation (not rapier) that I tracked down to relying on svd to provide singular values in sorted order. The doc for svd states:

Computes the Singular Value Decomposition using implicit shift. The singular values are guaranteed to be sorted in descending order. If this order is not required consider using svd_unordered.

But the singular values I was getting were not sorted. Here is a test that reproduces the bug:

#[test]
fn svd3_sorted() {
    let mat = Matrix3::<f32>::new(
        1.0,
        0.0,
        0.0,
        0.0,
        0.99999994,
        0.0,
        0.0,
        -0.00000054691924,
        1.0000005,
    );
    let svd = mat.svd(true, true);
    assert!(
        svd.singular_values[2] <= svd.singular_values[1]
            && svd.singular_values[1] <= svd.singular_values[0],
        "Singular values are not sorted in descending order: {}",
        svd.singular_values
    );
}

This results in:

thread 'svd3_sorted' panicked at 'Singular values are not sorted in descending order: 
  ┌            ┐
  │  1.0000005 │
  │ 0.99999994 │
  │          1 │
  └            ┘

I am using nalgebra version 0.32.1.

rasmusgo avatar Mar 04 '23 20:03 rasmusgo

There is a typo in svd3.rs:45 where rho0 is swapped with rho2 but it should be rho1 that is swapped with rho2. I don't think this is the cause of the bug though.

rasmusgo avatar Mar 04 '23 21:03 rasmusgo

Here is an updated test that looks more like the tests in tests/linalg/svd.rs

#[test]
// Exercises bug reported in issue #1215 of nalgebra (https://github.com/dimforge/nalgebra/issues/1215)
fn svd_regression_issue_1215() {
    let x = nalgebra::matrix![
        1.0f32, 0.0, 0.0;
        0.0, 0.99999994, 0.0;
        0.0, -0.00000054691924, 1.0000005
    ];
    let x_svd = x.svd(true, true);
    assert!(is_sorted_descending(x_svd.singular_values.as_slice()));
}

rasmusgo avatar Mar 04 '23 22:03 rasmusgo

I also encountered this; here's my matrix:

use nalgebra;

fn main() {
    let m = nalgebra::Matrix3::new(
        1.33333492279052734375f32,
        0.666667461395263671875,
        -1.3333332538604736328125,
        0.666667461395263671875,
        0.3333337306976318359375,
        -0.66666662693023681640625,
        -1.3333332538604736328125,
        -0.66666662693023681640625,
        1.333331584930419921875,
    );
    let s = nalgebra::SVD::new(m, false, false);
    assert!(
        s.singular_values.iter().rev().is_sorted(),
        "unsorted singular values {:?}",
        s.singular_values
    );
}

mkeeter avatar Jan 11 '25 15:01 mkeeter