Eigen::Tensor和Eigen::Matrix性能比较

Question

我想用一个3D的Eigen::Tensor代替我的代码中一系列的矩阵。为此，我试图比较Tensor和Matrix的性能。

下面的"tensorContractTest"函数对一个(n,n,n)三阶张量与大小为n的一阶张量进行收缩（n = 500）。这个收缩计算n ** 2个点积，因此在操作数方面，它等价于两个(n, n)矩阵的乘法（下面的"matrixProductTest"函数）。

当在Visual Studio 2013上运行时，"tensorContractTest"函数的运行速度比"matrixProductTest"慢了约40倍。 可能是我漏掉了什么。谢谢帮助。

#include <unsupported/Eigen/CXX11/Tensor>
using namespace Eigen;

// Contracts 3-dimensional (n x n x n) tensor with 1-dimensional (n) tensor. 
// By the number of operations, it's equivalent to multiplication of 
// two (n, n) matrices (matrixProdTest).
Tensor<double, 2>  tensorContractTest(int n)
{
  Tensor<double, 3> a(n, n, n);     a.setConstant(1.);  
  Tensor<double, 1> b(n);           b.setConstant(1.);
  auto indexPair            = array<IndexPair<int>, 1>{IndexPair<int>(2,0)}; 
  Tensor<double, 2> result  = a.contract(b, indexPair); 
  return result;
}

MatrixXd  matrixProductTest(int n)
{ 
  MatrixXd a = MatrixXd::Ones(n, n), result = a * a;
  return result;
}

- lbo

---

你启用了优化吗？请使用适当的基准测试工具，例如Google Benchmark，而不是仅仅随意给出数字。 - Henri Menke

我使用了/O2编译器标志。我的计时函数如下。#include <chrono></br> #include <utility><br/> #include <iostream><br>// 计时函数<br> typedef std::chrono::high_resolution_clock::time_point TimeVar;<br> #define duration(a) std::chrono::duration_caststd::chrono::nanoseconds(a).count() #define timeNow() std::chrono::high_resolution_clock::now()template<typename F, typename... Args> double funcTime(F func, Args&&... args) { TimeVar t1 = timeNow(); func(std::forward<Args>(args)...); return duration(timeNow() - t1); } - lbo

@HenriMenke。我使用了/O2编译器标志。 - lbo

Answer 1

尽管浮点运算的数量是同一级别的，但内存访问模式完全不同，因此这两个操作根本无法进行比较。一般来说，矩阵-矩阵操作将始终比矩阵-向量或向量-向量操作更快（以FLOPS为单位），因为前者可以更好地使用缓存，从而实现CPU ALU的近乎最优使用。在您的情况下，一方面您需要读取一个n^3张量与两个n^2矩阵进行比较，因此内存占用量根本无法进行比较。

在内部，Tensor::contract 在可能时会回退到Eigen的矩阵乘积内核，因此性能应该相当。

- ggael

Answer 2

因为张量缩并与矩阵乘法不同。

存在着专门用于矩阵乘法的算法，例如 Strassen 算法可以减少总操作数。此外，经过多年高度优化的矩阵库通常专门使用矢量化指令（SIMD 或 AVX），具体取决于平台（Intel、AMD、ARM）。对于小尺寸或稀疏模式的矩阵，与非专用代码相比，速度提升巨大。

相比之下，张量库往往不太优化。因此，如果您可以将张量数学转换为矩阵代数，则很有可能会提高速度。

- Mauricio Cele Lopez Belon

Answer 3

我可以使用Google基准测试程序重现缓慢的情况。经过测试，如果数据大小较小，Eigen的矩阵乘法比Eigen的张量收缩更快。尽管结果相同。

所以我不确定回退方案是什么

在内部，Tensor::contract在可能的情况下会回退到Eigen的矩阵乘积核心，因此性能应该是相当的。

@ggael 在这里提到了。

---

这是我的测试代码：

#include <Eigen/Dense>
#include <benchmark/benchmark.h>
#include <unsupported/Eigen/CXX11/Tensor>

constexpr std::size_t I = ;
constexpr std::size_t J = ;
constexpr std::size_t K = ;

