Hi,
I am trying to solve sparse & over-determined system in a least-squares manner (L * px = m). What I did with Eigen is, I apply SimplicialLDLt decomposition to normal matrix (L' * L) and solve L' * m using decomposition.
I could not manage to get right output with large L matrix (error is too high that can not be any numerical error). I also try to emulate Cholesky factorization & solution with Matlab which gives expected output. It is not logical to think that there is something wrong with Eigen's SimplicialLDLt but I don't know what I can do else...
I provide 'working' c++ code for demonstrating the problem (you should just check main method) which can read matrices from txt files, matlab code which is also reading matrices and emulating cholesky solution, also links to problematic L and px txt rars.

• pseudo C++ code for main:

Code: Select all

int main()
{
    Eigen::SparseMatrix<float> L;
    // Read from file...
   
    Eigen::VectorXf px;
    //Read from file...
   
    // Solve least squares problem
    // Equation :  L * px = m
    // L'*L * px = L' *m
    // apply cholesky factorization on (L'*L) and solve (L'*m)
   
    Eigen::SparseMatrix<float> Lt = L.transpose(); // transpose
    Eigen::VectorXf m = L * px;                    // construct problem
    Eigen::SparseMatrix<float> Ls = Lt * L;        // normal matrix

    Eigen::VectorXf mLt = Lt * m;

    Eigen::SimplicialLDLt<Eigen::SparseMatrix<float>> choSparse(Ls); //compute factorization etc.
    if(choSparse.info() != Eigen::Success) { /* FAILED*/ }
   
    Eigen::VectorXf pxd = choSparse.solve(mLt); // Solve
    if(choSparse.info() != Eigen::Success) { /* FAILED*/ }

    cout << (px - pxd) << endl; //differences between expected output and the solution.
}


• All of C++ code:

Code: Select all

#include <Eigen/Sparse>
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <algorithm>

using namespace std;

void readToVector(istream& stream, Eigen::VectorXf& px);

struct PairW{
    double value;
    int    row;
    int    column;
    PairW(double v, int r, int c) : row(r),
        column(c),value(v){}
};

class SparseMatrixWrapper
{
    friend ostream &operator<<(ostream &stream, SparseMatrixWrapper& A);
    friend istream &operator>>(istream &stream, SparseMatrixWrapper& A);
public:   
    SparseMatrixWrapper(int rowCount, int columnCount);
    ~SparseMatrixWrapper(void);
    void          addCell(double value, int row, int column);
    void          fillEigenSparseMatrix(Eigen::SparseMatrix<float>& sparse);
    vector<PairW> cells;
    int           rowCount, columnCount;
};

// THIS IS THE MAIN PART //
/////////////////////////////////////////////////////////////////////////////////////
int main()
{
    // Read L matrix from file
    ifstream inL("L.txt");
    SparseMatrixWrapper wrapper(0,0);
    inL >> wrapper;
    inL.close();
    Eigen::SparseMatrix<float> L(wrapper.rowCount, wrapper.columnCount);
    wrapper.fillEigenSparseMatrix(L);
   
    //Read px vector from file
    Eigen::VectorXf px;
    ifstream inPx("px.txt");
    readToVector(inPx, px);
    inPx.close();
   
    // Solve least squares problem
    // Equation :  L * px = m
    Eigen::SparseMatrix<float> Lt = L.transpose(); //transpose
    Eigen::VectorXf m = L * px; //
    Eigen::SparseMatrix<float> Ls = Lt * L; // normal matrix

    Eigen::VectorXf mLt = Lt * m;

    Eigen::SimplicialLDLt<Eigen::SparseMatrix<float>> choSparse;
    choSparse.compute(Ls);
   
    ofstream pxr("pxDif.txt");
    if(choSparse.info() != Eigen::Success) {
        pxr << "FAILED" << endl;
        pxr.close();
        return 0;
    }
   
    Eigen::VectorXf pxd = choSparse.solve(mLt); // Solve
    if(choSparse.info() != Eigen::Success) {
        pxr << "FAILED" << endl;
        pxr.close();
        return 0;
    }

    pxr << (px - pxd) << endl; //differences between expected output and the solution.
    pxr.close();
}
/////////////////////////////////////////////////////////////////////////////////////

// Utilities for reading/writing matrices to/from files.
void readToVector(istream& stream, Eigen::VectorXf& px)
{
    vector<float> values;
    float value;
    while(stream >> value)
    {
        values.push_back(value);
    }
    px = Eigen::VectorXf(values.size());
    int index = 0;
    for(vector<float>::iterator it = values.begin(); it != values.end(); ++it, ++index)
    {
        px(index) = (*it);
    }
}

 SparseMatrixWrapper::SparseMatrixWrapper(int rowCount, int columnCount)
    : rowCount(rowCount), columnCount(columnCount)
{ }

SparseMatrixWrapper::~SparseMatrixWrapper(void)
{ }

 void SparseMatrixWrapper::addCell(double value, int row, int column)
{
    cells.push_back(PairW(value,row,column));
}

bool comparePairsW (PairW& first, PairW& second)
{
    if(first.column > second.column)
        return false;
    if(first.column < second.column)
        return true;
    if(first.row > second.row)
        return false;
    return true;
}

bool comparePairsW2 (PairW& first, PairW& second)
{
    if(first.row > second.row)
        return false;
    if(first.row < second.row)
        return true;
    if(first.column > second.column)
        return false;
    return true;
}

void SparseMatrixWrapper::fillEigenSparseMatrix(Eigen::SparseMatrix<float>& sparse)
{
    sort(cells.begin(), cells.end(), comparePairsW);
    sparse.reserve(cells.size()); 
    int columnIndex = 0;
    sparse.startVec(0);
    for(vector<PairW>::iterator it = cells.begin(); it != cells.end(); ++it)
    {
        if(columnIndex < it->column)
        {
            columnIndex++;
            sparse.startVec(columnIndex);
        }
        sparse.insertBack(it->row,columnIndex) = it->value;
    }
    sparse.finalize(); 
}

ostream &operator<<(ostream& stream, SparseMatrixWrapper& A)
{
    sort(A.cells.begin(), A.cells.end(), comparePairsW2);
    vector<PairW>::iterator it = A.cells.begin();
    for(int r = 0; r < A.rowCount; ++r)
    {
        for(int c = 0; c < A.columnCount; ++c)
        {
            if(it != A.cells.end() && (*it).row == r && (*it).column == c)
            {
                stream << (*it).value;
                ++it;
            }
            else
            {
                stream << "0";
            }
            if(c != A.columnCount -1)
            {
                stream << " ";
            }
        }
        stream << endl;
    }
    return stream;
}

istream &operator>>(istream& stream, SparseMatrixWrapper& A)
{
    string line;
    float value;
    int column = 0;
    int row = 0;
    while(getline(stream, line))
    {
        istringstream lineStream(line);
        while(lineStream >> value)
        {
            if(value != 0)
            {
                A.addCell(value, row, column);
            } 
            column++;
        }
        A.columnCount = column;
        column = 0;
        row++;
    }
    A.rowCount = row;
    return stream;
}


• Matlab code:

Variable names and expression orders are similar to main method in C++ code (matlab counterpart of main method).

Code: Select all

L = dlmread('L.txt');
px = dlmread('px.txt');
L = sparse(L);
Lt = L';
m = L*px;
mLt = Lt * m;
Ls = Lt*L;
R = chol(Ls);
pxd = R'\mLt;
pxd = R\pxd;
a = pxd - px;


• Rar files for L matrix and px vector:

http://dl.dropbox.com/u/54262186/L.rar
http://dl.dropbox.com/u/54262186/px.rar

Note:
I append reading & writing utility (actually all of the code) recently (previous hour) so other people can try the code. I test it using Matlab, before solution step all matrices and vectors in Matlab and C++ are matching. So there should not be any problem related to reading & writing utilities.

I am very sorry for long question, any help will be appreciated. Thanks for reading.

there is nothing wrong, it's just the limitations of float (what you used) versus double (what Matlab uses).

For instance, with float, the infinite norm (Ls * pxd - mLt).cwiseAbs().maxCoeff() gives 3.68746e-07 which is the best accuracy you can get with floats, leading to:
(L * pxd - m).cwiseAbs().maxCoeff() = 0.00026805
(px - pxd).cwiseAbs().maxCoeff()/px.cwiseAbs().maxCoeff() = 0.300609

With double this drops down to:
(Ls * pxd - mLt).cwiseAbs().maxCoeff() = 8.68229e-16
(L * pxd - m).cwiseAbs().maxCoeff() = 5.95531e-14
(px - pxd).cwiseAbs().maxCoeff()/px.cwiseAbs().maxCoeff() = 6.48869e-11