Singular Value Decomposition (SVD) returns wrong decomposition [bug]

[Sollimann]

PROBLEM:

Hey, I'm trying to do a singular value decomposition of a 2x2 Matrix, and it seems that the decomposition performed by nalgebra when comparing to numpy and an online calculator for SVD is wrong. I get consistent results with numpy and a online calculator, however nalgebra produce a different decomposition. To me it seems that the decomposition performed by nalgebra produces the correct numerical values, just in the wrong places and with the wrong sign. Also the decomposition by numpy and the online calculator has u and v_t being identical matrices, which is not the case in nalgebra.

To reproduce my results, you can do:

SVD IN RUST WITH NALGEBRA

    let H = na::Matrix2::new(8.75, 10.75,
                             10.75, 14.75);

    // testing both approaches
    let svd1 = na::linalg::SVD::new(H, true, true);
    let svd2 = H.svd(true, true);

   // output from both svd1 and svd2 are:
  u: Matrix { data: [[0.7964919860908672, -0.6046490850840893],
                    [0.6046490850840898, 0.7964919860908671]] }

  v_t: Matrix { data: [[0.7964919860908662, 0.60464908508409],
                    [-0.60464908508409, 0.7964919860908662]] }

  s: Matrix { data: [[0.5892428572251418, 22.910757142774862]] }

SVD IN PYTHON WITH NUMPY

I compared this to numpy in Python:

   H = np.array([[ 8.75, 10.75],
                 [10.75, 14.75]])

   U, S, Vt = np.linalg.svd(H)

   // output
    H: [[ 8.75 10.75]
       [10.75 14.75]] 

    U: [[-0.60464909 -0.79649199]
       [-0.79649199  0.60464909]] 

    Vt: [[-0.60464909 -0.79649199]
        [-0.79649199  0.60464909]] 

    S: [22.91075714  0.58924286] 

SVD WITH ONLINE CALCULATOR see link to calculator

[Andlon]

You're running into an interesting problem here. First, it is true as @lebensterben says that the SVD is not in general unique. The situation would be slightly improved if nalgebra actually sorted its singular value in descending order, which is a common convention. If we don't already have an issue for this we should probably make one.

However, what is very confusing here is that the debug output for Matrix2 outputs the matrix in column-major format. So the u and v_t matrices are in fact the transposed of what they might look like given the debug output. With this in mind, I get the following in numpy:

>>> u
array([[ 0.79649199,  0.60464909],
       [-0.60464909,  0.79649199]])
>>> v_t
array([[ 0.79649199, -0.60464909],
       [ 0.60464909,  0.79649199]])
>>> s
array([[ 0.58924286,  0.        ],
       [ 0.        , 22.91075714]])
>>> u.dot(s).dot(v_t)
array([[ 8.75, 10.75],
       [10.75, 14.75]])

Furthermore, we can verify that the u and v_t matrices are in fact orthogonal:

>>> np.linalg.det(u)
1.0000000000000009
>>> np.linalg.det(v_t)
0.9999999999999998

So it would appear that the output of nalgebra is in fact correct, it just happens to be different from that of numpy. That said, there are two things I think we should consider addressing:

1. Always return descending singular values, as this is a common convention and users might expect this.
2. Consider a custom debug implementation of statically sized matrices so as to avoid confusion with storage order.

[Sollimann]

@Andlon @lebensterben Thank you for the quick answer :) I tried to compare the debug output of matrix u and v_t and then also print each entry one by one, and I see that they are inconsistent with each other as you said:

    println!("U: {}", u);
    let data = u.data.0;
    println!("m11: {}, m12: {}, m21: {}, m22: {}", data[0][0], data[0][1], data[1][0], data[1][1]);

    println!("Vt: {}", v_t);
    let data2 = v_t.data.0;
    println!("m11: {}, m12: {}, m21: {}, m22: {}", data2[0][0], data2[0][1], data2[1][0], data2[1][1]);

    // output
U: 
  ┌                                         ┐
  │  0.7964919860908672  0.6046490850840898 │
  │ -0.6046490850840893  0.7964919860908671 │
  └                                         ┘

m11: 0.7964919860908672, m12: -0.6046490850840893, m21: 0.6046490850840898, m22: 0.7964919860908671
Vt: 
  ┌                                       ┐
  │ 0.7964919860908662  -0.60464908508409 │
  │   0.60464908508409 0.7964919860908662 │
  └                                       ┘

m11: 0.7964919860908662, m12: 0.60464908508409, m21: -0.60464908508409, m22: 0.7964919860908662

However if I try to the same for the H matrix, it seems to be inconsistent the other way around with the debug output being correct, but if I print out each entry one by one it, then it seems to be transposed.

    let H: na::Matrix2<f64> = na::Matrix2::new(8.75, -10.75,
                                             10.75, 14.75);

    println!("H: {}", H);
    let d = H.data.0;
    println!("data: {:?}", d);
    println!("m11: {}, m12: {}, m21: {}, m22: {}", d[0][0], d[0][1], d[1][0], d[1][1]);

   // output
    H: 
      ┌               ┐
      │   8.75 -10.75 │
      │  10.75  14.75 │
      └               ┘

    data: [[8.75, 10.75], [-10.75, 14.75]]
    m11: 8.75, m12: 10.75, m21: -10.75, m22: 14.75

[Sollimann]

Now that I think of it, I believe it comes down to the issue you mentioned about storage order. I find this ordering very confusing though :thinking: I would expect to find m12=-10.75 at index [0][1], but in this case it is now found at index [1][0]

[Sollimann]

So, I've been doing some experimenting. I have implemented a best fit transform in both python and in rust to compute the transformation matrix between two pointclouds. I wanted to see if the difference in the svd decomposition in numpy and svd decomposition in nalgebra would cause any different behavior, and indeed it does.

