With a recent update from Ubuntu 14.04 to 15.10 one of the applications I maintain failed one of its tests. I have narrowed it down to Eigen and from my tests Eigen 3.2.1 still works as expected, whereas 3.2.2 fails.

The test that fails compares the output of my application to the output generated by R (http://www.r-project.org). What my application basically does is linear (least squares) regression. In the failing test I use a data set in which all points except one being at the same x position, so not a very robust linear system. In this particular test, I used to get values for my output parameters of the order of 6.5, whereas with Eigen > 3.2.2 I get values of the order of 76e6.
For other tests which use less singular data the outcomes are comparable to R's output.

Given that I use LDLT and the fact that the 3.2.2 changelog states:

Various numerical fixes in LDLT, including the case of semi-definite complex matrices.

I'm wondering if this could be the source of my problems.

The code I use is (roughly) the following:

Code: Select all

void linear_regression::LeastSquaredRegression(const <MatrixXd>& X, LDLT<MatrixXd>& Ch) {
    int m = X.cols();
    MatrixXd txx = MatrixXd(m, m).setZero().selfadjointView<Lower>().rankUpdate(X.adjoint());
    Ch = LDLT < MatrixXd > (txx.selfadjointView<Lower>());
    beta = Ch.solve(X.adjoint() * regression_data.Y);                                           // regression_data.Y is a vector
    sigma2 = (regression_data.Y - (X * beta)).squaredNorm();                                    // beta is a vector
}


Here beta is the vector with the problematic values.

Could you try the latest 3.2 branch to see if the issue have already been solved: https://bitbucket.org/eigen/eigen/get/3.2.zip.

thanks.

btw, LDLT is not rank-revealing, so it is not supposed to be stable for singular problems. Better use ColPivHouseholderQR, or SVD in this case.

ggael wrote:Could you try the latest 3.2 branch to see if the issue have already been solved: https://bitbucket.org/eigen/eigen/get/3.2.zip.

thanks.

Thanks for the quick reply!
I tried the following:
- 3.2.0: OK
- 3.2.1: OK
- 3.2.2: Broken
- 3.2.3: Broken
- 3.2.5: Broken
- 3.3-alpha1: Broken
- 3.2 "current" (e059accf541a): Broken

ggael wrote:btw, LDLT is not rank-revealing, so it is not supposed to be stable for singular problems. Better use ColPivHouseholderQR, or SVD in this case.

Strange, I just tried ColPivHouseholderQR and fullPivHouseholderQr with the latest 3.2 branch (e059accf541a) and get the same problem.
Using jacobiSvd() didn't help either. The value are slightly changed, but still incorrect.

hm, could you provide the matrix entries so that we can reproduce.

moreover, all rank-revealing solvers provide a setThreshold method that you can slightly increase to cancel smallest singularvalues, (http://eigen.tuxfamily.org/dox/classEig ... 823edacb2c). E.g.:

ColPivHouseholderQR<MatrixXd> qr;
qr.setThreshold(1e-14);
qr.compute(A),
x = qr.solve(b);

I guessed it worked for your use case with LDLT 3.2.1 because it used by default a rather high threshold which was broken for many other problems... Search for the "L-curve method" to see how to automatically find a best threshold. (not in Eigen).

ggael wrote:hm, could you provide the matrix entries so that we can reproduce.

The input data can be found on pastebin.com:

• y: http://pastebin.com/dM6DQkkb
• X: http://pastebin.com/wnzr8hv6

(I'm trying to solve y = X beta for beta)

ggael wrote:moreover, all rank-revealing solvers provide a setThreshold method that you can slightly increase to cancel smallest singularvalues, (http://eigen.tuxfamily.org/dox/classEig ... 823edacb2c). E.g.:

ColPivHouseholderQR<MatrixXd> qr;
qr.setThreshold(1e-14);
qr.compute(A),
x = qr.solve(b);

I guessed it worked for your use case with LDLT 3.2.1 because it used by default a rather high threshold which was broken for many other problems... Search for the "L-curve method" to see how to automatically find a best threshold. (not in Eigen).

That seems to be the case. Setting the threshold to 1e-14 gave me the same coefficient as R does.

I cannot reproduce using the following test program:

Code: Select all

#include <iostream>
#include <fstream>
#include <string>

#include <Eigen/Dense>
#include <Eigen/Cholesky>

using namespace Eigen;

template<class Matrix>
void init_from_file(Matrix& m, std::string filename) {
    std::ifstream fin(filename);
    for (int i = 0; i < m.rows(); ++i)
        for (int j = 0; j < m.cols(); ++j)
            fin >> m(i,j);
}

int main(int argc, char** argv) {
    MatrixXd D(182,5);
    init_from_file(D, argv[1]);
    VectorXd y(182);
    init_from_file(y, argv[2]);

    MatrixXd A = D.transpose() * D;
    VectorXd b = D.transpose() * y;

    JacobiSVD<MatrixXd> svd;
    svd.compute(D,ComputeThinU|ComputeThinV);

    std::cout << "singular-values:" << svd.singularValues().transpose() << "\n";
    std::cout << "rcond: " << svd.singularValues()(4)/svd.singularValues()(0) << "\n";

    VectorXd x_svd = svd.solve(y);

    VectorXd x_ldlt = (A).ldlt().solve(b);
    VectorXd x_qr = (D).colPivHouseholderQr().solve(y);

    std::cout << "SVD:  " << (D*x_svd-y).norm()/y.norm() << " , " << (A*x_svd-b).norm()/b.norm() << "  ;  " << x_svd.transpose() << "\n";
    std::cout << "LDLT: " << (D*x_ldlt-y).norm()/y.norm() << " , " << (A*x_ldlt-b).norm()/b.norm() << "  ;  " << x_svd.transpose() << "\n";
    std::cout << "QR:   " << (D*x_qr-y).norm()/y.norm() << " , " << (A*x_qr-b).norm()/b.norm() << "  ;  " << x_svd.transpose() << "\n";
   
    return 0;
}


I get:

Code: Select all

singular-values:    781.854      6.8432     4.58092     1.20859 8.28057e-16
rcond: 1.05909e-18
SVD:  0.0369822 , 1.32027e-15  ;     116.106    10.6605 -0.0547384    61.4203    54.6856
LDLT: 0.0369822 , 1.32034e-16  ;     116.106    10.6605 -0.0547384    61.4203    54.6856
QR:   0.0369822 , 5.28106e-16  ;     116.106    10.6605 -0.0547384    61.4203    54.6856


Sorry for not replying earlier. Work is quite busy this week. I hope to have a better look at your code and dig deeper in mine somewhere next week.

In any case: thanks a lot for investing your time into this!