Eigen::Ref 用于连接矩阵

Question

如果我想连接两个矩阵A 和B，我会这样做

using Eigen::MatrixXd;
const MatrixXd A(n, p);
const MatrixXd B(n, q);
MatrixXd X(n, p+q);
X << A, B;

现在如果n、p、q 很大，以这种方式定义X 将意味着创建A 和B 的副本。是否可以将X 定义为Eigen::Ref<MatrixXd>？

谢谢。

Answer 1

不，Ref 不是为此而设计的。我们/您需要为此定义一个新表达式，可以称为Cat。如果您只需要水平连接两个矩阵，在 Eigen 3.3 中，这可以在不到十几行代码中实现为空值表达式，请参阅一些示例 there。

编辑：这是一个独立的示例，表明可以混合矩阵和表达式：

#include <iostream>
#include <Eigen/Core>

using namespace Eigen;

template<typename Arg1, typename Arg2>
struct horizcat_helper {
  typedef Matrix<typename Arg1::Scalar,
    Arg1::RowsAtCompileTime,
    Arg1::ColsAtCompileTime==Dynamic || Arg2::ColsAtCompileTime==Dynamic
    ? Dynamic : Arg1::ColsAtCompileTime+Arg2::ColsAtCompileTime,
    ColMajor,
    Arg1::MaxRowsAtCompileTime,
    Arg1::MaxColsAtCompileTime==Dynamic || Arg2::MaxColsAtCompileTime==Dynamic
    ? Dynamic : Arg1::MaxColsAtCompileTime+Arg2::MaxColsAtCompileTime> MatrixType;
};

template<typename Arg1, typename Arg2>
class horizcat_functor
{
  const typename Arg1::Nested m_mat1;
  const typename Arg2::Nested m_mat2;

public:
  horizcat_functor(const Arg1& arg1, const Arg2& arg2)
    : m_mat1(arg1), m_mat2(arg2)
  {}

  const typename Arg1::Scalar operator() (Index row, Index col) const {
    if (col < m_mat1.cols())
      return m_mat1(row,col);
    return m_mat2(row, col - m_mat1.cols());
  }
};

template <typename Arg1, typename Arg2>
CwiseNullaryOp<horizcat_functor<Arg1,Arg2>, typename horizcat_helper<Arg1,Arg2>::MatrixType>
horizcat(const Eigen::MatrixBase<Arg1>& arg1, const Eigen::MatrixBase<Arg2>& arg2)
{
  typedef typename horizcat_helper<Arg1,Arg2>::MatrixType MatrixType;
  return MatrixType::NullaryExpr(arg1.rows(), arg1.cols()+arg2.cols(),
                                horizcat_functor<Arg1,Arg2>(arg1.derived(),arg2.derived()));
}

int main()
{
  MatrixXd mat(3, 3);
  mat << 0, 1, 2, 3, 4, 5, 6, 7, 8;

  auto example1 = horizcat(mat,2*mat);
  std::cout << example1 << std::endl;

  auto example2 = horizcat(VectorXd::Ones(3),mat);
  std::cout << example2 << std::endl;
  return 0;
}

Answer 2

我将添加 @ggaels horizcat 的 C++14 版本作为答案。实现有点草率，因为它没有考虑 Eigen 编译时常量，但作为回报，它只是一个两行：

auto horizcat = [](auto expr1, auto expr2)
{
    auto get = [expr1=std::move(expr1),expr2=std::move(expr2)](auto row, auto col)
            { return col<expr1.cols() ? expr1(row, col) : expr2(row, col - expr1.cols());};
    return Eigen::Matrix<decltype(get(0,0)), Eigen::Dynamic, Eigen::Dynamic>::NullaryExpr(expr1.rows(), expr1.cols() + expr2.cols(), get);
};

int main()
{
  Eigen::MatrixXd mat(3, 3);
  mat << 0, 1, 2, 3, 4, 5, 6, 7, 8;

  auto example1 = horizcat(mat,2*mat);
  std::cout << example1 << std::endl;

  auto example2 = horizcat(Eigen::MatrixXd::Identity(3,3), mat);
  std::cout << example2 << std::endl;
  return 0;
}

请注意，代码未经测试。

这应该适用于大多数应用程序。但是，如果您使用的是编译时矩阵维度并需要最高性能，请选择 ggaels 答案。在所有其他情况下，也更喜欢 ggaels 回答，因为他是 Eigen 的开发者:-)

Answer 3

我将ggael's answer 扩展为Array 类型、垂直连接和两个以上的参数：

#include <iostream>
#include <Eigen/Core>


namespace EigenCustom
{
using namespace Eigen;

constexpr Index dynamicOrSum( const Index& a, const Index& b ){
    return a == Dynamic || b == Dynamic ? Dynamic : a + b;
}

enum class Direction { horizontal, vertical };

template<Direction direction, typename Arg1, typename Arg2>
struct ConcatHelper {
    static_assert( std::is_same_v<
        typename Arg1::Scalar, typename Arg2::Scalar
    > );
    using Scalar = typename Arg1::Scalar;
    using D = Direction;
    static constexpr Index
        RowsAtCompileTime { direction == D::horizontal ?
            Arg1::RowsAtCompileTime :
            dynamicOrSum( Arg1::RowsAtCompileTime, Arg2::RowsAtCompileTime )
        },
        ColsAtCompileTime { direction == D::horizontal ?
            dynamicOrSum( Arg1::ColsAtCompileTime, Arg2::ColsAtCompileTime ) :
            Arg1::ColsAtCompileTime
        },
        MaxRowsAtCompileTime { direction == D::horizontal ?
            Arg1::MaxRowsAtCompileTime :
            dynamicOrSum( Arg1::MaxRowsAtCompileTime, Arg2::MaxRowsAtCompileTime )
        },
        MaxColsAtCompileTime { direction == D::horizontal ?
            dynamicOrSum( Arg1::MaxColsAtCompileTime, Arg2::MaxColsAtCompileTime ) :
            Arg1::MaxColsAtCompileTime
        };
    
