Eigen and Boost matrix inversions give different results

Question

I'm currently building a C++ code that uses matrix inversion. Looking for a library to do this, I noticed that Boost and Eigen3 give different results when inverting my 18x18 matrix ! The difference is not negligible, and therefore leads my algorithm to finally deliver very different results.

Did anyone have the same issue ? What is the most reliable library for matrix inversion ?

Please see below a piece of code that shows how I build the 18x18 matrix, and the difference between the inverted matrices. The infinity-norm of the difference matrix is about 130 with my data. Tested with Eigen-3.3.9 and Boost-1.73.

Note that I already read the following thread : Converting from boost::ublas to Eigen gives different results, but didn't find a real answer to my issue, since I would really like to make a code that is not so dependent on the matrix inversion library I use.

Any help would be greatly appreciated !

#include <Eigen/Dense>

#include <boost/numeric/ublas/lu.hpp>
using namespace boost::numeric;

#include <cmath>
#include <stdexcept>
#include <iostream>
#include <iomanip>

#define SIZE 18

Eigen::MatrixXd initMatrix_Eigen() {
    Eigen::MatrixXd matrix(SIZE, SIZE);
    for (size_t i = 0; i < SIZE / 2; i++) {
        for (size_t j = 0; j < SIZE; j++) {
            matrix(2 * i, j) = pow((double)i + 1, (double)j + 2);
            matrix(2 * i + 1, j) = (j + 2) * pow((double)i + 1, (double)j + 1);
        }
    }
    return matrix;
}

ublas::matrix<double> initMatrix_Boost() {
    ublas::matrix<double> matrix(SIZE, SIZE);
    for (size_t i = 0; i < SIZE / 2; i++) {
        for (size_t j = 0; j < SIZE; j++) {
            matrix(2 * i, j) = pow((double)i + 1, (double)j + 2);
            matrix(2 * i + 1, j) = (j + 2) * pow((double)i + 1, (double)j + 1);
        }
    }
    return matrix;
}

Eigen::MatrixXd invertMatrix_Eigen(const Eigen::MatrixXd& matrix) {
    return matrix.inverse();
}

ublas::matrix<double> invertMatrix_Boost(const ublas::matrix<double>& matrix) {
    if (matrix.size1() != matrix.size2()) throw("Input matrix must be square");
    ublas::matrix<double> tmpMatrix(matrix);
    ublas::permutation_matrix<size_t> permutationMatrix(tmpMatrix.size1());
    if (ublas::lu_factorize(tmpMatrix, permutationMatrix) != 0) throw std::invalid_argument("LU factorization failed");
    ublas::matrix<double> inverse(tmpMatrix.size1(), tmpMatrix.size2());
    inverse.assign(boost::numeric::ublas::identity_matrix<double>(tmpMatrix.size1()));
    ublas::lu_substitute(tmpMatrix, permutationMatrix, inverse);
    return inverse;
}

void main(int argc, char** argv)
{
    std::cout << std::setprecision(12);

    Eigen::MatrixXd matrix_Eigen = initMatrix_Eigen();
    Eigen::MatrixXd inverse_Eigen = invertMatrix_Eigen(matrix_Eigen);
    std::cout << "Eigen inverse matrix :" << std::endl;
    std::cout << inverse_Eigen << std::endl;

    ublas::matrix<double> matrix_Boost = initMatrix_Boost();
    ublas::matrix<double> inverse_Boost_tmp = invertMatrix_Boost(matrix_Boost);
    // Convert inverse_Boost_tmp into an Eigen matrix called inverse_Boost :
    const size_t size = inverse_Boost_tmp.size1();
    Eigen::MatrixXd inverse_Boost(size, size);
    for (size_t i = 0; i < size; i++)
        for (size_t j = 0; j < size; j++)
            inverse_Boost(i, j) = inverse_Boost_tmp(i, j);
    std::cout << "Boost inverse matrix :" << std::endl;
    std::cout << inverse_Boost << std::endl;

    Eigen::MatrixXd difference = inverse_Eigen - inverse_Boost;
    std::cout << "Difference between inverse matrices :" << std::endl;
    std::cout << difference << std::endl;
    std::cout << "NormInf(difference) = " << difference.lpNorm<Eigen::Infinity>() << std::endl;
}