The input to the algorithm is two identical point clouds. the best fit transform between two identical point clouds should output the identity matrix for the rotation and a vector of zero length for the transformation. The Python implementation achieves this results, however the rust implementation outputs something totally different due to the differences in implementation for the svd algorithm.

To recreate my tests, see code below:

FUNCTION

def best_fit_transform(A, B):
    '''
    Calculates the least-squares best-fit transform that maps corresponding points A to B in m spatial dimensions
    Input:
      A: Nxm numpy array of corresponding points
      B: Nxm numpy array of corresponding points
    Returns:
      T: (m+1)x(m+1) homogeneous transformation matrix that maps A on to B
      R: mxm rotation matrix
      t: mx1 translation vector
    '''

    assert A.shape == B.shape

    # get number of dimensions
    m = A.shape[1]

    # compute pointcloud centroids
    centroid_A = np.mean(A, axis=0)
    centroid_B = np.mean(B, axis=0)

    # compute cross-covariance
    AA = A - centroid_A
    BB = B - centroid_B
    H = np.dot(AA.T, BB)

    # rotation matrix
    U, S, Vt = np.linalg.svd(H)
    R = np.dot(Vt.T, U.T)

    # special reflection case
    if np.linalg.det(R) < 0:
       Vt[m-1,:] *= -1
       R = np.dot(Vt.T, U.T)

    # translation
    t = centroid_B.T - np.dot(R,centroid_A.T)

    # homogeneous transformation
    T = np.identity(m+1)
    T[:m, :m] = R
    T[:m, m] = t

    return T, R, t

MAIN

// function call
  A = np.array([[1.0, 2.0],
               [2.0, 2.0],
               [3.0, 3.0],
               [-1.0, -2.0]])

  B = np.copy(A)

  // homogeneous transformation matrix between two pointclouds
  _, R, t = best_fit_transform(A, B)

  print(f"R: {R} \n")
  print(f"t: {t} \n")

// output
>>> R
array([[ 1.0,  0.0],
       [0.0,  1.0]])
>>> t
array([ 0.0, 0.0])

Since I'm passing two identical point clouds to the best_fit_transform, It makes sense that the resulting rotation matrix is the identity matrix and the translation vector is zero. However, if I run this algorithm in rust with the nalgebra implementation of SVD.

FUNCTION

/// Calculates the least-squares best-fit transform that maps corresponding points pc_a to pc_b
/// in two-dimensions (x,y)
/// Input:
///     A: pointcloud in previous step
///     B: pointcloud in current step
/// Returns:
///     R: rotation matrix
///     t: translation vector
pub fn best_fit_transform(mut A: Vec<na::Point2<f64>>, mut B: Vec<na::Point2<f64>>) -> (na::Matrix2<f64>, na::Vector2<f64>) {

    // get centroid of each point cloud
    let centroid_A: na::Vector2<f64> = A.centroid();
    let centroid_B: na::Vector2<f64> = B.centroid();

    // compute cross-covariance
    let H: na::Matrix2<f64> = A
        .iter()
        .zip(B.iter())
        .map(|(a, b)| (a.clone().coords - centroid_A, b.clone().coords - centroid_B)) // center point clouds by subtracting centroid coordinate
        .map(|(aa, bb)| bb * aa.transpose())
        .fold(na::Matrix2::zeros(), |sum, m| sum + m);

    // compute SVD
    let svd = na::linalg::SVD::new(H, true, true);
    let u: na::Matrix2<f64> = svd.u.unwrap();
    let v_t: na::Matrix2<f64> = svd.v_t.unwrap();
    // let u = na::Matrix2::new(-0.60464909, -0.79649199,
    //                                     -0.79649199, 0.60464909);
    //
    // let v_t = na::Matrix2::new(-0.60464909, -0.79649199,
    //                                      -0.79649199, 0.60464909);

    // Check that decomposition produces orthogonal left and right singular vectors
    assert_eq!(u.is_orthogonal(0.01), true);
    assert_eq!(v_t.is_orthogonal(0.01), true);

    // get rotation matrix
    let mut R: na::Matrix2<f64> = v_t.transpose() * u.transpose();

    // special reflection case
    if R.determinant() < 0.0 {
        R.set_row(1, &(-1.0 * R.row(1)));
    }

    // compute translation offset
    let t: na::Vector2<f64> = centroid_B - R * centroid_A;

    return (R, t)
}

MAIN Now recreating the same test as I did in python

  let p0 = na::Point2::new(1.0, 2.0);
  let p1 = na::Point2::new(2.0, 2.0);
  let p2 = na::Point2::new(3.0, 3.0);
  let p3 = na::Point2::new(-1.0, -2.0);
  let A = vec![p0, p1, p2, p3];
  let B = A.clone();

  let (rot, t) = best_fit_transform(A, B);

  println!("R: {}", rot);
  println!("t: {}", t);
  // output
R: 
  ┌                                         ┐
  │ 0.26879896781394647  -0.963196301333304 │
  │  0.9631963013333045 0.26879896781394674 │
  └                                         ┘

t: 
  ┌                    ┐
  │  2.117996666899197 │
  │ -0.289994086434064 │
  └                    ┘

This result for rotation matrix and translation vector does not make sense for a best fit transform between two identical point clouds. If I hard-code the u and v_t matrix (see commented out code in the rust implementation) I end up with the same results as in Python:

R: 
  ┌                                       ┐
  │ 1.0000000121719883                  0 │
  │                  0 1.0000000121719883 │
  └                                       ┘

t: 
  ┌                            ┐
  │ -0.00000001521498527168319 │
  │ -0.00000001521498527168319 │
  └                            ┘

QUESTION

