[Eigen中文文档] 扩展/自定义Eigen（三）

文档总目录

添加新的表达式类型

英文原文(Adding a new expression type)

本页面针对非常高级的用户，他们不害怕处理一些Eigen的内部细节。在大多数情况下，可以通过使用自定义一元或二元函数避免使用自定义表达式，而极其复杂的矩阵操作可以通过零元函数(nullary-expressions)来实现，如前一页所述。

本页面通过示例介绍了如何在Eigen中实现新的轻量级表达式类型。它由三个部分组成：表达式类型本身、包含有关表达式编译时信息的特性类和评估器类，用于将表达式评估为矩阵。

设定

循环矩阵是每列与左边的列相同的矩阵，只不过它循环向下移位。例如，这是一个 4×4 循环矩阵：
$$\left[\begin{array}{cccc}1& 8& 4& 2\\ 2& 1& 8& 4\\ 4& 2& 1& 8\\ 8& 4& 2& 1\end{array}\right]$$
循环矩阵由其第一列唯一确定。我们希望编写一个函数 makeCirculant，在给定第一列的情况下，它返回一个表示循环矩阵的表达式。

为简单起见，我们将 makeCirculant 函数限制为稠密矩阵。也允许数组或稀疏矩阵可能是有意义的，但我们不会在这里这样做。我们也不希望支持矢量化。

开始

我们将分部分介绍实现makeCirculant函数的文件。首先，我们包含适当的头文件并向前声明我们将称之为Circulant的表达式类。makeCirculant函数将返回此类型的对象。实际上，类Circulant是一个类模板；模板参数ArgType是传递给makeCirculant函数的向量的类型。

#include <Eigen/Core>
#include <iostream>
 
template <class ArgType> class Circulant;

Traits类

对于每个表达式类 X，在 Eigen::internal 命名空间中应该有一个特征类 Traits<X>，其中包含有关 X 的信息（称为编译时）。

如“设定”中所解释的，我们设计了Circulant表达式类来引用密集矩阵。循环矩阵的条目与传递给makeCirculant函数的向量的条目具有相同的类型。用于索引条目的类型也相同。为简单起见，我们仅返回列主矩阵。最后，循环矩阵是一个方阵（行数等于列数），行数等于传递给makeCirculant函数的列向量的行数。如果这是一个动态大小向量，则循环矩阵的大小在编译时是未知的。

namespace Eigen {
  namespace internal {
    template <class ArgType>
    struct traits<Circulant<ArgType> >
    {
      typedef Eigen::Dense StorageKind;
      typedef Eigen::MatrixXpr XprKind;
      typedef typename ArgType::StorageIndex StorageIndex;
      typedef typename ArgType::Scalar Scalar;
      enum { 
        Flags = Eigen::ColMajor,
        RowsAtCompileTime = ArgType::RowsAtCompileTime,
        ColsAtCompileTime = ArgType::RowsAtCompileTime,
        MaxRowsAtCompileTime = ArgType::MaxRowsAtCompileTime,
        MaxColsAtCompileTime = ArgType::MaxRowsAtCompileTime
      };
    };
  }
}

表达式类

接下来的步骤是定义表达式类本身。要继承自MatrixBase，以公开密集矩阵的接口。在构造函数中，我们检查是否传递了列向量（参见断言），并将用于构建循环矩阵的向量存储在成员变量m_arg中。最后，表达式类应计算相应循环矩阵的大小。如上所述，这是一个方阵，其列数与用于构造矩阵的向量一样多。

template <class ArgType>
class Circulant : public Eigen::MatrixBase<Circulant<ArgType> >
{
public:
  Circulant(const ArgType& arg)
    : m_arg(arg)
  { 
    EIGEN_STATIC_ASSERT(ArgType::ColsAtCompileTime == 1,
                        YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX);
  }
 
  typedef typename Eigen::internal::ref_selector<Circulant>::type Nested; 
 
  typedef Eigen::Index Index;
  Index rows() const { return m_arg.rows(); }
  Index cols() const { return m_arg.rows(); }
 
  typedef typename Eigen::internal::ref_selector<ArgType>::type ArgTypeNested;
  ArgTypeNested m_arg;
};

评估器

最后一个大片段实现了Circulant表达式的评估器。评估器计算循环矩阵的条目。这是在.coeff()成员函数中完成的。通过找到构成循环矩阵的向量的相应条目来计算这些条目。当循环矩阵是从给定复杂表达式的向量构建时，获取此条目可能实际上并不容易，因此我们使用对应于该向量的评估器。

CoeffReadCost 常量记录计算循环矩阵一项的成本；我们忽略索引计算，并认为这与计算构建循环矩阵的向量条目的成本相同。

在构造函数中，我们保存定义循环矩阵的列向量的求值器。还保存该向量的大小；请记住，我们可以查询表达式对象来查找大小，但不能查询求值器。

namespace Eigen {
  namespace internal {
    template<typename ArgType>
    struct evaluator<Circulant<ArgType> >
      : evaluator_base<Circulant<ArgType> >
    {
      typedef Circulant<ArgType> XprType;
      typedef typename nested_eval<ArgType, XprType::ColsAtCompileTime>::type ArgTypeNested;
      typedef remove_all_t<ArgTypeNested> ArgTypeNestedCleaned;
      typedef typename XprType::CoeffReturnType CoeffReturnType;
 
      enum { 
        CoeffReadCost = evaluator<ArgTypeNestedCleaned>::CoeffReadCost,
        Flags = Eigen::ColMajor 
      };
      
      evaluator(const XprType& xpr)
        : m_argImpl(xpr.m_arg), m_rows(xpr.rows())
      { }
 
      CoeffReturnType coeff(Index row, Index col) const
      {
        Index index = row - col;
        if (index < 0) index += m_rows;
        return m_argImpl.coeff(index);
      }
 
      evaluator<ArgTypeNestedCleaned> m_argImpl;
      const Index m_rows;
    };
  }
}

入口

template <class ArgType>
Circulant<ArgType> makeCirculant(const Eigen::MatrixBase<ArgType>& arg)
{
  return Circulant<ArgType>(arg.derived());
}

完成这些之后，makeCirculant 函数就非常简单了。它只是创建一个表达式对象并返回它。

一个简单用于测试的 main 函数

最后，一个简短的主函数展示了如何调用 makeCirculant 函数。

int main()
{
    Eigen::VectorXd vec(4);
    vec << 1, 2, 4, 8;
    Eigen::MatrixXd mat;
    mat = makeCirculant(vec);
    std::cout << mat << std::endl;
}

如果组合所有代码片段，则会产生以下输出，表明程序按预期工作：

1 8 4 2
2 1 8 4
4 2 1 8
8 4 2 1