Hello, I'm trying to use the LeastSquare to do a linear regression, but it gives an innacurate (although close) answer.
The code I'm using is the exemple code on the LeastSquare module page :

http://eigen.tuxfamily.org/dox-2.0/grou ... odule.html

Code: Select all

#include <iostream>
#include <vector>
#include <../Eigen/LeastSquares>

using namespace Eigen;
int main()
{

    Vector3d points[5];
    points[0] = Vector3d( 3.02, 6.89, -4.32 );
    points[1] = Vector3d( 2.01, 5.39, -3.79 );
    points[2] = Vector3d( 2.41, 6.01, -4.01 );
    points[3] = Vector3d( 2.09, 5.55, -3.86 );
    points[4] = Vector3d( 2.58, 6.32, -4.10 );

    // create a vector of pointers to the points
    std::vector<Vector3d*> points_ptrs(5);
    for(int k=0; k<5; ++k) points_ptrs[k] = &points[k];
    Vector3d coeffs; // will store the coefficients a, b, c
    linearRegression(
                5,
                &(points_ptrs[0]),
            &coeffs,
            1 // the coord to express as a function of
            // the other ones. 0 means x, 1 means y, 2 means z.
            );

    std::cout << coeffs << std::endl;
    for(int k = 0; k<5; k++)
        std::cout << "y_incorrect(" << k+1 << ") = " << coeffs[0]*points[k][0]+coeffs[1]*points[k][2]+coeffs[2] << std::endl;

    Vector3d correct_coeffs = Vector3d(0.495, -1.927, -2.906);
    std::cout << std::endl << correct_coeffs << std::endl;
    for(int k = 0; k<5; k++)
        std::cout << " y_correct(" << k+1 << ") = " << correct_coeffs[0]*points[k][0]+correct_coeffs[1]*points[k][2]+correct_coeffs[2] << std::endl;
}


I just changed one detail in the declaration of points_ptr :
It is an array of pointer, so I think it should be std::vector<Vector3d*> and not std::vector<Vector3d>, like in the doc page.

This gives me the following results :

Code: Select all

-1.693
-6.189
-14.72
y_incorrect(1) = 6.901;
y_incorrect(2) = 5.331;
y_incorrect(3) = 6.015;
y_incorrect(4) = 5.628;
y_incorrect(5) = 6.284;

0.495
-1.927
-2.906
 y_correct(1) = 6.914;
 y_correct(2) = 5.392;
 y_correct(3) = 6.014;
 y_correct(4) = 5.567;
 y_correct(5) = 6.272;

And the correct results should be (0.495, -1.927, -2.906), just as the page said, as I also checked on Matlab, and lsqlin gives me the same result.

Any idea what could be wrong ?
(I tried in Eigen2 and Eigen3, so I dont think it's a version incompatibility problem).
Thank you in advance.

You may be aware of this, but the documentation page you cite pertains to Eigen 2.0. The current version is Eigen 3.2; the last release in the Eigen 2.0 series was done more than two years ago. It would probably be best for you to upgrade as I think it is unlikely that there will be any further releases.

The LeastSquares module was removed in Eigen 3; its C-style API was not consistent with the rest of the library. Have a look at http://eigen.tuxfamily.org/dox/group__T ... tml#title4 for an explanation on how to do least squares using the SVD decomposition.

Thank you, I just saw that.
I directly made a search for "least square" in the API, and that was the only result that came, I knew it was from 2.0, but I didn't know it was that obsolete.
I tested the SV Decomposition, it works. Thank you.

However, if it is an obsolete module, why keep it in Eigen3 ?

This is a simple self contained example how to do LS using the QR decomposition using the latest Eigen 3.2 version:

Code: Select all

VectorXd least_squares_solve(const MatrixXd& A, const VectorXd& b) {
   Eigen::ColPivHouseholderQR<MatrixXd> QR = A.colPivHouseholderQr();
   VectorXd x = QR.solve(b);
   return x;
}

This works because you try to solve:
x = (A^T*A)^-1*A^t*b or (A^T*A)*x = A^t*b
if you replace A by A=QR then you get
(R^T*Q^T*Q*R)*x=R^T*Q^T*b => simplifies to R*x=Q^T*b right solve which is exactly what the code above does. Note how here the solution is numerically very stable since we don't have to do A^T*A nor invert it.