从 boost::ublas 转换为 Eigen 会给出不同的结果

Question

我目前正在将算法从 boost::ublas 移植到 Eigen：

代码 1 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

//}

具有本征的代码 2：

#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

//}

这些功能是更大的代码和功能的一部分，但我认为这两个功能是造成差异的原因。与使用 boost 的代码的输出相比，Eigen 的函数为一些大矩阵提供了不同的输出，我无法理解这些错误。

如有任何帮助，我们将不胜感激。

Answer 1

您没有指定任何输入或您发现的差异是什么。

这让我构建了简单的测试器，我发现“差异”的一个明显来源是 [二进制] 浮点表示的不准确性。

您可以通过一些测试输入轻松确认它： 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";
        }
    } }

}

版画

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

确认所有差异都在（甚至 <）所选数据类型的机器 epsilon 附近。如果你更换那个：

using T = long double;

您得到以下增量：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

从这里去哪里

通过插入您的输入来确定这是否是您的问题。您可能会偶然发现以前未引起您注意的其他事情。如果没有，至少你现在有了工具来提出一个新的、更有针对性的问题。

如果您想了解更多关于浮点不准确的信息：

Why are floating point numbers inaccurate?
Is floating point math broken?

从 boost::ublas 转换为 Eigen 会给出不同的结果

Converting from boost::ublas to Eigen gives different results

c++

boost

matrix

linear-algebra

eigen

从这里去哪里