在 Eigen 中重载 + 运算符

Question

我需要像这样重载 + 操作 a + b = max(a,b)。 因此，我需要它来处理矩阵加法和矩阵乘法以及其他一些操作（跟踪、功率等）。这里 \bigoplus 是 ma​​x 操作 使用 Eigen 执行此操作的最佳方法是什么？我读到了特征扩展here，但我不明白如何为我的任务做到这一点。 目前我有这个：

    #include <iostream>
    
    #include <Eigen/Dense>
    #include <algorithm>
    using namespace Eigen;
    
    namespace MaxAlgebra {
    
        template <typename T>
        T operator+(const T& a,const T& b) {
            T c(a.rows(),a.cols());
            for (uint i = 0; i < a.rows(); ++i) {
                for (uint j = 0; j < a.cols(); ++j) {
                    c(i,j) = std::max(a(i,j),b(i,j));
                }
            }
    
            return c;
        }
        template <typename T>
        T operator*(const T& a,const T& b){
            T c(a.rows(),b.cols());
            for (uint i = 0; i < a.rows(); ++i) {
                for (uint j = 0; j < b.cols(); ++j) {
                    std::vector<uint> values;
                    for (uint k = 0; k < a.cols(); ++k) {
                        values.push_back(a(i,k) * b(k,j));
                    }
    
                    c(i,j) = *std::max_element(begin(values),end(values));
                }
            }
    
            return c;
    
        }
        template <typename T>
        uint trace(const T& a) {
            std::vector<uint> values;
            for (uint i = 0; i < a.rows(); ++i) {
                values.push_back(a(i,i));
            }
            return *std::max_element(begin(values),end(values));
        }
    }
    
    
    int main() {
    
        MatrixXd x(2,2);
        MatrixXd y(2,2);
        
        x(0,0) = 3;
        x(1,0) = 2;
        x(0,1) = 1;
        x(1,1) = 2;
    
    
        y(0,0) = 2;
        y(1,0) = 1;
        y(0,1) = 2;
        y(1,1) = 3;
        MatrixXd c = MaxAlgebra::operator*(x,y);
        std::cout << "Here is the matrix a:\n" << x << std::endl;
        std::cout << "Here is the matrix b:\n" << y << std::endl;
        std::cout << "Here is the matrix c:\n" << c << std::endl;
    
    
    
        return 0;
    }

Answer 1

如果我正确理解了您的代数，您可以简单地创建一个自定义标量类型。这似乎有效：

template<class T>
struct MaxAlg
{
  T scalar;

  MaxAlg() = default;
  MaxAlg(T scalar) noexcept // implicit conversion for convenience
    : scalar(scalar)
  {}
  explicit operator T() const noexcept
  { return scalar; }

  MaxAlg& operator+=(MaxAlg o) noexcept
  {
    scalar = std::max(scalar, o.scalar);
    return *this;
  }
  friend MaxAlg operator+(MaxAlg left, MaxAlg right) noexcept
  { left += right; return left; }

  MaxAlg& operator*=(MaxAlg o) noexcept
  {
    scalar *= o.scalar;
    return *this;
  }
  friend MaxAlg operator*(MaxAlg left, MaxAlg right) noexcept
  { left *= right; return left; }

  friend bool operator==(MaxAlg left, MaxAlg right) noexcept
  { return left.scalar == right.scalar; }

  friend bool operator!=(MaxAlg left, MaxAlg right) noexcept
  { return left.scalar != right.scalar; }

  friend bool operator<(MaxAlg left, MaxAlg right) noexcept
  { return left.scalar < right.scalar; }

  friend bool operator<=(MaxAlg left, MaxAlg right) noexcept
  { return left.scalar <= right.scalar; }

  friend bool operator>(MaxAlg left, MaxAlg right) noexcept
  { return left.scalar > right.scalar; }

  friend bool operator>=(MaxAlg left, MaxAlg right) noexcept
  { return left.scalar >= right.scalar; }

  friend std::ostream& operator<<(std::ostream& left, MaxAlg right)
  { return left << right.scalar; }
};

template<class T, Eigen::Index Rows, Eigen::Index Cols>
using MaxAlgMatrix = Eigen::Matrix<MaxAlg<T>, Rows, Cols>;

template<class T, Eigen::Index Rows, Eigen::Index Cols>
using MaxAlgArray = Eigen::Array<MaxAlg<T>, Rows, Cols>;

using MaxAlgMatrixXd = MaxAlgMatrix<double, Eigen::Dynamic, Eigen::Dynamic>;
using MaxAlgVectorXd = MaxAlgMatrix<double, Eigen::Dynamic, 1>;
using MaxAlgArrayXXd = MaxAlgArray<double, Eigen::Dynamic, Eigen::Dynamic>;
using MaxAlgArrayXd = MaxAlgArray<double, Eigen::Dynamic, 1>;


int main()
{
  Eigen::MatrixXd a = Eigen::MatrixXd::Random(10, 10);
  Eigen::MatrixXd b = Eigen::MatrixXd::Random(10, 10);
  MaxAlgMatrixXd maxalg_a = a.cast<MaxAlg<double> >();
  MaxAlgMatrixXd maxalg_b = b.cast<MaxAlg<double> >();
  std::cout << (maxalg_a * maxalg_b).cast<double>() << "\n\n";
  std::cout << (maxalg_a.array() + maxalg_b.array()).cast<double>() << "\n\n";
  std::cout << a.cwiseMax(b) << "\n\n";
}

这会禁用 Eigen 的矢量化，但是当您使用 -O3 编译时编译器仍然可以执行此操作，而且工作量要少得多。