从 boost::ublas 转换为 Eigen 会得到不同的结果

如何解决从 boost::ublas 转换为 Eigen 会得到不同的结果

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

带有 boost::ublas 的代码 1

#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()),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,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 的函数为某些大型矩阵提供了不同的输出，我无法理解这些错误。

任何帮助将不胜感激。

解决方法

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

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

您可以通过一些测试输入轻松确认：其逆是：

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()),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::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,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,3),(3,1),(2,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.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

确认所有差异都围绕（甚至

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

从这里去哪里

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

如果您想了解有关浮点误差的更多信息：

从 boost::ublas 转换为 Eigen 会得到不同的结果

如何解决从 boost::ublas 转换为 Eigen 会得到不同的结果

解决方法

从这里去哪里

相关推荐