Eigen 中的模板化协方差函数

Question

我的一个项目需要一个模板化的协方差函数，即将 MatrixXd、ArrayXXf and/or a .block() 作为输入并 returning 一个表达式用于进一步的计算。

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

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

q3。最后，当使用头文件将模板函数添加到 class 时，我如何找到需要使用哪些 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

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

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 作为 return 类型的 C++11，您需要一个尾随 return 类型：

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，我认为上述示例解决了该问题，因为没有更明确命名的类型。 return 类型和函数内部的 auto 都会处理这个问题。

最后，问题3。这取决于你到底想做什么。根据您所说的，上述功能似乎不适用于 MatrixXd 对象或通过调用 MatrixXd::block().

编辑的对象 return

对我来说，这并不真正表明您需要 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 不能 "see" 并且需要计算转置的转置。

另一方面，

SampleCovGeneric 接受任何类型，因此无需先评估 aTranspose。这（可能）允许 Eigen "see through" 一堆计算，从而以更优化的方式计算结果。参见例如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 作为 return 类型。由于函数定义放在单独的源文件中，编译器看不到函数体，所以调用函数时无法推导出return类型。所以你必须使用尾随 return 类型。

不过请注意以下几点：

An explicit instantiation declaration (an extern template) prevents implicit instantiations: the code that would otherwise cause an implicit instantiation has to use the explicit instantiation definition provided somewhere else in the program.

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

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

编辑：

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

Eigen 中的模板化协方差函数

Templated covariance function in Eigen

c++

templates

eigen