SVD: Occasional loss of accuracy

[robertsonj3]

Hi, I haven't been using nalgebra long, and I am not an expert with SVD functionality, so I'm not 100% sure if this is an issue with nalgebra or me.

I am using SVD for alignment of 3d points using random point clouds, and random transformations (the euler angle is random, the rotation matrix is valid). Occasionally, the alignment is poor.

I have traced this back to what appears to be low accuracy in the SVD result, that is, u s v_t does not yield the same matrix that the svd was performed on.

Example: H: │ 0.24777976656646927 0.5496207895780494 -0.7403583528727136 │ │ 0.5832953614864118 -0.38015695826419854 -0.03392366142812886 │ │ -0.3738110622901052 0.359160059894174 -0.09620293929130566 │

nalgebra svd produced: u: │ -0.873792828198889 0.47077673362678457 0.1218825683347777 │ │ 0.31813462431343165 0.7429469762516004 -0.5889143836684206 │ │ -0.3677994755313523 -0.4758140997850944 -0.7989521188685604 │

s: │ 0.9999999999999996 0 0 │ │ 0 0.8220531474729237 0 │ │ 0 0 0.00000000000000012761238385273053 │

v (this is the transpose of v_t produced by the svd): │ 0.1065457803497895 0.885789524773893 0.45169117158794986 │ │ -0.7332946769207143 -0.23679742707193852 0.6373428397117634 │ │ 0.6715109183694628 -0.3991289219449654 0.624314976736623 │

u s v_t: │ 0.24970470373979137 0.5491061976155459 -0.7412257125301849 │ │ 0.574884679370761 -0.3779085371334399 -0.030133882488828305 │ │ -0.3856591637197468 0.36232740366575733 -0.09086428962516223 │

The derived matrix is in error by the 3rd and 4th decimal.

I have tested this same matrix H in Octave (free Matlab) and the SVD results were different but yielded an exact match for the original H matrix when re-combined. I know that SVD can yield different results, but I understood that u s v_t will always re-produce the original matrix.

Note: although the copy from octave below only shows 6dp, it was loaded with the full precision presented above, as printed from rust debug.

Octave results for comparison, accurate to at least 6 figures: H =

0.247780 0.549621 -0.740358 0.583295 -0.380157 -0.033924 -0.373811 0.359160 -0.096203

[U,S,V] = svd(H) U =

0.87379 0.46823 -0.13133 -0.31813 0.75465 0.57384 0.36780 -0.45964 0.80837

S =

Diagonal Matrix

1.0000e+00 0 0 0 8.2189e-01 0 0 0 3.6928e-17

V =

-0.10655 0.88579 0.45169 0.73329 -0.23680 0.63734 -0.67151 -0.39913 0.62431

U S transpose(V) ans =

0.247780 0.549621 -0.740358 0.583295 -0.380157 -0.033924 -0.373811 0.359160 -0.096203

Thanks for any help with this.

[allan2]

Could you provide the code for this example?

[robertsonj3]

main.zip Sorry, the codes a mess. You'll probably be most interested in function find_transformation. It won't reproduce the exact sample I posted earlier, but usually gets several hits per run.

I am cleaning up the code, I'll post a new sample when I'm finished.

[robertsonj3]

Just re-reading my last post, and I realized it might not be clear that I did actually post a sample.

[robertsonj3]

Here is a slightly cleaner version of the code exhibiting the same behavior. three_d_operations.zip

[robertsonj3]

Hi Allan, did you find any issues with he SVD feature?

[robertsonj3]

``Hi All,

This issue seems to have gone dormant, but I have just tested with the latest nalgebra and as far as I can tell the behavior has not changed.

nalgebra = "0.32.3" Rust version: 1.73.0 (cc66ad468 2023-10-03)

I have created a new rust project, attached, to demonstrate the issue, and it includes clear samples of matrix that work well and work poorly. Along with this I have included copies of the projects command line output.

I have also included a MATLAB file showing the same matrix can be decomposed and recombined accurately. I ran the file in Octave rather than MATLAB as I don't have access to MATLAB. Again I have included a copy of the results.

svd_demo.zip

main.rs

// Demonstrate how the nalgebra SVD function can yield poor accuracy results
// This file contains a set of matrices that are known to pass and fail
//  (fail has been defined where any matrix element deviates by more than 0.0001 when reconstructed from the u, s & v decomposition)
// Results will be written to the command line for each matrix as: success - original matrix - reconstructed matrix - deviation matrix
// Note that all matrices are 3x3, I have not observed this in other dimensions
// Generally 0.08% of matricies I have tested result in poor accuracy

use nalgebra::DMatrix;
use std::num::ParseIntError;

fn main() {

    // Test matrices where SVD does NOT yield suitable accuracy
    let breaking_values = failing_mat();
    for matrix in breaking_values.into_iter(){
        test_svd(matrix);
    }

    // Test matrices where SVD does yield suitable accuracy
    let working_matices = working_mat();
    for matrix in working_matices.into_iter(){
        test_svd(matrix);
    }

}

