Null Space Calculation using Eigen is wrong?

Question

I'm using the Eigen library (http://eigen.tuxfamily.org) to do a Null Space calculation using the SVD function. I compared the output to the "Null" function in matlab, and got different results. Stepping through it with a debugger, and looking at the V matrix created by Eigen vs. the V matrix from matlab, there's an odd difference.

The left-singular vectors (left 3 columns in the example below) in the V matrix are almost the same, but the signs are switched. The right singular vectors (the null space; right 3 columns below) are not very similar at all.

Any idea what would cause this? Am I using the SVD Function incorrectly? Code and example results below.

Here's the code: "input" is a normal C++ array:

/* Create a matrix with the nessecary size */
MatrixXf A(inRows, inCols);

/* Populate the matrix from the input */
for (int i=0; i < inRows; i++)
{
  for(int j=0; j < inCols; j++)
  {
    A(i,j) = input[i*inCols + j];
  }
}

/* Do a singular value decomposition on the matrix */
JacobiSVD<MatrixXf> svd(A, Eigen::ComputeFullV);

/* Get the V matrix */
MatrixXf V((int)svd.matrixV().rows(), (int)svd.matrixV().cols());
V = svd.matrixV();

Here are some example results:

A(input) =

-0.5059    -0.0075   -0.0121   -0.3526   -0.3528   -0.0128
-0.0067     0.4915    0.0235   -0.3503    0.3559    0.0211
 0.0027     0.0010   -0.5015    0.0021   -0.0031   -0.4999

V(Matlab) =

 0.3120    0.6304    0.1115   -0.5031   -0.4895   -0.0027
 0.3628   -0.2761    0.5333    0.4955   -0.5121   -0.0018
 0.5180   -0.1804   -0.4480   -0.0002    0.0000   -0.7060
-0.0353    0.6404   -0.2953    0.7081    0.0074   -0.0023
 0.4859    0.2283    0.4623    0.0032    0.7057    0.0048
 0.5151   -0.1775   -0.4489    0.0014   -0.0080    0.7082

V(Eigen) =

-0.3120   -0.6304   -0.1115   -0.5040   -0.4886   -0.0038
-0.3628    0.2761   -0.5333    0.4638   -0.4832    0.2432
-0.5180    0.1804    0.4480    0.1693   -0.1736   -0.6630
 0.0353   -0.6404    0.2953    0.6878    0.0257    0.1666
-0.4859   -0.2283   -0.4623    0.0258    0.6851   -0.1677
-0.5151    0.1775    0.4489   -0.1689    0.1665    0.6674

Thank you for your assistance!

Answer 1

Firstly, there is more than one way to form an orthonormal basis for a space. (For example [1 0; 0 1] and 1/sqrt(2) * [ 1 -1; 1 1 ] both describe the same 2D Euclidean space). So we wouldn't necessarily expect two alternative implementations to pick the same basis set.

If we take the right-hand three columns in each case, we learn the following:

> Vmat = Vmat(:,4:end);
> Veig = Veig(:,4:end);
> Vmat' * V_mat

ans =

 1.0000e+00   8.8800e-06  -1.4120e-05
 8.8800e-06   9.9999e-01  -5.1830e-05
-1.4120e-05  -5.1830e-05   1.0000e+00

> Veig' * Veig

ans = 

 1.0001e+00  -1.4050e-05   2.4200e-06
-1.4050e-05   1.0001e+00  -4.8310e-05
 2.4200e-06  -4.8310e-05   1.0000e+00

> A * Vmat

ans =

 7.7612e-17   7.8916e-17   0.0000e+00
-4.1193e-17   4.8139e-17   0.0000e+00
 6.6136e-18  -6.0715e-18   1.1102e-16

> A * Veig

ans = 

-1.2030e-05   1.1000e-05  -6.0000e-07
-4.8600e-06   3.8750e-05   1.5490e-05
-3.4400e-06  -4.5210e-05  -3.6090e-05

So these are both orthonormal basis sets, and they're both basically null spaces. However, the error level in the Eigen case appears to correspond to the fact that it was done in single-precision. Try it again in double-precision, and see how the results compare this time (I'm not claiming this will definitely help, merely that this is one obvious difference with Matlab.)