• Is there anyway I can get past this issue in nalgebra atm? I'd like to achieve the same behavior as in numpy.
• If not, do you plan on changing this behavior to be consistent with what you would expect from numpy?

Ps! thank you for your time :smile: And keep up the good work with this library! :100:

[Sollimann]

@lebensterben I changed the variable names to make it a bit clearer :)

Basically, this step here:

  # compute cross-covariance
  AA = A - centroid_A
  BB = B - centroid_B
  H = np.dot(AA.T, BB)

corresponds to:

  // compute cross-covariance
  let centroid_A: na::Vector2<f64> = A.centroid();
  let centroid_B: na::Vector2<f64> = B.centroid();
  let H: na::Matrix2<f64> = A
      .iter()
      .zip(B.iter())
      .map(|(a, b)| (a.clone().coords - centroid_A, b.clone().coords - centroid_B)) // center point clouds by subtracting centroid coordinate
      .map(|(aa, bb)| bb * aa.transpose())
      .fold(na::Matrix2::zeros(), |sum, m| sum + m);

with the same input above, they both produce the same matrix cross-covariance matrix, H:

H: 
  ┌             ┐
  │  8.75 10.75 │
  │ 10.75 14.75 │
  └             ┘

I stepped through both algorithms, and the only behavior that is different is the SVD decomposition. The decomposition produced by numpy produces an expected rotation and translation matrix in the end, while the same implementation with nalgebra produces the wrong rotation and translation matrices due to a different decomposition for SVD. If I hardcode the decomposition instead of using the nalgebra decomposition I get the same rotation and translation as in my python implementation.

[Sollimann]

@Andlon @lebensterben do you think the nalgebra SVD decomposition would produce the same decomposition as numpy if nalgebra sorted its singular value in descending order?

[Sollimann]

@lebensterben I'm now able to successfully compute the rotation matrix and translation matrix by not transposing the vectors and matrices like I do in Python.

So in Python I transpose the v_t and u matrices and then compute the product to get the rotation matrix. I also transpose the centroids which are 2-dim vectors

U, S, Vt = np.linalg.svd(H)
R = np.dot(Vt.T, U.T)

# translation
t = centroid_B.T - np.dot(R,centroid_A.T)

in rust I had to skip the transpose to achieve the same result

let u: na::Matrix2<f64> = svd.u.unwrap();
let v_t: na::Matrix2<f64> = svd.v_t.unwrap();
let mut R: na::Matrix2<f64> = v_t * u;

// compute translation offset
let t: na::Vector2<f64> = centroid_B - R * centroid_A;

Thank you for your time and patience :+1: I hope though you could have a look into the debug format of the matrices. I found it hard to debug with the print showing matrices in column-major format.

[Andlon]

@Sollimann: I'm very confused. I ran both your numpy code and Rust code (before your latest change to transposition), and they already return comparable results, i.e. here's the numpy output:

H:  [[ 8.75 10.75]
 [10.75 14.75]] 

R: [[1.00000000e+00 0.00000000e+00]
 [2.22044605e-16 1.00000000e+00]] 

t: [ 0.00000000e+00 -2.22044605e-16] 

and here's the output from the Rust version:

H: 
  ┌             ┐
  │  8.75 10.75 │
  │ 10.75 14.75 │
  └             ┘

R: 
  ┌                                                                       ┐
  │                 1.0000000000000004  0.0000000000000010547118733938987 │
  │ -0.0000000000000007771561172376096                                  1 │
  └                                                                       ┘

t: 
  ┌                                    ┐
  │ -0.0000000000000017763568394002505 │
  │  0.0000000000000008881784197001252 │
  └                                    ┘

Apart from some numerical round-off error, they are identical. So I'm left a little confused as to what the problem is? Also, note that in the Rust version you compute aa * bb.transpose() where as in Python you compute np.dot(AA.T, B). The result appears to be the same, however.

The internal ordering of the elements is an internal implementation detail that you rarely have to worry about. Therefore you should index entries as matrix[[i, j]] and not access .data directly for this purpose. I'll make an issue for the Debug representation, though I am not sure what the solution should be like (in the end the purpose of Debug is kind of to display the internal data representation of the object for debugging purposes. It's just a coincidence that the internal data representation "looks" like a matrix).

The SVD of a matrix A is defined by the matrices U, S, V such that A = U * S * V^T, where U and V are orthogonal matrices and S is a diagonal matrix with non-negative entries. nalgebra satisfies all these properties and therefore returns a correct SVD. Given that the singular vectors U and V^T are not unique, it cannot be a goal for nalgebra to return the same results as numpy, as this is just an internal implementation detail, that might even change with different versions of numpy.

Finally, although I don't have the time to dig into what you are doing, I gather that you're computing something akin to the "best affine transformation" and then extracting the rotation associated with that affine transformation. In that case, you want the polar decomposition of A:

A = R * P

where R is a rotation matrix and P is a (usually defined as positive definite) symmetric matrix. Let V S V^T be the Eigenvalue decomposition of P. Then

A = R V S V^T = U S V^T

is the SVD of A, and

R V = U       =>                R = U V^T

So if this is what you're doing, it looks slightly off to me, but that might also be because of different conventions.

[Sollimann]

@Andlon Thanks for the answer! Hmm, I'm using bb * aa.tranpose() in the examples above so I don't think that's the problem :thinking:

But this is my new code now:

