c - 3D vector 加速欧氏距离

Question

我需要执行一个非常常见且简单的矩阵运算。
然而我需要它快，非常快......

我已经在考虑多线程实现，但现在我只想看看在单个处理器上实现它的速度有多快。

矩阵运算如下:

我正在计算点 vector (A) 和引用点 (B) 之间的欧几里德距离。
这些点位于 3D 空间中，每个点都有一组 X、Y 和 Z 坐标。
因此，点 vector 由三个 float 组来描述，其中保存每个点的 X、Y、Z 坐标。
输出是另一个长度为 N 的 vector ，保存数组中每个点与引用点之间的距离。

三个 XYZ 数组排列为 Nx3 矩阵的列。

x[0]      y[0]      z[0]
x[1]      y[1]      z[1]
x[2]      y[2]      z[2]
x[3]      y[3]      z[3]
.         .         .
.         .         .
.         .         .
x[N-1]    y[N-1]    z[N-1]

在内存中，矩阵按行主序排列为一维数组，依次包含 X、Y 和 Z 列的值。

x[0], x[1], x[2], x[3] . . . x[N-1], y[0], y[1], y[2], y[3] . . . y[N-1], z[0], z[1], z[2], z[3] . . . z[N-1]

整个事情有点复杂，因为我们需要在求平方根之前向矩阵的每个成员添加一个标量。

以下是简单 C 代码中的例程:

void calculateDistances3D(float *matrix, float Bx, float By, float Bz, float scalar, float *distances, int N)
{
    float *Ax = matrix;
    float *Ay = Ax + N;
    float *Az = Ay + N;
    int i;

    for (i = 0; i < N; i++) {

        float dx = Ax[i] - Bx;
        float dy = Ay[i] - By;
        float dz = Az[i] - Bz;

        float dx2 = dx * dx;
        float dy2 = dy * dy;
        float dz2 = dz * dz;

        float squaredDistance = dx2 + dy2 + dz2;
        float squaredDistancePlusScalar = squaredDistance + scalar;

        distances[i] = sqrt(squaredDistancePlusScalar);
    }
}

…这是简单的 Accelerate 实现(使用 vDSP 和 VecLib):
(请注意，所有处理都是就地执行的)

void calculateDistances3D_vDSP(float *matrix, float Bx, float By, float Bz, float scalar, float *distances, int N)
{
    float *Ax = matrix;
    float *Ay = Ax + N;
    float *Az = Ay + N;

    // for each point in the array take the difference with the reference point
    Bx = -Bx;
    By = -By;
    Bz = -Bz;
    vDSP_vsadd(Ax, 1, &Bx, Ax, 1, N);
    vDSP_vsadd(Ay, 1, &By, Ay, 1, N);
    vDSP_vsadd(Az, 1, &Bz, Az, 1, N);

    // square each coordinate
    vDSP_vsq(Ax, 1, Ax, 1, N);
    vDSP_vsq(Ay, 1, Ay, 1, N);
    vDSP_vsq(Az, 1, Az, 1, N);

    // reduce XYZ columns to a single column in Ax (reduction by summation)
    vDSP_vadd(Ax, 1, Ay, 1, Ax, 1, N);
    vDSP_vadd(Ax, 1, Az, 1, Ax, 1, N);

    // add scalar
    vDSP_vsadd(Ax, 1, &scalar, Ax, 1, N);

    // take sqrt
    vvsqrtf(distances, Ax, &N);
}

在 vDSP 库中，唯一可用于计算 vector 之间距离的函数是:

vDSP_vdist()
vDSP_distancesq()
vDSP_vpythg()

也许我遗漏了一些东西，但据我所知，它们都不支持计算 3D 距离所需的三个输入 vector 。

需要注意的几件事:
(1) 我不是在比较距离，所以我不能接受平方距离。我需要真实距离，因此计算平方根是绝对必要的。
(2) 如果您确实认为取倒数平方根可以使代码显着加快，那么这是一种可能。

