1

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

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

矩阵运算如下：

我正在计算点向量 (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 中进行开发。

谢谢。

4

1 回答 1

1

试试这个。

  • N = 2^20在我的 iMac 和 iPhone 上，对于大型 @1000 次重复，它的性能提高了大约 20%
  • 此外，matrix可以被视为严格只读，因为 onlydistances被操纵
  • 它在代数上等价，但在数字上与您的实现不等价；期望输出差异在附近10^-6

在我看来，vDSP对于进一步“有针对性”的优化来说，水平太高了。相反，您可以查看Ray Wenderlich 的 iOS 组装教程作为使用 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);
}
于 2016-09-26T06:09:04.647 回答