#[allow(non_snake_case)]
pub fn best_fit_transform(A: PointCloud, B: PointCloud) -> (na::Matrix2<f64>, na::Vector2<f64>) {

    // make sure dimensions are the same
    assert_eq!(A.size(), B.size());

    // convert my type of point to nalgebra point
    let to_na_p2: fn(Point) -> na::Point2<f64> = |p: Point| na::Point2::new(p.x, p.y);
    let to_na_v2: fn(Point) -> na::Vector2<f64> = |p: Point| na::Vector2::new(p.x, p.y);

    // get centroid of each point cloud
    let centroid_A: na::Vector2<f64> = to_na_v2(A.centroid());
    let centroid_B: na::Vector2<f64> = to_na_v2(B.centroid());

    // convert to nalgebra library and center point clouds
    let A_: Vec<na::Point2<f64>> = A.iter().map(|p: &Point| to_na_p2(*p)).collect();
    let B_: Vec<na::Point2<f64>> = B.iter().map(|p: &Point| to_na_p2(*p)).collect();

    // compute cross-covariance
    let H: na::Matrix2<f64> = A_
        .iter()
        .zip(B_.iter())
        .map(|(a, b)| (a.clone().coords - centroid_A, b.clone().coords - centroid_B))
        .map(|(aa, bb)| bb * aa.transpose())
        .fold(na::Matrix2::zeros(), |sum, m| sum + m)
        .transpose();

    // compute SVD
    // https://medium.com/machine-learning-world/linear-algebra-points-matching-with-svd-in-3d-space-2553173e8fed
    let svd = na::linalg::SVD::new(H, true, true);
    let U: na::Matrix2<f64> = svd.u.unwrap();
    let Vt: na::Matrix2<f64> = svd.v_t.unwrap();

    // Check that decomposition produces orthogonal left and right singular vectors
    assert_eq!(U.is_orthogonal(0.01), true);
    assert_eq!(Vt.is_orthogonal(0.01), true);

    // get rotation matrix
    let mut R: na::Matrix2<f64> = Vt.transpose() * U.transpose();

    // special reflection case
    if R.determinant() < 0.0 {

        // create column vector from index m01 and m11
        let u: na::Vector2<f64>= get_col_vec(U.transpose(), 1);
        // create row vector from index m10 and m11
        let v: na::Vector2<f64> = get_row_vec(Vt.transpose(), 1)*2.0;

        // outer product
        let m = u * v.transpose();
        R -= m;
        assert!(R.determinant() >= 0.0)
    }

    // compute translation offset
    let t: na::Vector2<f64> = centroid_B - R * centroid_A;

    return (R, t)
}

I have debugged a lot, and comparing with the python code I am now able to compute the right rotation matrix R and translation vector t in all cases, expect for the cases where the R.determinant() < 0.0:

the python code, which btw I copyed from here: link solves the case of improper rotation matrix by doing:

m = A.shape[1]
# special reflection case
if np.linalg.det(R) < 0:
   Vt[m-1,:] *= -1
   R = np.dot(Vt.T, U.T)

which I created an equivalent rust code:

if R.determinant() < 0.0 {
  Vt.set_row(1, &(-1.0 * Vt.row(1)));
  R = Vt.transpose() * U.transpose();
}

which does not produce the same rotation matrix and translation vector. so I tried another approach which I found here:

m = A.shape[1]
if np.linalg.det(R) < 0.0:
    # R does not constitute right handed system
    R -= np.outer(U[:, m-1], Vt[m-1, :]*2.0)

My rust implementation of this was:

// special reflection case
if R.determinant() < 0.0 {

    // create column vector from index m01 and m11
    let u: na::Vector2<f64>= get_col_vec(U.transpose(), 1);
    // create row vector from index m10 and m11
    let v: na::Vector2<f64> = get_row_vec(Vt.transpose(), 1)*2.0;

    // outer product
    let m = u * v.transpose();
    R -= m;
}

Neither this code snippet for handling improper rotation matrix produces the same 'corrected' rotation matrix and translation vector.

question @lebensterben @Andlon

Is there always a unique rotation matrix and translation between two vectors, or could another combination of rotation and translation produce the same transformation between the two vectors? I curious if the rotation matrices and translation vector I get in rust with nalgebra could still be valid even though it's not the same as I get in numpy

Do you know of any other ways of handling improper rotation matrices (det(R) is negative) ?

[Andlon]

First, the way you calculate the cross-variance matrix is super weird.

@lebensterben Can I please ask you to refrain from such harsh language? Here we have a user who has problems reproducing his Python results in nalgebra. Surely we want to improve the user experience when porting code to nalgebra from other languages, and experience reports like these represent invaluable data that we may use to improve both documentation and the codebase. I think calling users' attempts "super weird" might alienate them from persisting in their efforts - and moreover I don't think it's a nice way of communication in the first place. The Rust community has an image of being open and welcoming to new users, and I think we very much should extend this to nalgebra. With that in mind, it is important to look at the interaction in places like GitHub and Discord and make sure that we use a friendly tone that encourages users to stick around. So, perhaps instead of claiming that the way they calculate something is "super weird", you might instead propose a simpler way of thinking about it, or a more "standard" way to approach the problem, if you think this helps?

For the record, I do not think the way the cross-covariance matrix is computed is "super weird" at all. There are many lenses through which you can view this particular problem. Just because the way @Sollimann has approached the problem is different from the way you might look at it, please consider that in mathematics there are often many ways to reach the same goal. What is clear to one person may be foreign to another.

@Sollimann: regarding my comment on transposition in the computation of H. I missed that you have a .transpose() at the end of the chain of methods. With the "reversed" order that I pointed out earlier this effectively transposes twice, so you end up with the same as in Python.

I'd be happy to help you further investigate this, but my time is very limited. Since you say now that it seems to work for proper rotations, can I suggest that you split out a function that only fixes the improper rotation and then provide input to this function that gives an unexpected result? That way the problem is much more isolated and I think it would be easier to help you. Flipping the sign of any row in v_t should fix the problem, as far as I can tell, since you are anyway not using the singular values for anything. Here are two alternative approaches for flipping the sign of a row in nalgebra:

