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:
whose inverse is
:
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: