Eigen 中的模板化协方差函数

Question

我的一个项目需要一个模板化协方差函数，即将 MatrixXd、ArrayXXf 和/或它们的 .block() 作为输入并返回一个表达式以用于进一步的计算。

q1。作为验证，我在下面的尝试似乎有效，但它实际上返回了一个 Eigen 表达式吗？ （编辑：@Darhuuk 确认没有；幸运的是，C++14 或更高版本不是问题！）

q2。显然，我需要一些“内部矩阵”，而不是我在 Eigen 文档中遇到的关于模板化函数的内部 RowVector。因此我的问题是，应该如何创建一个内部 MatrixType z = x.rowwise() - x.colwise().mean()，以便函数可以返回 (z.transpose()*z)/(x.rows() -1)?内部行向量用于 Eigen 手册中的协方差示例（https://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html）

第三季度。最后，当将模板化函数添加到带有头文件的类时，我如何找到需要使用哪些 Eigen 类型来显式实例化模板，例如MatrixXd、ArrayXXd 的 .block() 还是 Map？我能找到的只是使用简单数据类型的示例（例如Storing C++ template function definitions in a .CPP file）

template <typename Derived>
Matrix<typename Derived::Scalar, Derived::ColsAtCompileTime, Derived::ColsAtCompileTime>  
SampleCov( const DenseBase<Derived>& x )
{

  typedef typename internal::plain_row_type<Derived>::type RowVectorType;
  const RowVectorType x_mean = x.colwise().mean();

  return ((x.rowwise() - x_mean).matrix().transpose() * (x.rowwise() - x_mean).matrix()) / (x.rows()-1);
}

Answer 1

您的函数显式返回 Matrix 对象，所以不，它不返回表达式对象。解决此问题的一种方法是让编译器为您推断返回类型（因为表达式对象将是一些巨大的模板化混乱）：

template <typename Derived>
auto SampleCov (DenseBase<Derived> const & x) {
  auto const x_mean = x.colwise().mean();
  return ((x.rowwise() - x_mean).matrix().transpose()
    * (x.rowwise() - x_mean).matrix()) / (x.rows() - 1);
}

这是假设您使用的是 C++14 或更高版本。对于仅使用 auto 作为返回类型的 C++11，您需要一个尾随返回类型：

template <typename Derived>
auto SampleCov (DenseBase<Derived> const & x)
    -> decltype(((x.rowwise() - x.colwise().mean()).matrix().transpose()
      * (x.rowwise() - x.colwise().mean()).matrix()) / (x.rows() - 1)) {
  auto const x_mean = x.colwise().mean();
  return ((x.rowwise() - x_mean).matrix().transpose()
    * (x.rowwise() - x_mean).matrix()) / (x.rows() - 1);
}

请注意，我删除了RowVectorType typedef，因为我没有看到重点。

关于问题 2，我认为上述示例解决了该问题，因为没有更明确命名的类型。返回类型和函数内部的 auto 负责这一点。

最后，问题 3。这取决于您到底想做什么。根据您的说法，上述函数似乎不适用于MatrixXd 对象或调用MatrixXd::block() 返回的对象。

对我来说，这并不表示您需要explicit template instantiation，这就是您要问的。

相反，您可以简单地使SampleCov 的参数类型更通用：

template <typename T>
auto SampleCovGeneric (T const & x) {
  auto const x_mean = x.colwise().mean();
  return ((x.rowwise() - x_mean).matrix().transpose()
    * (x.rowwise() - x_mean).matrix()) / (x.rows() - 1);
}

只要您可以在 T 类型的对象上调用 colwise()、matrix()、...，就可以了。无论如何，这可能是一种更好的方法，因为现在可以将 Eigen 表达式传递给SampleCovGeneric。例如：

Eigen::MatrixXd a(2,2);

// aTranspose is NOT a matrix object!
// Its type is Eigen::Transpose<Eigen::Matrix<double, Dynamic, Dynamic>>.
auto aTranspose = a.transpose();

/* Note use of auto and eval() calls!
 * See https://eigen.tuxfamily.org/dox/TopicPitfalls.html.
 */
auto b = SampleCov(aTranspose.eval()).eval();
auto c = SampleCovGeneric(aTranspose).eval();

在上面的例子中，SampleCov 需要一个 DenseBase<Derived> 类型的对象。但是aTranspose不是这样的类型。所以我们必须首先明确评估（即实际计算）转置。根据SampleCov 内部发生的情况，这是一个无用的计算。例如。想象计算的第一件事是它的转置参数 (x)。在这种情况下，那将只是原来的 a 。但是如果转置已经被评估，那么 Eigen 就无法“看到”它，需要计算转置的转置。

另一方面，SampleCovGeneric 接受任何类型，因此无需先评估aTranspose。这（可能）允许 Eigen “看穿”一堆计算，从而以更优化的方式计算结果。参见例如https://eigen.tuxfamily.org/dox/TopicLazyEvaluation.html.

如果你真的需要，你可以explicitly instantiate templates，这样它们就可以被分割成头文件和源文件。会是这样的：

/* Header */
#pragma once

#include <utility>

// Reduce trailing return type mess a little
template <class T>
using SampleCovResultType = decltype(((std::declval<T const &>().rowwise()
      - std::declval<T const &>().colwise().mean()).matrix().transpose()
    * (std::declval<T const &>().rowwise()
      - std::declval<T const &>().colwise().mean()).matrix())
  / (std::declval<T const &>().rows() - 1));

template <typename T>
auto SampleCov (T const & x) -> SampleCovResultType<T>;

// Explicit template declaration
extern template
auto SampleCov<Eigen::MatrixXd> (Eigen::MatrixXd const & x)
    -> SampleCovResultType<Eigen::MatrixXd>;

/* Source */
#include "SampleCov.h"

template <class T>
auto SampleCov (T const & x) -> SampleCovResultType<T> {
  auto const x_mean = x.colwise().mean();
  return ((x.rowwise() - x_mean).matrix().transpose()
    * (x.rowwise() - x_mean).matrix()) / (x.rows() - 1);
}

// Explicit template definition
template
auto SampleCov<Eigen::MatrixXd> (Eigen::MatrixXd const & x)
    -> SampleCovResultType<Eigen::MatrixXd>;

我认为你不能只使用 auto 作为返回类型。由于函数定义放在单独的源文件中，编译器看不到函数体，所以在调用函数时无法推断返回类型。所以你必须使用尾随返回类型。

请注意以下几点：

显式实例化声明（外部模板）可防止 隐式实例化：否则会导致 隐式实例化必须使用显式实例化 在程序的其他地方提供的定义。

即如果您尝试使用未显式实例化 SampleCov 的类型 T 调用 SampleCov，编译将失败。

注意：以上所有代码均未经测试。拼写错误的可能性很高:)。

编辑：

• 修复缺少的模板参数。
• 删除了关于调用matrix() 的备注。我想知道它是否评估了您不希望它执行的表达式。 AFAIK，它没有，基于例如ArrayBase's documentation。
• 添加了有关显式实例化的信息。