3D 中的向量加速欧氏距离

Question

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

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

矩阵运算如下：

我正在计算点向量 (A) 和参考点 (B) 之间的欧式距离。
这些点是 3D space 并且每个点都有一组 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()

也许我遗漏了什么，但据我所知 none 它们支持计算 3D 距离所需的三个输入向量。

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

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

我也乐于接受有关在 Intel 处理器上工作的其他高度 optimized/vectorized 线性代数库的建议。我不在乎它们是商业解决方案还是开源解决方案，只要它们的性能快速且稳健即可。

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

我正在 Xcode 7 中开发 MacBook Pro（视网膜，15 英寸，2014 年中）运行 Mac OS X埃尔卡皮坦。

谢谢。

Answer 1

试试这个。

它在我的 iMac 和 iPhone

N = 2^20

此外，matrix可以严格地只读，因为只有distances被操作
它与您的实现在代数上等价，但在数值上不等价；期望输出差异在 10^-6

在我看来，vDSP 级别太高，无法进一步 "targeted" 优化。相反，您可以查看 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);
}

3D 中的向量加速欧氏距离

Vector Accelerated Euclidean Distance in 3D

c

matrix

3d

euclidean-distance

vdsp