1

I am currently porting an algorithm from boost::ublas to Eigen:

Code 1 with boost::ublas

#ifndef KHACH_H
#define KHACH_H

#include <set>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>

#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/lu.hpp>

#include <iostream>
#include <boost/numeric/ublas/io.hpp>

//namespace Minim {

  namespace ublas=boost::numeric::ublas;

  template<class T>
  bool InvertMatrix(const ublas::matrix<T> &input,
                    ublas::matrix<T> &inverse)
  {
    using namespace boost::numeric::ublas;
    typedef permutation_matrix<std::size_t> pmatrix;
    matrix<T> A(input);
    pmatrix pm(A.size1());
    int res = lu_factorize(A,pm);
    if( res != 0 ) return false;
    inverse.assign(ublas::identity_matrix<T>(A.size1()));
    lu_substitute(A, pm, inverse);
    return true;
  }


  inline void Lift(const ublas::matrix<double> &A,
            ublas::matrix<double> &Ap)
  {
    Ap.resize(A.size1()+1,
              A.size2());
    ublas::matrix_range<ublas::matrix<double> >
      sub(Ap,
          ublas::range(0, A.size1()),
          ublas::range(0, A.size2()));
    sub.assign(A);
    ublas::row(Ap, Ap.size1()-1)=ublas::scalar_vector<double>(A.size2(),1.0);

  }

#endif

//}

Code 2 with Eigen:

#ifndef KHACH_H
#define KHACH_H

#include <set>
#include <iostream>
#include <Eigen/Eigen>

//namespace Minim {

  template <class NT>
  using MT = Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic>;

  template <class NT>
  using VT = Eigen::Matrix<NT, Eigen::Dynamic, 1>;

  template<typename Derived>
  inline bool is_nan(const Eigen::MatrixBase<Derived>& x)
  {
  return ((x.array() == x.array())).all();
  }

  template<class T>
  bool InvertMatrix(const MT<T> &input,
                    MT<T> &inverse)
  {
    inverse.setIdentity(input.rows(), input.cols());
    inverse = input.inverse();
    return !is_nan(inverse);
  }

  inline void Lift(const MT<double> &A, MT<double> &Ap)
  {
    Ap.resize(A.rows()+1, A.cols());
    Ap.topLeftCorner(A.rows(), A.cols()) = A;
    Ap.row(Ap.rows()-1).setConstant(1.0); 
  }
  

#endif

//}

These functions are part of the bigger code and functionality, but I think these two functions are the ones creating the difference. The functions with Eigen are giving a different output for some large matrices compared to the output of the code using boost, I am not able to understand the bugs.

Any help would be appreciated.

  • What is an example input that gives different output? In fact, please add a reproducer with those inputs, so you don't have to wait for someone experienced in both libraries to make up the code – sehe Feb 01 '21 at 13:50

1 Answers1

1

You didn't specify any inputs or what the discrepancy is you're finding.

This lead me to build simple testers, in which I find that an obvious source of "differences" is the inaccuracy of [binary] floating point representations.

You can easily confirm it with some test input: enter image description here whose inverse is enter image description here:

Live On Compuler Explorer

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <set>

#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/triangular.hpp>

#include <boost/numeric/ublas/io.hpp>
#include <iostream>

namespace Minim1 {

namespace ublas = boost::numeric::ublas;

template <class T> using MT = ublas::matrix<T>;

template <class T> bool InvertMatrix(const MT<T>& input, MT<T>& inverse)
{
    using namespace boost::numeric::ublas;
    typedef permutation_matrix<std::size_t> pmatrix;
    matrix<T> A(input);
    pmatrix pm(A.size1());
    int res = lu_factorize(A, pm);
    if (res != 0)
        return false;
    inverse.assign(ublas::identity_matrix<T>(A.size1()));
    lu_substitute(A, pm, inverse);
    return true;
}

template <class T>
inline void Lift(const ublas::matrix<T>& A, ublas::matrix<T>& Ap)
{
    Ap.resize(A.size1() + 1, A.size2());
    ublas::matrix_range<ublas::matrix<T>> sub(
        Ap, ublas::range(0, A.size1()), ublas::range(0, A.size2()));
    sub.assign(A);
    ublas::row(Ap, Ap.size1() - 1) = ublas::scalar_vector<T>(A.size2(), 1.0);
}

}

#include <Eigen/Eigen>
#include <iostream>
#include <set>

namespace Minim2 {

template <class T>
using MT = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;

static_assert(Eigen::RowMajor == 1);
template <class T>
using VT = Eigen::Matrix<T, Eigen::Dynamic, Eigen::RowMajor>;

template <typename Derived>
inline bool is_nan(const Eigen::MatrixBase<Derived>& x)
{
    return ((x.array() == x.array())).all();
}

template <class T> bool InvertMatrix(const MT<T>& input, MT<T>& inverse)
{
    inverse.setIdentity(input.rows(), input.cols());
    inverse = input.inverse();
    return !is_nan(inverse);
}

template <typename T>
inline void Lift(const MT<T>& A, MT<T>& Ap)
{
    Ap.resize(A.rows() + 1, A.cols());
    Ap.topLeftCorner(A.rows(), A.cols()) = A;
    Ap.row(Ap.rows() - 1).setConstant(1.0);
}

}

