[Python] C++ equivalent of Python's numpy.linalg.lstsq

I have a Python program that solves a equation system using the least squares method using the following call to linalg.lstsq:

np.linalg.lstsq(left_hand, right_hand, rcond=None)

I need to write a function in C++ that does the same as that Python program. My first approach for this was using OpenCV's cv::solve like this:

cv::solve(equationSystem, equationsRightSide, solutions, cv:ECOMP_SVD);

However, though their results are very similar most of the time, sometimes they differ by quite a lot. For example, for a system with 3 variables:

C++ yields this solution: [0.0008230780686749673; 0.001646156137349933; 8.470329472543003e-20]

Python yields this solution: [-8.76551484e-05 1.41121756e-03 8.62337999e-04]

I understand that this is most likely because cv::solve and np.linalg.lstsq don't use the same algorithm and/or parameters, but cv::solve only allows you to set the matrix decomposition method, while lstsq has many more arguments. Thus, I'm trying to look for alternative functions to implement the C++ version. I have tried Eigen with the following code:

Eigen::JacobiSVD<Eigen::Matrix<double, Eigen:ynamic, Eigen:ynamic, Eigen::RowMajor>> svd(
ES_Eigen, Eigen::ComputeThinU | Eigen::ComputeThinV);
// In Python, if rcond = None, linalg.lstsq calculates the rcond value as the machine's precision
// times the greatest between the equations right side's columns and rows
double rcond = std::numeric_limits<double>::epsilon() * std::max(equationsRightSide.rows, equationsRightSide.cols);
svd.setThreshold(rcond);

Eigen::VectorXd x = svd.solve(RS_Eigen);
std::cout << "Solution with Eigen:\n" << x << std::endl;


Still, this code's solutions are more similar to OpenCV's than lstsq's:

Solution with Eigen:
0.000823078
0.00164616
1.69407e-20


Am I doing something wrong? If neither OpenCv nor Eigen can be used to solve my equation system just like lstsq, is there a C/C++ library that does it? I decided to try Eigen out because I read that lstsq used LAPACK under the hood, so I'm not sure if I should use different parameters for svd.solve().

Also, I made sure to use the same input matrices for C++ and Python, so I don't think the cause of the issue lies there.

Continue reading...