{
    // Using the *= operator unfortunately does not permit an expression on the left-hand side, so we need to create
    // an auxiliary variable
    let mut row = matrix.row_mut(1);
    row *= -1.0;
}

{
    // But we can call `mul_assign` directly instead of using the *= operator
    matrix.row_mut(1).mul_assign(-1.0);
}

Vt[m-1,:] *= -1 it's changing the sign of m-1 row of Vt, which is the m-1 column of V.

This relies on the invariant that result of SVD is ordered by the magnitude of singular values.

So you need to find the m-1 largest singular value, or second smallest singular value. The change the corresponding line in Vt.

I don't think fixing the rotation here relies on any particular ordering of the singular values. It's conventional to "fix" rotations by flipping the sign of the row/col corresponding to the smallest singular value, but this is not a requirement. Flipping the sign of any column in V, for example, will flip the sign of det(R).

EDIT: I should say that flipping the sign corresponding to the smallest singular value may be a good idea for purposes of reducing numerical error, but I think that is not necessary to resolve the current issue at hand.

[Sollimann]

@Andlon @lebensterben Thank you for your patience! :) So first and foremost, I'd like to give a big huge thanks to both of you for even bothering trying to investigate this with me. I do not expect anyone to take all this time helping me out. I had a look at the ndarray library which is more similar to numpy, whereas I find nalgebra more similar to Eigen in C++. However, I found the documentation and usability of nalgebra to be much better than ndarray, so I wanted to give a shot on implementing this in your library :)

So I have set up two scenarios where I have two identical functions and two identical datasets.

functions

numpy

def handle_improper_rotation(R, U, Vt):

    print(f"R det before: {np.linalg.det(R)} \n \n")
    print(f"R before: {R} \n \n")
    print(f"U: {U} \n \n")
    print(f"Vt before: {V} \n \n")

    m = R.shape[0] # (2,2)
    Vt[m-1,:] *= -1
    R = np.dot(Vt.T, U.T)

    print(f"R det before: {np.linalg.det(R)} \n \n")
    print(f"Vt after: {Vt} \n")
    print(f"R after: {R} \n \n")
    return R

nalgebra

type M2x2 = na::Matrix2<f64>;
#[allow(non_snake_case)]
pub fn handle_improper_rotation(mut R: M2x2, U: M2x2, mut Vt: M2x2) -> M2x2 {

    println!("R det before: {}", R.determinant());
    println!("R before: {}", R);
    println!("U: {}", U);
    println!("Vt before: {}", Vt);

    Vt.row_mut(1).mul_assign(-1.0);
    R = Vt.transpose() * U.transpose();
    assert!(R.determinant() >= 0.0);

    println!("Vt after: {}", Vt);
    println!("R det after: {}", R.determinant());
    println!("R after: {}", R);

    return R
}

Test case 1 (Success :heavy_check_mark: )

Test setup

numpy

U = np.array([[-0.0099995, -0.99995],
              [-0.99995, 0.0099995]])

Vt = np.array([[0.0099995, 0.99995],
               [-0.99995,0.0099995]])

R = np.dot(Vt.T, U.T)

if np.linalg.det(R) < 0.0:
    print("RESULT: \n \n")
    R = handle_improper_rotation(R, U, Vt)

nalgebra

#[test]
#[allow(non_snake_case)]
fn improper_rotation_case_1()  {
    let U = na::Matrix2::new(-0.0099995, -0.99995,
                             -0.99995, 0.0099995);

    let Vt = na::Matrix2::new( 0.0099995, 0.99995,
                              -0.99995,0.0099995);

    let mut R = Vt.transpose() * U.transpose();

    if R.determinant() < 0.0 {
        println!("RESULT: \n \n ");
        R = handle_improper_rotation(R, U, Vt);
    }
}

Test Result ( equal results in both numpy and nalgebra :heavy_check_mark: )

numpy

RESULT: 

R det before: -0.9999999850005001 

R before: [[ 0.99980001 -0.019998  ]
 [-0.019998   -0.99980001]] 

U: [[-0.0099995 -0.99995  ]
 [-0.99995    0.0099995]] 

Vt before: [[ 0.0099995  0.99995  ]
 [-0.99995    0.0099995]] 

Vt after: [[ 0.0099995  0.99995  ]
 [ 0.99995   -0.0099995]] 

R det after: 0.9999999850005 

R after: [[-9.99999993e-01  2.21337554e-19]
 [ 2.21337554e-19 -9.99999993e-01]] 

nalgebra

RESULT: 

R det before: -0.9999999850005001
R before: 
  ┌                                     ┐
  │  0.99980001249975    -0.01999800005 │
  │    -0.01999800005 -0.99980001249975 │
  └                                     ┘

U: 
  ┌                       ┐
  │ -0.0099995   -0.99995 │
  │   -0.99995  0.0099995 │
  └                       ┘

Vt before: 
  ┌                     ┐
  │ 0.0099995   0.99995 │
  │  -0.99995 0.0099995 │
  └                     ┘

Vt after: 
  ┌                       ┐
  │  0.0099995    0.99995 │
  │    0.99995 -0.0099995 │
  └                       ┘

R det after: 0.9999999850005
R after: 
  ┌                                     ┐
  │ -0.99999999250025                 0 │
  │                 0 -0.99999999250025 │
  └                                     ┘

Test case 2 (Success :heavy_check_mark: )

Test setup

numpy

U = np.array([[-0.9998181379392826, 0.01907068555730712],
              [-0.01907068555730701, -0.9998181379392826]])

Vt = np.array([[0.0007255423829884326, 0.9999997367940907],
               [0.9999997367940907, -0.0007255423829884326]])

R = np.dot(Vt.T, U.T)

if np.linalg.det(R) < 0.0:
    print("RESULT: \n \n")
    R = handle_improper_rotation(R, U, Vt)

