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

文档总目录

[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;
}

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

注意：在上面的讨论中，术语 Matrix 和 Array 以及 MatrixBase 和 ArrayBase 可以互换，并且所有参数仍然有效。

总结

• 总之，采用不可写（const 引用）对象作为函数参数的实现并不是一个大问题，也不会在编译和运行程序方面导致问题。然而，一个简单的实现可能会在代码中引入不必要的临时对象。为了避免将参数转换为临时对象，将它们作为（const）引用传递给 MatrixBase 或 ArrayBase（因此将函数模板化）。
• 需要接受可写的（非 const）参数的函数必须采用 const 引用，并在函数体内强制转换 const 属性。
• 接受 MatrixBase（或 ArrayBase）对象作为参数的函数，如果这些对象可以调整大小（即重置大小），则必须在派生类上调用 resize()函数，由 derived() 返回。