// Retrieve matrix that yield an accurate SVD - as hex strings
// Stored matrices - (r:3, c:3) row major - f64 hex big endian
fn working_mat() -> Vec<DMatrix<f64>> {
    let matrices = vec![
        vec!["bfd006ac63b40e10","bfc3ab314a6ec61d","bfc899de29ec8743","bfa917a53710fed7","bfadeb4364115031","bfb8f8c624152c92","bfe795f55efb725e","bfd889dc052eea54","bfdae40a6f5173fe"],
        vec!["3fbcd57cb673f3ac","3fb49223b6e644a5","bf93c6b5d66aa0bf","bfe3b36f41a123b8","bfdda90174dffdf2","3fbaccedbe038737","3fdba4c532c5ae6e","3fdbd8b357eb9272","bfb1c9491e77042e"],
    ];

    matrices.into_iter().map(|x| get_matrix(&x)).collect()
}

// Retrieve matrix that yield an accurate SVD - as hex strings
// Stored matrices - (r:3, c:3) row major - f64 hex big endian
fn failing_mat() -> Vec<DMatrix<f64>> {
    let matrices = vec![
        vec!["bfdfb839adcb0d65","bfd321a3cac82310","3fb2c901e3881e76","3fb959e84cffb53c","bfd923574c548c76","3fe981e5460e628b","bfe7593e647c9500","bfd55218cb48782c","bfb57515d9d721ea"],
        vec!["bfda0bf3eb4e5e1f","bfd02ec349ca53a6","bfe1c0d9dc202bd9","3fe1a377131e0780","bfe7ab8801039f07","bfd1b76e89bc11be","bfdb627dc0bf7d18","3fb5af7472630848","bfd0037068afd465"],
    ];

    matrices.into_iter().map(|x| get_matrix(&x)).collect()
}

// Test nalgebra SVD function using the supplied matrix
fn test_svd(mat: DMatrix<f64>) -> bool {

    // Perform SVD function on the matrix
    let svd = mat.clone().svd(true, true);

    let v_t = svd.v_t.expect("");
    let u = svd.u.expect("");
    let s = DMatrix::from_fn(3, 3, |r, c| if r == c { svd.singular_values[(r,0)] } else { 0.0 });

    // Reconstruct the h matrix from the SVD result
    let mat2 = u * s * v_t;

    // Compare the original and reconstructed matrix
    let result = mat2.relative_eq(&mat, 0.0001, 0.0001);

    // Print the result
    println!("Recreated original matrix: {}\nMatrix 1:{}\nMatrix 2:{}Difference:{}", result, mat.clone(), mat2.clone(), mat - mat2);

    result
}

// build a nalgebra matrix from a vector of f64 values in hex string representation
fn get_matrix(values: &Vec<&str>) -> DMatrix<f64> {

    let mut mat = DMatrix::zeros(3,3);

    for c in 0..3 {
        for r in 0..3 {
            mat[(r,c)] = decode_hex(values[c*3+r]);
        }
    }

    mat
}

// Decode a hex string into an f64
pub fn decode_hex(s: &str) -> f64 {
    let val: Result<Vec<u8>, ParseIntError> = (0..s.len())
        .step_by(2)
        .map(|i| u8::from_str_radix(&s[i..i + 2], 16))
        .collect();

    let val = val.expect("");
    let val: [u8; 8] = val.try_into().expect("");
    f64::from_be_bytes(val)
}

MATLAB script

X = [{'bfdfb839adcb0d65'; 'bfd321a3cac82310'; '3fb2c901e3881e76'; '3fb959e84cffb53c'; 'bfd923574c548c76'; '3fe981e5460e628b'; 'bfe7593e647c9500'; 'bfd55218cb48782c'; 'bfb57515d9d721ea'},{'bfda0bf3eb4e5e1f';'bfd02ec349ca53a6';'bfe1c0d9dc202bd9';'3fe1a377131e0780';'bfe7ab8801039f07';'bfd1b76e89bc11be';'bfdb627dc0bf7d18';'3fb5af7472630848';'bfd0037068afd465'},{'bfd006ac63b40e10';'bfc3ab314a6ec61d';'bfc899de29ec8743';'bfa917a53710fed7';'bfadeb4364115031';'bfb8f8c624152c92';'bfe795f55efb725e';'bfd889dc052eea54';'bfdae40a6f5173fe'},{'3fbcd57cb673f3ac';'3fb49223b6e644a5';'bf93c6b5d66aa0bf';'bfe3b36f41a123b8';'bfdda90174dffdf2';'3fbaccedbe038737';'3fdba4c532c5ae6e';'3fdbd8b357eb9272';'bfb1c9491e77042e'}]

Y =zeros(3,3);
for N = 0:3
  for col = 0:2
    for row = 0:2
        Y(row + 1, col + 1)=hex2num(X{N*9+col*3+row+1});
    end
  end

  Original = Y

  [U, S, V] = svd (Y);

  Reconstructed = U * S * V'

  Diff = Original - Reconstructed

end

[robertsonj3]

I have reduced the example for this issue. I can update the test cases to flag this issue, but I don't know the SVD algorithm well enough to find the issue. Can anyone help with this issue or point me toward the documentation for the algorithm.

#[macro_use]
extern crate approx;

fn main() {
    let m = nalgebra::dmatrix![
        -0.49561922046651424, 0.09902812843184433, -0.7296440088872771;
        -0.29892821123989766, -0.3927820439732338, -0.3331357941805233;
        0.07337962918256288,  0.7971063965115531, -0.0838178307377643
    ];

    let svd = m.clone().svd(true, true);
    let m2 = svd.recompose().unwrap();
    assert_relative_eq!(&m, &m2, epsilon = 1e-5);
}

The epsilon value was selected to match the test cases in nalgebra.