[Eigen中文文档] 扩展/自定义Eigen(三)

简介: 本页面针对非常高级的用户,他们不害怕处理一些Eigen的内部细节。在大多数情况下,可以通过使用自定义一元或二元函数避免使用自定义表达式,而极其复杂的矩阵操作可以通过零元函数(nullary-expressions)来实现,如前一页所述。本页面通过示例介绍了如何在Eigen中实现新的轻量级表达式类型。它由三个部分组成:表达式类型本身、包含有关表达式编译时信息的特性类和评估器类,用于将表达式评估为矩阵。

文档总目录

添加新的表达式类型

英文原文(Adding a new expression type)

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

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

设定

循环矩阵是每列与左边的列相同的矩阵,只不过它循环向下移位。例如,这是一个 4×4 循环矩阵:
$$ \begin{bmatrix} 1&8&4&2 \\ 2&1&8&4 \\ 4&2&1&8 \\ 8&4&2&1 \end{bmatrix} $$
循环矩阵由其第一列唯一确定。我们希望编写一个函数 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
相关文章
|
存储 编译器
[Eigen中文文档] 深入了解 Eigen - 类层次结构
本页面介绍了Eigen类层次结构中 Core 类的设计及其相互关系。一般用户可能不需要关注这些细节,但对于高级用户和Eigen开发人员可能会有用。
291 0
|
存储 C语言 C++
|
存储 编译器
|
存储 算法 NoSQL
[Eigen中文文档] 稀疏矩阵操作
在许多应用中(例如,有限元方法),通常要处理非常大的矩阵,其中只有少数系数不为零。在这种情况下,可以通过使用仅存储非零系数的特殊表示来减少内存消耗并提高性能。这样的矩阵称为稀疏矩阵。
463 0
[Eigen中文文档] 编写以特征类型为参数的函数(一)
Eigen使用表达式模板的方式导致每个表达式的类型可能都不同。如果将这样的表达式传递给一个需要Matrix类型参数的函数,则表达式将隐式地被评估为一个临时Matrix,然后再传递给函数。这意味着失去了表达式模板的好处。
144 0
|
存储 NoSQL API
[Eigen中文文档] Matrix类
在Eigen中,所有矩阵和向量都是Matrix模板类的对象。向量只是行数或者列数为1的特殊矩阵。
399 1
[Eigen中文文档] 在 CMake 项目中使用 Eigen
Eigen提供了CMake(CMake 3.0或更高版本)支持,使得该库可以轻松地在CMake项目中使用。
702 1
|
存储 对象存储 索引
[Eigen中文文档] 扩展/自定义Eigen(二)
CwiseNullaryOp 类的主要目的是定义过程矩阵,例如由Ones()、Zero()、Constant()、Identity()和Random()方法返回的常量或随机矩阵。然而,通过一些想象力,可以用最小的努力实现非常复杂的矩阵操作,因此很少需要实现新的表达式。
118 0
|
测试技术 API C++
[Eigen中文文档] 扩展/自定义Eigen(一)
在本节中,将介绍如何向MatrixBase添加自定义方法。由于所有表达式和矩阵类型都继承自MatrixBase,因此将方法添加到MatrixBase会立即使其对所有表达式可用!一个典型的用例是,使Eigen与另一个API兼容。
271 0
|
存储 编译器 调度