Matlab repelem 的 Eigen 等效项是啥？

Posted 2023-02-16

【中文标题】Matlab repelem 的 Eigen 等效项是啥？

【英文标题】：What is the Eigen equivalent of Matlab repelem?Matlab repelem 的 Eigen 等效项是什么？

【发布时间】：2018-10-23 09:46:35

【问题描述】：

在 Matlab 中，有repmat 和repelem。对于大小为[d, n1] 和[d, n2] 的矩阵x1 和x2，您可以在Matlab 中执行此操作：

n1 = size(x1, 2);
n2 = size(x2, 2);

M = repmat(x1, 1, n2) - repelem(x2, 1, n1);

什么是等效的特征码？以下是我不太满意的四种变体。我想知道它是否可以制成更快的单线？

TL;DR：变体 2 是最好的，但这取决于编译器和其他因素。

int d = x1.rows();
int n1 = x1.cols();
int n2 = x2.cols();

Eigen::MatrixXd M(d, n1*n2);

// Variant 1:
int idx = 0;
for (int c = 0; c != n2; c++) 
  for (int r = 0; r != n1; r++) 
    M.col(idx) = x1.col(r) - x2.col(c);
    idx += 1;
  


// Variant 2:
for (int c = 0, idx = 0; c != n2; c += 1, idx += n1)
  M.block(0, idx, d, n1) = x1.colwise() - x2.col(c);

// Variant 3:
M = - x2.replicate(n1, 1);
M.resize(d, n1*n2);
M += x1.replicate(1, n2);

// Variant 5:
M = x1.rowwise().replicate(n2) - x2(all,VectorXi::LinSpaced(n1*n2,0,n2-1));

这是我目前的时间安排，完整代码见下文。 Matlab 时序适用于多线程和单线程 Matlab R2017b。使用标志 -O3 -DNDEBUG -march=native -mtune=native 编译的 C++ 版本。所有都在i5-6500 上运行，所以我有AVX。

以秒为单位的时间 代码 gcc-7 gcc-8 clang-6 ---------------------------------- Matlab mt 51 Matlab st 84 五、1 38 37 57 五、2 36 34 23 五、3 598 599 187 五、5 94 172 107

Matlab 代码：

ds = 1:10;
n1s = 5:5:500;
n2s = 5:5:500;

z1 = randn(max(ds), max(n1s));
z2 = randn(max(ds), max(n2s));

tic;
s = 0;
for idx = 1:numel(ds)
    for jdx = 1:numel(n1s)
        for kdx = 1:numel(n2s)
            K = MFdiff(z1(1:ds(idx), 1:n1s(jdx)),...
                z2(1:ds(idx), 1:n2s(kdx)));
            s = s + K(1,1);
        end
    end
end
toc

C++ 代码：

#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;

template <typename Derived1, typename Derived2>
MatrixXd MFdiff1(
  DenseBase<Derived1> const & x1,
  DenseBase<Derived2> const & x2)

  int d = x1.rows();
  int n1 = x1.cols();
  int n2 = x2.cols();

  MatrixXd out(d, n1*n2);

  int idx = 0;
  for (int c = 0; c != n2; c++) 
    for (int r = 0; r != n1; r++) 
      out.col(idx) = x1.col(r) - x2.col(c);
      idx += 1;
    
  

  return out;


template <typename Derived1, typename Derived2>
MatrixXd MFdiff2(
  DenseBase<Derived1> const & x1,
  DenseBase<Derived2> const & x2)

  int d = x1.rows();
  int n1 = x1.cols();
  int n2 = x2.cols();

  MatrixXd out(d, n1*n2);

  for (int c = 0, idx = 0; c != n2; c+=1, idx += n1)
    out.block(0, idx, d, n1) = x1.colwise() - x2.col(c);

  return out;


template <typename Derived1, typename Derived2>
MatrixXd MFdiff3(
  DenseBase<Derived1> const & x1,
  DenseBase<Derived2> const & x2)

  int d = x1.rows();
  int n1 = x1.cols();
  int n2 = x2.cols();

  MatrixXd out;
  out = - x2.replicate(n1, 1);
  out.resize(d, n1*n2);
  out += x1.replicate(1, n2);

  return out;


template <typename Derived1, typename Derived2>
MatrixXd MFdiff5(
  DenseBase<Derived1> const & x1,
  DenseBase<Derived2> const & x2)

  int n1 = x1.cols();
  int n2 = x2.cols();

  return x1.rowwise().replicate(n2) - x2(all,VectorXi::LinSpaced(n1*n2,0,n2-1));