nalgebra

#[test]
#[allow(non_snake_case)]
fn improper_rotation_case_2()  {
    let U = na::Matrix2::new(-0.9998181379392826, 0.01907068555730712,
                              -0.01907068555730701, -0.9998181379392826);

    let Vt = na::Matrix2::new( 0.0007255423829884326, 0.9999997367940907,
                               0.9999997367940907, -0.0007255423829884326);

    let mut R = Vt.transpose() * U.transpose();

    if R.determinant() < 0.0 {
        println!("RESULT: \n \n ");
        R = handle_improper_rotation(R, U, Vt);
    }
}

Test Result ( equal results in both numpy and nalgebra :heavy_check_mark: )

numpy

RESULT: 

R det before: -1.0000000000000002 

R before: [[ 0.01834527 -0.99983171]
 [-0.99983171 -0.01834527]] 

U: [[-0.99981814  0.01907069]
 [-0.01907069 -0.99981814]] 

Vt before: [[ 0.0099995  0.99995  ]
 [-0.99995    0.0099995]] 

Vt after: [[ 7.25542383e-04  9.99999737e-01]
 [-9.99999737e-01  7.25542383e-04]] 

R det after: 1.0000000000000002 

R after: [[-0.01979609  0.99980404]
 [-0.99980404 -0.01979609]] 

nalgebra

RESULT: 

R det before: -1.0000000000000002
R before: 
  ┌                                             ┐
  │  0.018345270103434466    -0.999831711371885 │
  │    -0.999831711371885 -0.018345270103434355 │
  └                                             ┘

U: 
  ┌                                           ┐
  │  -0.9998181379392826  0.01907068555730712 │
  │ -0.01907068555730701  -0.9998181379392826 │
  └                                           ┘

Vt before: 
  ┌                                               ┐
  │  0.0007255423829884326     0.9999997367940907 │
  │     0.9999997367940907 -0.0007255423829884326 │
  └                                               ┘

Vt after: 
  ┌                                             ┐
  │ 0.0007255423829884326    0.9999997367940907 │
  │   -0.9999997367940907 0.0007255423829884326 │
  └                                             ┘

R det after: 1.0000000000000004
R after: 
  ┌                                             ┐
  │ -0.019796090972145512     0.999804038190596 │
  │    -0.999804038190596   -0.0197960909721454 │
  └                                             ┘

My conclusion from these tests

So from testing, it seems to me that given the same input U and Vt, the two methods in numpy and nalgebra produce the same output R. The problem I have occurs due to the fact that when I have the cross-covariance matrix H, and I feed that to the SVD decomposition function in numpy and in nalgebra they produce the different decompositions of U and Vt. With the two different decompositions I am still able to produce the same R matrix for both numpy and nalgebra, but that is only true when R.determinant() >= 0.0. In the case that R.determinant() < 0.0 and I have different decompositions of U and Vt, I end up with different results in many cases (not all)

My methods

So, I was thinking. Instead of pasting all my tests here, which I think could be very long and mabye also confusing. Mabye you would like the run the tests yourself? :) Here is a link to all my tests on best fit transform: tests . I have ran the same tests in Python, so the values to I expect and that I compare with in end are my results from numpy (which I consider correct and ground truth). The methods I use for my tests are found here . Ps! these are not found on the master branch, but particle-filter branch

In case you don't want to run the tests or have a look at my methods. Here are the most important once that I use:

#[derive(Debug, Copy, Clone)]
pub struct Point {
    pub x: f64,
    pub y: f64
}

#[derive(Debug, Clone)]
pub struct PointCloud {
    points: Vec<Point>,
}

impl PointCloud {
    pub fn size(&self) -> usize {
        self.points.len()
    }
    pub fn centroid(&self) -> Point {
        let x_avg: Scalar = self.iter().map(|p: &Point| p.x).sum::<f64>() / self.size();
        let y_avg: Scalar = self.iter().map(|p: &Point| p.y).sum::<f64>() / self.size();
        Point::new(x_avg, y_avg)
    }
}

#[allow(non_snake_case)]
pub fn best_fit_transform(A: PointCloud, B: PointCloud) -> (na::Matrix2<f64>, na::Vector2<f64>) {

    // make sure dimensions are the same
    assert_eq!(A.size(), B.size());

    // convert my type of point to nalgebra point
    let to_na_p2: fn(Point) -> na::Point2<f64> = |p: Point| na::Point2::new(p.x, p.y);
    let to_na_v2: fn(Point) -> na::Vector2<f64> = |p: Point| na::Vector2::new(p.x, p.y);

    // get centroid of each point cloud
    let centroid_A: na::Vector2<f64> = to_na_v2(A.centroid()); // A.centroid computes the mean point 
    let centroid_B: na::Vector2<f64> = to_na_v2(B.centroid());

    // convert to nalgebra library and center point clouds
    let A_: Vec<na::Point2<f64>> = A.iter().map(|p: &Point| to_na_p2(*p)).collect();
    let B_: Vec<na::Point2<f64>> = B.iter().map(|p: &Point| to_na_p2(*p)).collect();

    // compute cross-covariance
    let H: na::Matrix2<f64> = A_
        .iter()
        .zip(B_.iter())
        .map(|(a, b)| (a.clone().coords - centroid_A, b.clone().coords - centroid_B))
        .map(|(aa, bb)| bb * aa.transpose())
        .fold(na::Matrix2::zeros(), |sum, m| sum + m)
        .transpose();

    // compute SVD
    // https://medium.com/machine-learning-world/linear-algebra-points-matching-with-svd-in-3d-space-2553173e8fed
    let svd = na::linalg::SVD::new(H, true, true);
    let U: na::Matrix2<f64> = svd.u.unwrap();
    let mut Vt: na::Matrix2<f64> = svd.v_t.unwrap();

    // Check that decomposition produces orthogonal left and right singular vectors
    assert_eq!(U.is_orthogonal(0.01), true);
    assert_eq!(Vt.is_orthogonal(0.01), true);

    // get rotation matrix
    let mut R: na::Matrix2<f64> = Vt.transpose() * U.transpose();

    // special reflection case
    if R.determinant() < 0.0 {
        R = handle_improper_rotation(R, U, Vt);
    }

    // compute translation offset
    let t: na::Vector2<f64> = centroid_B - R * centroid_A;

    return (R, t)
}