static void bmForLoop(benchmark::State &state) {
  // Perform setup here
  Eigen::Matrix<double, I, J> m_0 = Eigen::Matrix<double, I, J>::Random();
  Eigen::Matrix<double, J, K> m_1 = Eigen::Matrix<double, J, K>::Random();
  Eigen::Matrix<double, I, K> m_2;

  for (auto _ : state) {
    // This code gets timed
    m_2 = Eigen::Matrix<double, I, K>::Zero();
    for (std::size_t i = 0; i < I; ++i) {
      for (std::size_t k = 0; k < K; ++k) {
        for (std::size_t j = 0; j < J; ++j) {
          benchmark::DoNotOptimize(m_2.coeffRef(i, k) +=
                                   m_0.coeff(i, j) * m_1.coeff(j, k));
        }
      }
    }
  }
}

static void bmEigenMatrix(benchmark::State &state) {
  // Perform setup here
  Eigen::Matrix<double, I, J> m_0 = Eigen::Matrix<double, I, J>::Random();
  Eigen::Matrix<double, J, K> m_1 = Eigen::Matrix<double, J, K>::Random();
  Eigen::Matrix<double, I, K> m_2;

  for (auto _ : state) {
    // This code gets timed
    benchmark::DoNotOptimize(m_2 = m_0 * m_1);
  }
}

static void bmEigenMatrixDynamic(benchmark::State &state) {
  // Perform setup here
  Eigen::MatrixXd m_0(I, J);
  m_0 = Eigen::Matrix<double, I, J>::Random();
  Eigen::MatrixXd m_1(J, K);
  m_1 = Eigen::Matrix<double, J, K>::Random();
  Eigen::MatrixXd m_2;

  for (auto _ : state) {
    // This code gets timed
    benchmark::DoNotOptimize(m_2 = m_0 * m_1);
  }
}

static void bmEigenTensorFixed(benchmark::State &state) {
  // Perform setup here
  Eigen::Matrix<double, I, J> m_0 = Eigen::Matrix<double, I, J>::Random();
  Eigen::Matrix<double, J, K> m_1 = Eigen::Matrix<double, J, K>::Random();
  Eigen::TensorMap<Eigen::TensorFixedSize<double, Eigen::Sizes<I, J>>> t_0(
      m_0.data(), I, J);
  Eigen::TensorMap<Eigen::TensorFixedSize<double, Eigen::Sizes<J, K>>> t_1(
      m_1.data(), J, K);
  Eigen::TensorFixedSize<double, Eigen::Sizes<I, K>> t_2;

  for (auto _ : state) {
    // This code gets timed
    benchmark::DoNotOptimize(
        t_2 = t_0.contract(
            t_1, Eigen::array<Eigen::IndexPair<Eigen::Index>, 1>{{{1, 0}}}));
  }
}

static void bmEigenTensor(benchmark::State &state) {
  // Perform setup here
  Eigen::Matrix<double, I, J> m_0 = Eigen::Matrix<double, I, J>::Random();
  Eigen::Matrix<double, J, K> m_1 = Eigen::Matrix<double, J, K>::Random();
  Eigen::TensorMap<Eigen::Tensor<double, 2>> t_0(m_0.data(), I, J);
  Eigen::TensorMap<Eigen::Tensor<double, 2>> t_1(m_1.data(), J, K);
  Eigen::Tensor<double, 2> t_2;

  for (auto _ : state) {
    // This code gets timed
    benchmark::DoNotOptimize(
        t_2 = t_0.contract(
            t_1, Eigen::array<Eigen::IndexPair<Eigen::Index>, 1>{{{1, 0}}}));
  }
}

BENCHMARK(bmForLoop);
BENCHMARK(bmEigenMatrix);
BENCHMARK(bmEigenMatrixDynamic);
BENCHMARK(bmEigenTensorFixed);
BENCHMARK(bmEigenTensor);
BENCHMARK_MAIN();

其中，我、J和K

constexpr std::size_t I = ;
constexpr std::size_t J = ;
constexpr std::size_t K = ;

这些袋子里填充着从1到8的相同数字。

我使用Bazel进行编译，编译命令是build --cxxopt=-std=c++17 --cxxopt=-march=native --cxxopt=-O3 --cxxopt=-DNDEBUG。Eigen的版本是3.4.0。

- Yukun Xia