c++matrixboostlinear-algebraeigen
Converting from boost::ublas to Eigen gives different results
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.
Solution
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 
:
#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:
- C typedef name conflict
- How to use 32-bit pointers in 64-bit application?
- problems with scanf input stream
- Is it guaranteed to be safe to perform memcpy(0,0,0)?
- Need help understanding pointer to a type of n elements in C (int (*p)[10])?
- Local variables performance vs global in example
- How to Print SOAP Message Request Body generated using GSOAP when calling service
- Big Number Subtraction in C
- Docker custom program stop
- GTK4 Windows Multiple File Choosing
- Cannot get lldb to read file input through redirect
- Modeling nested lists in C
- Linker and Memory map for shared memory in autosar embedded environment
- Why is only the numerator cast to float to obtain a float quotient while dividing two integers?
- $gcc problemMatcher doesn't show all errors or warnings
- What registers to save in the ARM C calling convention?
- What does ** do in C language?
- How to install pywin32 module in windows 7
- How to differentiate a Cypher clause from an SQL clause in C?
- How do I check if my executable is running for a service or for a normal program?
- How to fix a simple error in Atmel studio
- Given day, month, year calculate day in year
- Definition of int64_t
- Can we have multiple c files to define functions for one header file?
- Check if array is ASCII
- What does C or GCC do when implicit function declaration exists
- How can I programmatically find the CPU frequency with C
- Can I use a binary literal in C or C++?
- Want to Convert Integer to String without itoa function
- How to permanently configure VSCode to use cl.exe without launching it from the Developer Command Prompt?