template <typename T>
static inline std::string compare(Minim1::MT<T> const& a, Minim2::MT<T> const& b) {
    if (a.size1() != static_cast<size_t>(b.rows())) return "rows do not match";
    if (a.size2() != static_cast<size_t>(b.cols())) return "cols do not match";
    for (size_t r = 0; r < a.size1(); r++) {
        for (size_t c = 0; c < a.size2(); c++) {
            auto va = a(r,c);
            auto vb = b(r,c);
            auto delta = std::abs(va-vb);
            if (va != vb) {
                std::ostringstream oss;
                oss 
                    << "mismatch at (" << r << ", " << c << "): " 
                    << va << " != " << vb 
                    << " delta:" << std::abs(va-vb)
                    << " significant? " << std::boolalpha
                        << (std::numeric_limits<T>::epsilon() < delta) << "\n";
                return oss.str();
            }
        }
    }
    return "equivalent";
}

template <typename T>
auto convert(Minim1::MT<T> const& a) {
    Minim2::MT<T> b(a.size1(), a.size2());
    for (size_t r = 0; r < a.size1(); r++) {
    for (size_t c = 0; c < a.size2(); c++) {
        b(r, c) = a(r, c);
    } }

    return b;
}

int main() {
    using T = double;
    using M1 = Minim1::MT<T>;
    using M2 = Minim2::MT<T>;

    auto report = [](auto const& a, auto const& b) {
        std::cout << "\na: ------\n" << a;
        std::cout << "\nb: ------\n" << b;
        std::cout << "\n" << compare(a, b) << "\n";
    };

    M1 a(3, 3);
    a(0, 0) = 1; a(0, 1) = 2; a(0, 2) = 3;
    a(1, 0) = 3; a(1, 1) = 2; a(1, 2) = 1;
    a(2, 0) = 2; a(2, 1) = 1; a(2, 2) = 3;

    M2 b(3, 3);
    b << 1, 2, 3,
        3, 2, 1,
        2, 1, 3;
    report(a, b);

    std::cout << "\nINVERSIONS";
    M1 ai(a.size1(), a.size2());
    M2 bi(b.rows(), b.cols());
    Minim1::InvertMatrix(a, ai);
    Minim2::InvertMatrix(b, bi);

    report(ai, bi);

    M2 deltas = (convert(ai) - bi).cwiseAbs();
    constexpr auto eps = std::numeric_limits<T>::epsilon();
    std::cout << "deltas:\n" << deltas << "\n";
    for (int r = 0; r < deltas.rows(); r++) {
    for (int c = 0; c < deltas.cols(); c++) {
        auto d = deltas(r,c);
        if (d > eps) {
            std::cout << "The delta at (" << r << ", " << c << ") (" << d << " is > ε (" << eps << ")\n";
        }
    } }

}

Prints

a: ------
[3,3]((1,2,3),(3,2,1),(2,1,3))
b: ------
1 2 3
3 2 1
2 1 3
equivalent

INVERSIONS
a: ------
[3,3]((-0.416667,0.25,0.333333),(0.583333,0.25,-0.666667),(0.0833333,-0.25,0.333333))
b: ------
-0.416667      0.25  0.333333
 0.583333      0.25 -0.666667
0.0833333     -0.25  0.333333
mismatch at (0, 0): -0.416667 != -0.416667 delta:5.55112e-17 significant? false

deltas:
5.55112e-17           0           0
          0 2.77556e-17           0
          0 2.77556e-17           0

Confirming that all differences are around (even <) the machine epsilon for the chosen data type. If you replace that one:

using T = long double;

You get the following deltas: Compiler Explorer

mismatch at (0, 0): -0.416667 != -0.416667 delta:2.71051e-20 significant? false

deltas:
2.71051e-20 1.35525e-20           0
5.42101e-20           0           0
6.77626e-21           0           0

Where To Go From Here

Find out whether this is your problem by plugging in your inputs. You might stumble on other things that escaped your attention before. If not, at least you now have the tools to make a new, more focused question.

If you want to learn more about floating point inaccuracy:

sehe
  • 374,641
  • 47
  • 450
  • 633
  • Come to think of it I should probably have demonstrated how to express `compare()` in terms of tolerances: this means that the inverses will actually be reported as _equivalent_: https://godbolt.org/z/WjK5ar. In practice you may want to increase tolerances. – sehe Feb 01 '21 at 15:19
  • Thanks for the help @sehe. I will test the code on my system and let you know my results. But I have a problem in one thing, the functions gave same results in many cases, but only 2 cases gave significant different outputs, like one was 6e+20 and other gave 8e+20. How can I test and validate such error ? It is possible that any of the inputs that I test on gives same results. – Vaibhav Thakkar Feb 01 '21 at 19:56
  • "but only 2 cases gave significant different outputs" - just find out what the inputs are there. Plug them in if you can. I accept ig you have to put it in a file on a pastebin/gist or something. – sehe Feb 01 '21 at 20:41