    static_assert(
        (std::is_base_of_v<MatrixBase<Arg1>, Arg1> &&
         std::is_base_of_v<MatrixBase<Arg2>, Arg2> ) ||
        (std::is_base_of_v<ArrayBase<Arg1>, Arg1> &&
         std::is_base_of_v<ArrayBase<Arg2>, Arg2> )
    );
    using DenseType = std::conditional_t<
        std::is_base_of_v<MatrixBase<Arg1>, Arg1>,
        Matrix<
            Scalar,   RowsAtCompileTime,    ColsAtCompileTime,
            ColMajor, MaxRowsAtCompileTime, MaxColsAtCompileTime
        >,
        Array<
            Scalar,   RowsAtCompileTime,    ColsAtCompileTime,
            ColMajor, MaxRowsAtCompileTime, MaxColsAtCompileTime
        >
    >;
};
 
template<Direction direction, typename Arg1, typename Arg2>
class ConcatFunctor
{
    using Scalar = typename ConcatHelper<direction, Arg1, Arg2>::Scalar;
    const typename Arg1::Nested m_mat1;
    const typename Arg2::Nested m_mat2;

public:
    ConcatFunctor(const Arg1& arg1, const Arg2& arg2)
        : m_mat1(arg1), m_mat2(arg2)
    {}

    const Scalar operator() (Index row, Index col) const {
        if constexpr (direction == Direction::horizontal){
            if (col < m_mat1.cols())
                return m_mat1(row,col);
            return m_mat2(row, col - m_mat1.cols());
        } else {
            if (row < m_mat1.rows())
                return m_mat1(row,col);
            return m_mat2(row - m_mat1.rows(), col);
        }
    }
};

template<Direction direction, typename Arg1, typename Arg2>
using ConcatReturnType = CwiseNullaryOp<
    ConcatFunctor<direction,Arg1,Arg2>,
    typename ConcatHelper<direction,Arg1,Arg2>::DenseType
>;

template<Direction direction, typename Arg1, typename Arg2>
ConcatReturnType<direction, Arg1, Arg2>
concat(
    const Eigen::DenseBase<Arg1>& arg1,
    const Eigen::DenseBase<Arg2>& arg2
){
    using DenseType = typename ConcatHelper<direction,Arg1,Arg2>::DenseType;
    using D = Direction;
    return DenseType::NullaryExpr(
        direction == D::horizontal ? arg1.rows() : arg1.rows() + arg2.rows(),
        direction == D::horizontal ? arg1.cols() + arg2.cols() : arg1.cols(),
        ConcatFunctor<direction,Arg1,Arg2>( arg1.derived(), arg2.derived() )
    );
}

template<Direction direction, typename Arg1, typename Arg2, typename ... Ts>
decltype(auto)
concat(
    const Eigen::DenseBase<Arg1>& arg1,
    const Eigen::DenseBase<Arg2>& arg2,
    Ts&& ... rest
){
    return concat<direction>(
        concat<direction>(arg1, arg2),
        std::forward<Ts>(rest) ...
    );
}

template<typename Arg1, typename Arg2, typename ... Ts>
decltype(auto)
concat_horizontal(
    const Eigen::DenseBase<Arg1>& arg1,
    const Eigen::DenseBase<Arg2>& arg2,
    Ts&& ... rest
){
    return concat<Direction::horizontal>(
        arg1, arg2, std::forward<Ts>(rest) ...
    );
}

template<typename Arg1, typename Arg2, typename ... Ts>
decltype(auto)
concat_vertical(
    const Eigen::DenseBase<Arg1>& arg1,
    const Eigen::DenseBase<Arg2>& arg2,
    Ts&& ... rest
){
    return concat<Direction::vertical>(
        arg1, arg2, std::forward<Ts>(rest) ...
    );
}
    

} // namespace EigenCustom

int main()
{
    using namespace Eigen;
    using namespace EigenCustom;
    
    MatrixXd mat(3, 3);
    mat << 0, 1, 2, 3, 4, 5, 6, 7, 8;

    auto example1 = concat_horizontal(mat,2*mat);
    std::cout << "example1:\n" << example1 << '\n';

    auto example2 = concat_horizontal(VectorXd::Ones(3),mat);
    std::cout << "example2:\n" << example2 << '\n';
    
    auto example3 = concat_vertical(mat,RowVectorXd::Zero(3));
    std::cout << "example3:\n" << example3 << '\n';
    
    ArrayXXi arr (2,2);
    arr << 0, 1, 2, 3;
    
    auto example4 = concat_vertical(arr,Array2i{4,5}.transpose());
    std::cout << "example4:\n" << example4 << '\n';
    
    /* concatenating more than two arguments */
    auto example5 = concat_horizontal(mat, mat, mat);
    std::cout << "example5:\n" << example5 << '\n';

    using RowArray2i = Array<int, 1, 2>;
    
    auto example6 = concat_vertical( arr, RowArray2i::Zero(), RowArray2i::Ones() );
    std::cout << "example6:\n" << example6 << '\n';

    return 0;
}