我的印象是我没有充分利用 Accelerate 框架的潜力。
我正在寻找一些更智能、可能更简洁的东西，在更少的函数调用中完成更多的工作。以其他方式重新排列内存也可以，但我认为内存布局就这样就很好了。

我也愿意接受有关在英特尔处理器上工作的其他高度优化/矢量化线性代数库的建议。我不在乎它们是商业解决方案还是开源解决方案，只要它们的性能快速且强大即可。

问题是:Accelerate 框架中实现比上述更快的代码的最佳函数或函数组合是什么？

我正在运行 Mac OS X El Capitan 的 MacBook Pro(Retina，15 英寸，2014 年中)上的 Xcode 7 中进行开发。

谢谢。

Answer 1

试试这个。

• 在 iMac 和 iPhone 上重复 1000 次时，N = 2^20 的性能提高约 20%
• 此外，matrix 可以严格视为只读，因为仅操作距离
• 它在代数上等价，但在数值上不等价于您的实现；预计输出差异约为 10^-6

在我看来，vDSP 对于进一步“有针对性”的优化而言级别太高。相反，您可以查看 Ray Wenderlich’s iOS Assembly Tutorial作为使用 NEON 的起点，并针对此特定问题编写您自己的 SIMD 指令。

根据问题的大小 N，您还可以通过使用 GPU 获得进一步的加速，例如使用 Metal .

void calculateDistances3D_vDSP(float *matrix, float Bx, float By, float Bz, float scalar, float *distances, int N)
{
  float *Ax = matrix;
  float *Ay = Ax + N;
  float *Az = Ay + N;

  float constants = Bx*Bx + By*By + Bz*Bz + scalar;

  Bx = -2.0f*Bx;
  By = -2.0f*By;
  Bz = -2.0f*Bz;

  vDSP_vsq(Ax, 1, distances, 1, N);                      // Ax^2
  vDSP_vma(Ay, 1, Ay, 1, distances, 1, distances, 1, N); // Ax^2 + Ay^2
  vDSP_vma(Az, 1, Az, 1, distances, 1, distances, 1, N); // Ax^2 + Ay^2 + Az^2
  vDSP_vsma(Ax, 1, &Bx, distances, 1, distances, 1, N);  // Ax^2 + Ay^2 + Az^2 - 2*Bx
  vDSP_vsma(Ay, 1, &By, distances, 1, distances, 1, N);  // Ax^2 + Ay^2 + Az^2 - 2*Bx - 2*By
  vDSP_vsma(Az, 1, &Bz, distances, 1, distances, 1, N);  // Ax^2 + Ay^2 + Az^2 - 2*Bx - 2*By - 2*Bz
  vDSP_vsadd(distances, 1, &constants, distances, 1, N); // ... + constants = (Ax-Bx)^2 + (Ay-By)^2 + (Az-Bz)^2 + scalar

  /*
  vDSP_vsq(Ax, 1, distances, 1, N);                      // Ax^2
  vDSP_vsma(Ax, 1, &Bx, distances, 1, distances, 1, N);  // Ax^2 - 2*Ax*Bx
  vDSP_vma(Ay, 1, Ay, 1, distances, 1, distances, 1, N); // Ax^2 - 2*Ax*Bx + Ay^2
  vDSP_vsma(Ay, 1, &By, distances, 1, distances, 1, N);  // Ax^2 - 2*Ax*Bx + Ay^2 - 2*Ay*By
  vDSP_vma(Az, 1, Az, 1, distances, 1, distances, 1, N); // Ax^2 - 2*Ax*Bx + Ay^2 - 2*Ay*By + Az^2
  vDSP_vsma(Az, 1, &Bz, distances, 1, distances, 1, N);  // Ax^2 - 2*Ax*Bx + Ay^2 - 2*Ay*By + Az^2 - 2*Az*Bz
  vDSP_vsadd(distances, 1, &constants, distances, 1, N); // ... + constants = (Ax-Bx)^2 + (Ay-By)^2 + (Az-Bz)^2 + scalar
   */

  // take sqrt
  vvsqrtf(distances, distances, &N);
}