用数组进行特征向量化

Question

我正在处理点云数据（每个云 150k 点）。我想，对于每个 (x,y) 点，计算到参考点 O 的距离和方位角：

for each point p in points
    dx = p.x - ox
    dy = p.y - oy
    d = hypot(dx, dy)
    az = atan2(dy, dx)

我有一个手动 SSE 实施。我希望使用 eigen 使代码更清晰：

ArrayXf x(points.size()), y(points.size());
for(unsigned i=0; i<points.size(); ++i) {
    x[i] = points[i].x;
    y[i] = points[i].y;
}
const ArrayXf d = (dx.square() + dy.square()).sqrt();
// implement a polynomial approximation to atan (same as the SSE)

但是，从我的计时实验来看，这似乎根本没有矢量化，因为时间与基线实现相同。而且我知道 SSE2 已启用，因为我在同一个文件中编译了一些 SSE2 代码。

但是，根据文档，Eigen 在支持时确实利用了 SSE2（以及 3.3 中的 AVX）。是不是只针对向量和矩阵运算？

编辑：我研究了生成的汇编代码，它确实包含一些 SSE 指令。但是还是很慢

编辑：这里有更多时间信息。我正在循环超过 100 帧，每帧大约 150k 点。

• 没有 atan2 的简单实现：150 毫秒
• sse 实现（处理 4 x 4 点并丢弃最后几个未填充完整数据包的点）：30 毫秒
• 使用特征图的特征实现：90 毫秒（差异：36 毫秒，假设：16 毫秒，索引：17 毫秒）

这是我的特征码：

const Eigen::Map<const Eigen::ArrayXf, Eigen::Unaligned, Eigen::InnerStride<4> > px(&(points[0].x), points.size());
const Eigen::Map<const Eigen::ArrayXf, Eigen::Unaligned, Eigen::InnerStride<4> > py(&(points[0].y), points.size());

// difference with the origin (ox and oy are floats)
const Eigen::ArrayXf dx = px - ox, dy = py - oy;

// distance and index
const Eigen::ArrayXf d = sqrt(dx.square() + dy.square());

static const float r_res_mult = 1.0f / r_res; //2x faster than div
const Eigen::ArrayXi didx = (d * r_res_mult).cast<int>();

Answer 1

您的主要问题是您的数据格式不适合 SIMD。您正在使用结构数组 (xyxyxyxyxyxy...)，然后对代码进行矢量化处理

for(unsigned i=0; i<points.size(); ++i) {
    x[i] = points[i].x;
    y[i] = points[i].y;
}

转换为数组结构（xxxxxxxx....yyyyyyy...）。这种转换成本很高。

更好的解决方案是已经将您的点存储为数组结构。一个更好的解决方案是使用数组的混合结构，也就是数组结构的数组。对于 SSE，假设您使用的是单浮点，那么您将执行 xxxxyyyyxxxxyyyy....

接下来我建议您使用 SIMD 数学库。英特尔提供昂贵且封闭源代码的 SVML。 AMD 提供免费但封闭源代码的libm。 但是这些库在竞争对手的硬件上都不能很好地发挥作用。 最好的 SIMD 库是 Agner Fog 的 Vector Class Library (VCL) 。它是开源的、免费的，并且针对英特尔和 AMD 处理器进行了优化。它也与 Eigen 一样，只是头文件，因此，与 Eigen 一样，您不必编译和链接库。您刚刚包含了头文件。以下是 SSE 或 AVX 的浮动方式（VLC 将在没有 AVX 的系统上模拟 AVX）。

//    g++ -O3 -Ivectorclass -msse4.2 foo.cpp
// or g++ -O3 -Ivectorclass -mavx foo.cpp
#include <vectorclass.h>
#include <vectormath_trig.h>

struct Point2DBlock {
    float x[8];
    float y[8];
};

int main(void) {
    const int nblocks = 10; //each block contains eight points
    Point2DBlock aosoa[nblocks]; //xxxxxxxxyyyyyyyy xxxxxxxxyyyyyyyy ...
    float ox = 0.0f, oy = 0.0f;
    Vec8f vox = ox, voy = oy;
    for(int i=0; i<nblocks; i++) {
        Vec8f dx = Vec8f().load(aosoa[i].x) - vox;
        Vec8f dy = Vec8f().load(aosoa[i].y) - voy;
        Vec8f d  = sqrt(dx*dx + dy*dy);
        Vec8f az = atan2(dy,dx);
    } 
}