double test(VectorXi const & ds,
      VectorXi const & n1s,
      VectorXi const & n2s)

  MatrixXd z1 = MatrixXd::Random(ds.maxCoeff(), n1s.maxCoeff());
  MatrixXd z2 = MatrixXd::Random(ds.maxCoeff(), n2s.maxCoeff());

  double s = 0;
  for (int idx = 0; idx!=ds.rows(); idx++) 
    for (int jdx = 0; jdx!=n1s.rows(); jdx++) 
      for (int kdx = 0; kdx!=n2s.rows(); kdx++) 
        MatrixXd K = MFdiff5(z1.block(0, 0, ds(idx), n1s(jdx)),
                 z2.block(0, 0, ds(idx), n2s(kdx)));
        s += K(0,0);
      
    
  

  return s;
  

int main() 
  VectorXi ds = VectorXi::LinSpaced(10, 1, 10);
  VectorXi n1s = VectorXi::LinSpaced(100, 5, 500);
  VectorXi n2s = VectorXi::LinSpaced(100, 5, 500);

  std::cout << test(ds, n1s, n2s) << '\n';

【问题讨论】：

this 或 this 是否回答您的问题？

@AviGinsburg 不，但第二个问题启发了上面的变体 3，实际上速度较慢。

你能发布时间和全套编译器标志吗？另外，你说的是什么尺寸？并尝试像 M = x1.replicate(1, n2) - x2.replicate(n1, 1); 这样的 Eigen 版本。中间的调整大小只会破坏性能。

@AviGinsburg 声明 M = x1.replicate(1, n2) - x2.replicate(n1, 1); 只比上面的 V3 稍微快一点，但它计算不正确！ （试试看！）我的编译器标志是g++-8 -std=c++1z -Wall -Wextra -pedantic -O3 -march=native -mtune=native，但预期的问题更多的是算法类型。

大小如下：d在1-10范围内，n1'和n2在1-1000范围内。

【参考方案1】：

用 Eigen 的头你可以写：

M = x1.rowwise().replicate(n2) - x2(Eigen::all,VectorXi::LinSpaced(n1*n2,0,n2-1));

与您的变体 2 具有相同的速度。

独立的基准测试（BenchTimer.h 需要一个 repo 的克隆），使用 gcc 7 和 clang 6 使用 -O3 -DNDEBUG 进行测试：

#include <iostream>
#include <Eigen/Dense>
#include <bench/BenchTimer.h>
using namespace Eigen;
using namespace std;

EIGEN_DONT_INLINE
void foo1(const MatrixXd& x1, const MatrixXd& x2, MatrixXd& M)

  int d = x1.rows();
  int n1 = x1.cols();
  int n2 = x2.cols();

  int idx = 0;
  for (int c = 0; c != n2; c++) 
    M.block(0, idx, d, n1) = x1.colwise() - x2.col(c);
    idx += n1;
  


EIGEN_DONT_INLINE
void foo2(const MatrixXd& x1, const MatrixXd& x2, MatrixXd& M)

  int n1 = x1.cols();
  int n2 = x2.cols();

  M = x1.rowwise().replicate(n2) - x2(all,VectorXi::LinSpaced(n1*n2,0,n2-1));


int main()

  int tries = 2;
  int rep = 1;

  int d = 100;
  int n1 = 100;
  int n2 = 100;

  MatrixXd x1(d,n1); x1.setRandom();
  MatrixXd x2(d,n2); x2.setRandom();
  MatrixXd M(d, n1*n2);

  BenchTimer t;
  BENCH(t, tries, rep, foo1(x1, x2, M));
  std::cout << "Time: " << t.best() << "s" << std::endl;

  BENCH(t, tries, rep, foo2(x1, x2, M));
  std::cout << "Time: " << t.best() << "s" << std::endl;

【讨论】：

从马的嘴里。奇怪的是，当我运行它时，它比 v.2 慢了很多。我必须制定一个可靠的基准，然后再回到这个问题。

我在我的问题中将其添加为变体 5。在我的基准测试中它要慢得多。我还在我的机器上运行了你的代码，得到了这个：时间：0.00131418s 时间：0.0111729s

hm... 在 2.6GHz Haswell 上，使用 clang 6 和 -O3 -DNDEBUG，我得到 = Time: 0.000690089s Time: 0.000568082s。 gcc 5、6 或 7 给了我类似的结果。

是的，我之前报告的时间是异常的，运行它很多次我看到大约 1:1 在 clang 6 上，1:1,5 在 gcc 7 上和 1:8 在 gcc 8 上（在赞成第 2 节）。似乎依赖于编译器。在问题中查看我的基准。

我看到在您的基准测试中 d 非常小，介于 1 到 10 之间。在这种情况下，我同意 v2 更快。