type M2x2 = na::Matrix2<f64>;
#[allow(non_snake_case)]
pub fn handle_improper_rotation(mut R: M2x2, U: M2x2, mut Vt: M2x2) -> M2x2 {

    // println!("R det before: {}", R.determinant());
    // println!("R before: {}", R);
    // println!("U: {}", U);
    // println!("Vt before: {}", Vt);

    Vt.row_mut(1).mul_assign(-1.0);
    R = Vt.transpose() * U.transpose();
    assert!(R.determinant() >= 0.0);

    // println!("Vt after: {}", Vt);
    // println!("R det after: {}", R.determinant());
    // println!("R after: {}", R);

    return R
}

An important observation!

So something I have noticed in my tests, is that if any of the tests (that I linked to here: tests) fails, then I can just switch this this row here Vt.row_mut(1).mul_assign(-1.0) to be Vt.row_mut(0).mul_assign(-1.0) instead in function handle_improper_rotation. The test that first failed, will now succeed with the correct rotation matrices. However, the tests that were originally succeeding are now failing again :thinking: When running the same examples in numpy I don't get this behavior.

[Sollimann]

if you need me to provide more test cases or more context to some on any of the things above, please let me know! :)

[Andlon]

Ah, I see that I misunderstood you earlier, and I hadn't thought it through what it is you're doing. So, when "fixing" the rotation matrix in order to obtain a proper rotation, there are indeed different choices: you can flip any of the columns in V, or any of the columns in U. Each of these choices lead to a different rotation matrix.

The implementation you presented using numpy will indeed, as @lebensterben pointed out, always decide to flip the sign of the column in V corresponding to the smallest singular value, because numpy returns singular values in descending order. Since nalgebra (currently) does not impose any ordering on its singular values, you may or may not get the same result as numpy. If you want the same result as numpy, you might want to first find the index of the singular value that is smallest, and then flip the corresponding row in v_t. This should give you consistently the same results as with numpy.

To summarize my understanding of what is happening here:

1. The algorithm first finds a "best fit" affine transformation
2. This "best fit" affine transformation may correspond to a reflection.
3. By (conditionally) flipping a row in v_t, you ensure that you still get a proper rotation. However, if you were to multiply your matrices U S V^T, you will no longer have the same affine transformation.
4. I think that by modifying the right singular vector (column in V) that corresponds to the smallest singular value, you might perhaps minimize the deviation from the original "best fit" transformation in some sense. I am unsure of what mathematical guarantees you get here, if any.

[Sollimann]

@Andlon The ordering seem to have solved my issue :) All tests are passing now. Thanks for the help! I'll close this now

[Andlon]

@Sollimann: I created the issues #897 and #898 in an attempt to improve the situation around eigen/singular value ordering and confusing debug output of matrices. Thank you for your persistence in this matter; the issue has revealed some pain points in nalgebra that we can hopefully improve in the future to make the transition for new users more painless!

[ChristopherRabotin]

Hi there,

I'm also encountering differences in the SVD decomposition between nalgebra 0.32 and numpy (version 2). In short, with a 9x9 covariance matrix, numpy returns a distribution matrix where only the last three rows and columns are zeros, but nalgebra returns zeros on the last five rows and three columns. I wonder if the source of the issue isn't me misunderstanding the mathematics of a singular value decomposition.

In my specific use case, I'm trying to implement a multivariate normal distribution, using numpy's multivariate_normal function.

I've implemented it as follows. There are a few transpose calls in there because that's the closest I've been to numpy's results.

        let svd = cov.svd(false, true); // Compute u, but not v, and get the ordered SVDs.
        if svd.v_t.is_none() {
            return Err(Box::new(NyxError::CovarianceMatrixNotPsd));
        }

        let sqrt_s = svd.singular_values.map(|x| x.sqrt());
        let mut sqrt_s_v_t = svd.v_t.unwrap().transpose();

        for (i, mut col) in sqrt_s_v_t.column_iter_mut().enumerate() {
            col *= -sqrt_s[i]; // np.sqrt(s) is an item-by-item square root
        }

        println!("{:.6e}", sqrt_s_v_t.transpose());

I have a 9x9 covariance matrix as follows:

  ┌                                                                                                                               ┐
  │   0.024406137   0.412504544  -1.811501993   0.000031159   0.000526633  -0.002312628   0.000000000   0.000000000   0.000000000 │
  │   0.412504544   6.972016745 -30.617413550   0.000526621   0.008900761  -0.039087329   0.000000000   0.000000000   0.000000000 │
  │  -1.811501993 -30.617413550 134.455502158  -0.002312632  -0.039087329   0.171650920   0.000000000   0.000000000   0.000000000 │
  │   0.000031159   0.000526621  -0.002312632   0.000000266   0.000003145  -0.000002399   0.000000000   0.000000000   0.000000000 │
  │   0.000526633   0.008900761  -0.039087329   0.000003145   0.000038338  -0.000043939   0.000000000   0.000000000   0.000000000 │
  │  -0.002312628  -0.039087329   0.171650920  -0.000002399  -0.000043939   0.000221444   0.000000000   0.000000000   0.000000000 │
  │   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000 │
  │   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000 │
  │   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000   0.000000000 │
  └                                                                                                                               ┘

