[Eigen中文文档] 编写以特征类型为参数的函数(二)

简介: Eigen使用表达式模板的方式导致每个表达式的类型可能都不同。如果将这样的表达式传递给一个需要Matrix类型参数的函数,则表达式将隐式地被评估为一个临时Matrix,然后再传递给函数。这意味着失去了表达式模板的好处。

文档总目录

[Eigen中文文档] 编写以特征类型为参数的函数(一)

英文原文(Writing Functions Taking Eigen Types as Parameters)

采用普通矩阵或数组参数的函数在哪些情况下有效?

如果不使用模板函数,并且没有 Ref 类,先前 cov 函数的简单实现可能如下所示:

MatrixXf cov(const MatrixXf& x, const MatrixXf& y)
{
   
  const float num_observations = static_cast<float>(x.rows());
  const RowVectorXf x_mean = x.colwise().sum() / num_observations;
  const RowVectorXf y_mean = y.colwise().sum() / num_observations;
  return (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

与一开始想的相反,除非需要一个可以与双矩阵一起使用的通用实现,否则这个实现是可以的,除非不关心临时对象。为什么会这样?哪些地方涉及到临时变量?下面给出的代码是如何编译的?

MatrixXf x,y,z;
MatrixXf C = cov(x,y+z);

在这种特殊情况下,示例是可以的并且能够工作,因为两个参数都声明为const引用。编译器会创建一个临时变量,将表达式x+z计算到这个临时变量中。一旦函数运行完成,临时变量就会被释放,并将结果分配给C

注意:对 Matrix(或 Array)进行 const 引用的函数可以以临时变量为代价来处理表达式。

在哪些情况下,采用普通矩阵或数组参数的函数会失败?

在这里,我们考虑上面给出的函数的稍微修改版。这一次,我们不想返回结果,而是传递一个额外的非常量参数,以便我们可以存储结果。一个第一次的朴素实现可能如下所示:

// Note: This code is flawed!
void cov(const MatrixXf& x, const MatrixXf& y, MatrixXf& C)
{
   
  const float num_observations = static_cast<float>(x.rows());
  const RowVectorXf x_mean = x.colwise().sum() / num_observations;
  const RowVectorXf y_mean = y.colwise().sum() / num_observations;
  C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

当尝试执行以下代码时:

MatrixXf C = MatrixXf::Zero(3,6);
cov(x,y, C.block(0,0,3,3));

编译器会失败,因为将 MatrixXf::block() 返回的表达式转换为非常量 MatrixXf& 是不可能的。这是因为编译器希望避免将结果写入临时对象的风险。在这种特殊情况下,这种保护并不适用 - 如果想要写入临时对象。那么我们如何解决这个问题呢?

目前首选的解决方案是基于一点技巧。需要将常量引用传递给矩阵,并且需要在内部丢弃常量。 C98 兼容编译器的正确实现是:

template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> const & C)
{
   
  typedef typename Derived::Scalar Scalar;
  typedef typename internal::plain_row_type<Derived>::type RowVectorType;

  const Scalar num_observations = static_cast<Scalar>(x.rows());

  const RowVectorType x_mean = x.colwise().sum() / num_observations;
  const RowVectorType y_mean = y.colwise().sum() / num_observations;

  const_cast< MatrixBase<OtherDerived>& >(C) =
    (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

上面的实现现在不仅可以使用临时表达式,而且还允许使用具有任意浮点标量类型的矩阵的函数。

注意:const cast 仅适用于模板化函数。它不适用于 MatrixXf实现,因为无法将 Block 表达式转换为 Matrix引用!

如何在通用实现中调整矩阵大小?

为了使协方差函数具有普遍适用性,有如下代码:

MatrixXf x = MatrixXf::Random(100,3);
MatrixXf y = MatrixXf::Random(100,3);
MatrixXf C;
cov(x, y, C);

当我们使用以 MatrixBase 为参数的实现时,情况就不同了。一般来说,Eigen 支持自动调整大小,但不能在表达式上这样做。为什么要允许重调矩阵块大小?它是一个子矩阵的引用,我们绝对不想重调它的大小。那么,如果我们不能在 MatrixBase 上调整大小,如何合并调整大小呢?解决方法是调整派生对象的大小,如此实现:

template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> const & C_)
{
   
  typedef typename Derived::Scalar Scalar;
  typedef typename internal::plain_row_type<Derived>::type RowVectorType;

  const Scalar num_observations = static_cast<Scalar>(x.rows());

  const RowVectorType x_mean = x.colwise().sum() / num_observations;
  const RowVectorType y_mean = y.colwise().sum() / num_observations;

  MatrixBase<OtherDerived>& C = const_cast< MatrixBase<OtherDerived>& >(C_);

  C.derived().resize(x.cols(),x.cols()); // resize the derived object
  C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

这个实现现在可以处理参数为表达式和参数为矩阵且尺寸大小不正确的情况。在这种情况下,调整表达式的大小不会有任何问题,除非它们确实需要调整。这意味着,传递一个尺寸不正确的表达式将导致运行时错误(仅在调试模式下),而传递尺寸正确的表达式将正常工作。

注意:在上面的讨论中,术语 MatrixArray 以及 MatrixBaseArrayBase 可以互换,并且所有参数仍然有效。

总结

  • 总之,采用不可写(const 引用)对象作为函数参数的实现并不是一个大问题,也不会在编译和运行程序方面导致问题。然而,一个简单的实现可能会在代码中引入不必要的临时对象。为了避免将参数转换为临时对象,将它们作为(const)引用传递给 MatrixBaseArrayBase(因此将函数模板化)。
  • 需要接受可写的(非 const)参数的函数必须采用 const 引用,并在函数体内强制转换 const 属性。
  • 接受 MatrixBase(或 ArrayBase)对象作为参数的函数,如果这些对象可以调整大小(即重置大小),则必须在派生类上调用 resize()函数,由 derived() 返回。
相关文章
|
存储 编译器
[Eigen中文文档] 深入了解 Eigen - 类层次结构
本页面介绍了Eigen类层次结构中 Core 类的设计及其相互关系。一般用户可能不需要关注这些细节,但对于高级用户和Eigen开发人员可能会有用。
671 0
|
编译器 索引
[Eigen中文文档] 块操作
本文介绍了块操作。块是matrix或array的部分矩形元素。块表达式既可以用作右值也可以用作左值。与Eigen表达式一样,如果让编译器进行优化,则块操作的运行时间成本为零。
498 0
|
存储 算法 NoSQL
[Eigen中文文档] 存储顺序
矩阵和二维数组有两种不同的存储顺序:列优先和行优先。本节解释了这些存储顺序以及如何指定应该使用哪一种。
716 0
|
存储 算法 NoSQL
[Eigen中文文档] 稀疏矩阵操作
在许多应用中(例如,有限元方法),通常要处理非常大的矩阵,其中只有少数系数不为零。在这种情况下,可以通过使用仅存储非零系数的特殊表示来减少内存消耗并提高性能。这样的矩阵称为稀疏矩阵。
1065 0
|
存储 C语言 C++
|
存储 编译器
|
存储 缓存
[Eigen中文文档] 深入了解 Eigen - 惰性求值与混叠(Aliasing)
Eigen具有智能的编译时机制,可以实现惰性求值并在适当的情况下删除临时变量。它会自动处理大多数情况下的混叠问题,例如矩阵乘积。自动行为可以通过使用MatrixBase::eval()和MatrixBase::noalias()方法手动覆盖。
710 0
|
存储 安全 编译器
[Eigen中文文档] 常见的陷阱
本文将介绍一些Eigen常见的陷阱
844 0
|
安全 编译器 C++
[Eigen中文文档] 矩阵与向量运算
本文章旨在提供有关如何使用 Eigen 在矩阵、向量和标量之间执行算术操作的概述和一些详细信息。
1076 0
|
编译器 API 索引
[Eigen中文文档] 切片和索引
本文介绍了如何使用操作运算符operator()索引行和列的子集。该 API 在 Eigen 3.4 中引入。它支持 block API 提供的所有功能。特别是,它支持切片,即获取一组行、列或元素,以及等间隔的从矩阵或者数组中提取元素。
865 0