如果你真的需要hypot。您可以使用 pseudo-code from wikipedia 从 VCL 构建一个。

static inline Vec8f hypot(Vec8f const &x, Vec8f const &y) {
    Vec8f t;
    Vec8f ax = abs(x), ay = abs(y);
    t  = min(ax,ay);
    ax = max(ax,ay);
    t  = t/ax;
    return ax*sqrt(1+t*t);
}

编辑：

这是一个使用结构数组的方法。这需要一些改组，但与其他计算相比，这可能可以忽略不计。 VLC 使用模板元编程来确定一种有效的洗牌方法。

#include <vectorclass.h>
#include <vectormath_trig.h>

int main(void) {
    const int npoints=80;
    float points[2*npoints]; //xyxyxyxyxyxy...
    float ox = 0.0, oy = 0.0;
    Vec8f vox = ox, voy = oy;
    for(int i=0; i<npoints; i+=16) {
        Vec8f l1 = Vec8f().load(&points[i+0]);
        Vec8f l2 = Vec8f().load(&points[i+8]);
        Vec8f dx = blend8f<0, 2, 4, 6, 8, 10, 12, 14>(l1,l2) - vox;
        Vec8f dy = blend8f<1, 3, 5, 7, 9, 11, 13, 15>(l1,l2) - voy;
        Vec8f d  = sqrt(dx*dx + dy*dy);
        Vec8f az = atan2(dy,dx);
    } 
}

Answer 2

复制需要很长时间。与计算本身相同或更长。您不必像那样复制数据。它很冗长，可能更慢。您可以改用地图，甚至可以直接使用地图进行计算。我写了一个快速演示：

int sz = 15000000;
std::vector<Point> points(sz);

Eigen::Map<ArrayXd, Unaligned, InnerStride<2>> mapX(&(points[0].x), sz);
Eigen::Map<ArrayXd, Unaligned, InnerStride<2>> mapY(&(points[0].y), sz);

mapX = ArrayXd::Random(sz);
mapY = ArrayXd::Random(sz);

auto cpstart = std::chrono::high_resolution_clock::now();
ArrayXd x = mapX;
ArrayXd y = mapY;
ArrayXd d;
auto cpend = std::chrono::high_resolution_clock::now();

auto mpSumstart = std::chrono::high_resolution_clock::now();

d = (mapX.square() + mapY.square()).sqrt().eval();

auto mpSumend = std::chrono::high_resolution_clock::now();

std::cout << d.mean() << "\n";

auto arStart = std::chrono::high_resolution_clock::now();

d = (x.square() + y.square()).sqrt().eval();

auto arEnd = std::chrono::high_resolution_clock::now();

std::cout << d.mean() << "\n";

auto elapsed = cpend - cpstart;
std::cout << "Copy: " <<  elapsed.count() << '\n';
std::cout << "Map: " <<  (mpSumend - mpSumstart).count() << '\n';
std::cout << "Array: " <<  (arEnd - arStart).count() << '\n';

我得到的时间是你的数组长度的 100 倍，我只是懒得写一个循环来更好地测试。复制在我的系统上大约需要 90ms（VS2012 /Ox in Release (-DNDEBUG)），映射版本需要 185ms，复制的数组也需要大约 90ms。 2 的近似因子对于 SIMD 操作是有意义的，因为映射的版本每隔一个 double 会跳过。如果你有一个数组结构而不是结构数组，那么映射的性能应该与复制数组的性能相当。

编辑：我定义了EIGEN_DONT_VECTORIZE，并且复制的数组（几乎）加倍了它的时间（如预期的那样）。然而，地图保持不变。好奇的。可能与未对齐的地图有关。或者只是因为只有两个双打的空间，而其他所有的双打都属于错误的地图。

EDIT 2 关于问题中提出的具体问题，我想到了一个愚蠢的想法。您可以将x、y 值视为std::complex<double>，然后将其作为没有内存副本的单个块加载：

Eigen::Map<ArrayXcd> mapC((std::complex<double>*)(&(points[0].x)), sz);
//...
cd = mapC.cwiseAbs2().sqrt().eval();

时间只比我电脑上预先复制的数组稍长。您还可以将原点作为复数减去

cd = (mapC - std::complex<double>(ox, oy)).cwiseAbs2().sqrt().eval();