The output of sqrt_s_v_t.transpose() is the following:

┌                                                                                                                               ┐
  │  -1.562246e-1   -2.640458e0    1.159549e1  -1.994423e-4  -3.370907e-3   1.480324e-2   -0.000000e0   -0.000000e0   -0.000000e0 │
  │   2.676776e-6   4.625904e-6   1.089447e-6   4.759677e-4   5.189239e-3   1.188076e-3   -0.000000e0   -0.000000e0   -0.000000e0 │
  │   1.628157e-8   2.752071e-7  -1.208560e-6  -1.275238e-5  -2.155714e-4   9.466745e-4   -0.000000e0   -0.000000e0   -0.000000e0 │
  │   -0.000000e0 -5.186472e-13 -1.181092e-13  -2.102264e-9  1.894320e-10  1.481738e-11   -0.000000e0   -0.000000e0   -0.000000e0 │
  │    0.000000e0   -0.000000e0   -0.000000e0    0.000000e0    0.000000e0    0.000000e0   -0.000000e0   -0.000000e0   -0.000000e0 │
  │   -0.000000e0   -0.000000e0   -0.000000e0    0.000000e0    0.000000e0    0.000000e0   -0.000000e0   -0.000000e0   -0.000000e0 │
  │   -0.000000e0   -0.000000e0    0.000000e0   -0.000000e0   -0.000000e0    0.000000e0    0.000000e0   -0.000000e0   -0.000000e0 │
  │   -0.000000e0   -0.000000e0    0.000000e0   -0.000000e0   -0.000000e0    0.000000e0   -0.000000e0    0.000000e0   -0.000000e0 │
  │   -0.000000e0   -0.000000e0    0.000000e0   -0.000000e0   -0.000000e0    0.000000e0   -0.000000e0   -0.000000e0    0.000000e0 │
  └                                                                                                                               ┘

When I plug-in that same matrix in numpy, and run through the instructions manually, the sqrt_s_v result only matches nalgebra on the first row:

In [34]: cov = np.array([[0.024406, 0.412505, -1.811502, 0.000031, 0.000527, -0.002313, 0.000000, 0.000000, 0.000000, ],[0.412505, 6.972017, -30.617414, 0.000527, 0.008901, -0.039087, 0.000000, 0.000000, 0.000000, ],[-1.8115
    ...: 02, -30.617414, 134.455502, -0.002313, -0.039087, 0.171651, 0.000000, 0.000000, 0.000000, ],[0.000031, 0.000527, -0.002313, 0.000000, 0.000003, -0.000002, 0.000000, 0.000000, 0.000000, ],[0.000527, 0.008901, -0.0390
    ...: 87, 0.000003, 0.000038, -0.000044, 0.000000, 0.000000, 0.000000, ],[-0.002313, -0.039087, 0.171651, -0.000002, -0.000044, 0.000221, 0.000000, 0.000000, 0.000000, ],[0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
    ...: 0.000000, 0.000000, 0.000000, 0.000000, ],[0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, ],[0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
    ...:  0.000000, ],])

In [35]: (u, s, v) = svd(cov)

In [36]: sqrt_s_v = (np.sqrt(s)[:, None] * v)

In [37]: sqrt_s_v[0]
Out[37]:
array([-1.56224647e-01, -2.64045772e+00,  1.15954949e+01, -1.99479347e-04,
       -3.37088489e-03,  1.48032430e-02, -0.00000000e+00, -0.00000000e+00,
       -0.00000000e+00])

In [38]: sqrt_s_v[1]
Out[38]:
array([-5.33121957e-05, -7.83152256e-05, -1.85619949e-05, -4.65192515e-04,
       -5.15722661e-03, -1.17260911e-03, -0.00000000e+00, -0.00000000e+00,
       -0.00000000e+00])

In [39]: sqrt_s_v[2]
Out[39]:
array([ 3.24550472e-04, -1.81497622e-04, -3.58872037e-05, -3.69571736e-04,
        2.13468887e-04, -7.94303697e-04, -0.00000000e+00, -0.00000000e+00,
       -0.00000000e+00])

In [40]: sqrt_s_v
Out[40]:
array([[-1.56224647e-01, -2.64045772e+00,  1.15954949e+01,
        -1.99479347e-04, -3.37088489e-03,  1.48032430e-02,
        -0.00000000e+00, -0.00000000e+00, -0.00000000e+00],
       [-5.33121957e-05, -7.83152256e-05, -1.85619949e-05,
        -4.65192515e-04, -5.15722661e-03, -1.17260911e-03,
        -0.00000000e+00, -0.00000000e+00, -0.00000000e+00],
       [ 3.24550472e-04, -1.81497622e-04, -3.58872037e-05,
        -3.69571736e-04,  2.13468887e-04, -7.94303697e-04,
        -0.00000000e+00, -0.00000000e+00, -0.00000000e+00],
       [-6.02972796e-04,  4.89822403e-04,  1.03687430e-04,
        -3.07427774e-04,  7.17752362e-05, -2.00652929e-04,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
       [ 1.55549924e-04, -2.04330859e-05, -2.98084653e-06,
        -5.52415966e-04, -2.39768318e-05,  3.18943941e-04,
        -0.00000000e+00, -0.00000000e+00, -0.00000000e+00],
       [ 3.73260269e-04,  4.86667635e-04,  1.15855589e-04,
         8.30414366e-05, -1.75042727e-05, -7.26591222e-06,
        -0.00000000e+00, -0.00000000e+00, -0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00]])

The SVD decomposition also does not match numpy's result if I grab a fixed_view of my covariance to restrict it to the non-zero upper 6x6.

Thanks for any help on this.