3D 中的矢量加速欧几里得距离

Question

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

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

矩阵运算如下：

我正在计算点向量 (A) 和参考点 (B) 之间的欧几里得距离。
这些点位于 3D 空间中，每个点都有一组 X、Y 和 Z 坐标。
因此，点的向量由三个浮点数组来描述，其中包含每个点的 X、Y、Z 坐标。
输出是另一个长度为 N 的向量，其中包含数组中每个点与参考点之间的距离。

三个 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 库中，唯一可用于计算向量之间距离的函数是：

vDSP_vdist()
vDSP_distancesq()
vDSP_vpythg()

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

有几点需要注意：
(1) 我不是在比较距离，所以我不能忍受平方距离。我需要实际距离，因此绝对有必要计算平方根。
(2) 如果您真的认为这样做会显着加快代码速度，那么取平方根倒数是可能的。

我的印象是我没有充分利用 Accelerate 框架。
我正在寻找更智能、更简洁的东西，在更少的函数调用中做更多的工作。以其他方式重新排列内存也可以，但是我认为内存布局还是不错的。

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

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

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

谢谢。

Answer 1

试试这个。

• 在我的 iMac 和 iPhone 上，对于大型 N = 2^20 @ 1000 次重复，它的性能提高了大约 20%
• 此外，matrix 可以被视为严格只读，因为只有distances 被操纵
• 它在代数上等价，但在数值上与您的实现不等价；预计输出差异在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);
}