SVD fails with non-square matrices (Bug #1370)

I noticed some weird behaviour when using non-square matrices as input to OpenCV's singular value decomposition (cv::SVD). Note that FULL_UV has to be computed.

An example symbolises the problem. Let's have P = [1, 0, 0]^T^. This can be decomposed as P = U D V^T^ with U = I, and D is [1, 0, 0]^T^ and V = r1.

As one of many tests that a SVD is correct we have U^T^ P = D V^T^. Why this is so, can be seen by having P = UDV^T and prepending U^T^ on both sides ( and knowing U^T^ U = I).

In the case where two eigenvalues are zero, the last two rows of D V^T^ must be zero.

So, what does cv::SVD do? It finds the eigenvalues correctly, V is also correct, but U is wrong.

[1, -0.70710677, -0.70710677;
  0, -0.70710677, 0.29325801;
  0, 0, 0.29325682]

Hence, U^T^ P is [1; -0.70710677; -0.70710677].

This problem occurs not only with this matrix but also with 6x4 matrices and the matrix [1, 2, 3]^T^. I postulate that it happens with all non-square matrices.

I looked at Eigen and they get the correct result (just as commercial software Mathematica and Matlab). However, they changed their SVD implementation to JacobiSVD which "[ensures] optimal reliability and accuracy" r2

r1
r2 http://eigen.tuxfamily.org/api/classEigen_1_1JacobiSVD.html

Minimal Example:

#include <iostream>

#include <opencv/cv.h>
#include <opencv/highgui.h>
using cv::Mat;
using cv::Matx;

#include <Eigen/Core>
#include <Eigen/SVD>

int main()
{
    // The results are not valid using [[OpenCV]]
    Matx<float, 3, 1> A(1,0,0);
    cv::SVD svdA(A, cv::SVD::FULL_UV);
    std::cout << svdA.u.t() * Mat(A) << std::endl;
    std::cout << svdA.u << std::endl;
    std::cout << svdA.w << std::endl;
    std::cout << svdA.vt << std::endl << std::endl;

    // Showing that this works in Eigen3 NOTE: Eigen3 required!
    Eigen::Matrix<float, 3, 1> B(1,0,0);
    Eigen::JacobiSVD<Eigen::MatrixXf> svdEigen(B, Eigen::ComputeFullU | Eigen::ComputeFullV);
    std::cout << svdEigen.matrixU().transpose() * B << std::endl;
    std::cout << svdEigen.matrixU() << std::endl;
    std::cout << svdEigen.matrixV() << std::endl;
    std::cout << svdEigen.singularValues() << std::endl;

    return 0;
}