无分支K-means（或其他优化）

6 回答

• 3

  太糟糕了，我们可以对此进行测试，只是编译它以查看是否存在语法错误并查看程序集是否有意义（它是大多数情况下，虽然GCC溢出 min_index 甚至还有一些未使用的 xmm 寄存器，但不确定为什么会发生这种情况）

  int find_closest(float *x, float *y, float *z,
                   float pt_x, float pt_y, float pt_z, int n) {
      __m128i min_index = _mm_set_epi32(3, 2, 1, 0);
      __m128 xdif = _mm_sub_ps(_mm_set1_ps(pt_x), _mm_load_ps(x));
      __m128 ydif = _mm_sub_ps(_mm_set1_ps(pt_y), _mm_load_ps(y));
      __m128 zdif = _mm_sub_ps(_mm_set1_ps(pt_z), _mm_load_ps(z));
      __m128 min_dist = _mm_add_ps(_mm_add_ps(_mm_mul_ps(xdif, xdif), 
                                              _mm_mul_ps(ydif, ydif)), 
                                              _mm_mul_ps(zdif, zdif));
      __m128i index = min_index;
      for (int i = 4; i < n; i += 4) {
          xdif = _mm_sub_ps(_mm_set1_ps(pt_x), _mm_load_ps(x + i));
          ydif = _mm_sub_ps(_mm_set1_ps(pt_y), _mm_load_ps(y + i));
          zdif = _mm_sub_ps(_mm_set1_ps(pt_z), _mm_load_ps(z + i));
          __m128 dist = _mm_add_ps(_mm_add_ps(_mm_mul_ps(xdif, xdif), 
                                              _mm_mul_ps(ydif, ydif)), 
                                              _mm_mul_ps(zdif, zdif));
          index = _mm_add_epi32(index, _mm_set1_epi32(4));
          __m128i mask = _mm_castps_si128(_mm_cmplt_ps(dist, min_dist));
          min_dist = _mm_min_ps(min_dist, dist);
          min_index = _mm_or_si128(_mm_and_si128(index, mask), 
                                   _mm_andnot_si128(mask, min_index));
      }
      float mdist[4];
      _mm_store_ps(mdist, min_dist);
      uint32_t mindex[4];
      _mm_store_si128((__m128i*)mindex, min_index);
      float closest = mdist[0];
      int closest_i = mindex[0];
      for (int i = 1; i < 4; i++) {
          if (mdist[i] < closest) {
              closest = mdist[i];
              closest_i = mindex[i];
          }
      }
      return closest_i;
  }
  

  像往常一样，它希望指针是16对齐的 . 此外，填充应该是无限远的点（因此它们永远不会最接近目标） .

  SSE 4.1会让你替换它

  min_index = _mm_or_si128(_mm_and_si128(index, mask), 
                           _mm_andnot_si128(mask, min_index));
  

  这样

  min_index = _mm_blendv_epi8(min_index, index, mask);
  

  这是一个asm版本，为vsyasm制作，测试了一下（似乎工作）

  bits 64
  
  section .data
  
  align 16
  centroid_four:
      dd 4, 4, 4, 4
  centroid_index:
      dd 0, 1, 2, 3
  
  section .text
  
  global find_closest
  
  proc_frame find_closest
      ;
      ;   arguments:
      ;       ecx: number of points (multiple of 4 and at least 4)
      ;       rdx -> array of 3 pointers to floats (x, y, z) (the points)
      ;       r8 -> array of 3 floats (the reference point)
      ;
      alloc_stack 0x58
      save_xmm128 xmm6, 0
      save_xmm128 xmm7, 16
      save_xmm128 xmm8, 32
      save_xmm128 xmm9, 48
  [endprolog]
      movss xmm0, [r8]
      shufps xmm0, xmm0, 0
      movss xmm1, [r8 + 4]
      shufps xmm1, xmm1, 0
      movss xmm2, [r8 + 8]
      shufps xmm2, xmm2, 0
      ; pointers to x, y, z in r8, r9, r10
      mov r8, [rdx]
      mov r9, [rdx + 8]
      mov r10, [rdx + 16]
      ; reference point is in xmm0, xmm1, xmm2 (x, y, z)
      movdqa xmm3, [rel centroid_index]   ; min_index
      movdqa xmm4, xmm3                   ; current index
      movdqa xmm9, [rel centroid_four]     ; index increment
      paddd xmm4, xmm9
      ; calculate initial min_dist, xmm5
      movaps xmm5, [r8]
      subps xmm5, xmm0
      movaps xmm7, [r9]
      subps xmm7, xmm1
      movaps xmm8, [r10]
      subps xmm8, xmm2
      mulps xmm5, xmm5
      mulps xmm7, xmm7
      mulps xmm8, xmm8
      addps xmm5, xmm7
      addps xmm5, xmm8
      add r8, 16
      add r9, 16
      add r10, 16
      sub ecx, 4
      jna _tail
  _loop:
      movaps xmm6, [r8]
      subps xmm6, xmm0
      movaps xmm7, [r9]
      subps xmm7, xmm1
      movaps xmm8, [r10]
      subps xmm8, xmm2
      mulps xmm6, xmm6
      mulps xmm7, xmm7
      mulps xmm8, xmm8
      addps xmm6, xmm7
      addps xmm6, xmm8
      add r8, 16
      add r9, 16
      add r10, 16
      movaps xmm7, xmm6
      cmpps xmm6, xmm5, 1
      minps xmm5, xmm7
      movdqa xmm7, xmm6
      pand xmm6, xmm4
      pandn xmm7, xmm3
      por xmm6, xmm7
      movdqa xmm3, xmm6
      paddd xmm4, xmm9
      sub ecx, 4
      ja _loop
  _tail:
      ; calculate horizontal minumum
      pshufd xmm0, xmm5, 0xB1
      minps xmm0, xmm5
      pshufd xmm1, xmm0, 0x4E
      minps xmm0, xmm1
      ; find index of the minimum
      cmpps xmm0, xmm5, 0
      movmskps eax, xmm0
      bsf eax, eax
      ; index into xmm3, sort of
      movaps [rsp + 64], xmm3
      mov eax, [rsp + 64 + rax * 4]
      movaps xmm9, [rsp + 48]
      movaps xmm8, [rsp + 32]
      movaps xmm7, [rsp + 16]
      movaps xmm6, [rsp]
      add rsp, 0x58
      ret
  endproc_frame
  

  在C：

  extern "C" int find_closest(int n, float** points, float* reference_point);
  

  回复于 2024-05-18T23:42:23+08:00
• 6

  您可以使用无分支三元运算符，有时也称为bitselect（condition？true：false） .
  只需将它用于2名成员，默认为什么都不做 .
  不要担心额外的操作，它们与if语句分支相比没什么 .

  bitselect实现：

  inline static int bitselect(int condition, int truereturnvalue, int falsereturnvalue)
  {
      return (truereturnvalue & -condition) | (falsereturnvalue & ~(-condition)); //a when TRUE and b when FALSE
  }
  
  inline static float bitselect(int condition, float truereturnvalue, float falsereturnvalue)
  {
      //Reinterpret floats. Would work because it's just a bit select, no matter the actual value
      int& at = reinterpret_cast<int&>(truereturnvalue);
      int& af = reinterpret_cast<int&>(falsereturnvalue);
      int res = (at & -condition) | (af & ~(-condition)); //a when TRUE and b when FALSE
      return  reinterpret_cast<float&>(res);
  }
  

  你的循环应该是这样的：

  for (int i=0; i < num_centroids; ++i)
  {
    const ClusterCentroid& cent = centroids[i];
    const float dist = ...;
    bool isSmaeller = dist < pt.min_dist;
  
    //use same value if not smaller
    pt.min_index = bitselect(isSmaeller, i, pt.min_index);
    pt.min_dist = bitselect(isSmaeller, dist, pt.min_dist);
  }
  

  回复于 2024-05-18T23:42:23+08:00
• 10

  C是一种高级语言 . 您假设C源代码中的控制流转换为分支指令是有缺陷的 . 我没有从您的示例中定义某些类型，因此我创建了一个具有类似条件赋值的简单测试程序：

  int g(int, int);
  
  int f(const int *arr)
  {
      int min = 10000, minIndex = -1;
      for ( int i = 0; i < 1000; ++i )
      {
          if ( arr[i] < min )
          {
              min = arr[i];
              minIndex = i;
          }
      }
      return g(min, minIndex);
  }
  

  请注意，使用undefined "g"仅仅是为了防止优化器删除所有内容 . 我将带有-O3和-S的G 4.9.2翻译成x86_64程序集（甚至不需要更改-march的默认值），并且（不过于令人惊讶）结果是循环体不包含分支

  movl    (%rdi,%rax,4), %ecx
  movl    %edx, %r8d
  cmpl    %edx, %ecx
  cmovle  %ecx, %r8d
  cmovl   %eax, %esi
  addq    $1, %rax
  

  除此之外，无分支必然更快的假设也可能存在缺陷，因为新距离的概率正在减少你所观察到的元素越多 . 这不是硬币投掷 . 当编译器生成"as-if"组件时，它们比现在更不具有攻击性，因此发明了"bitselect"技巧 . 我更愿意在尝试重新编写代码之前先看看编译器实际生成的汇编类型，以便编译器能够更好地优化它，或者将结果作为手写汇编的基础 . 如果您想查看SIMD，我建议尝试使用减少数据依赖性的"minimum of minimums"方法（在我的示例中，"min"上的依赖可能是瓶颈） .

  回复于 2024-05-18T23:42:23+08:00
• 21

  首先，我建议在尝试任何代码更改之前，先查看优化构建中的反汇编 . 理想情况下，您希望在汇编级别查看分析器数据 . 这可以显示各种事物，例如：

  • 编译器可能没有生成实际的分支指令 .

  • 具有瓶颈的代码行可能会有比您想象的更多的与之关联的指令 - 例如，dist计算 .

  除此之外还有标准技巧，当你谈论距离计算时，它们通常需要一个平方根 . 你应该在最小平方值的过程结束时做那个平方根 .

  SSE可以使用_mm_min_ps一次处理四个值，而不使用任何分支 . 如果你真的需要速度，那么你想要使用SSE（或AVX）内在函数 . 这是一个基本的例子：

  float MinimumDistance(const float *values, int count)
    {
      __m128 min = _mm_set_ps(FLT_MAX, FLT_MAX, FLT_MAX, FLT_MAX);
      int i=0;
      for (; i < count - 3; i+=4)
      {
          __m128 distances = _mm_loadu_ps(&values[i]);
          min = _mm_min_ps(min, distances);
      }
      // Combine the four separate minimums to a single value
      min = _mm_min_ps(min, _mm_shuffle_ps(min, min, _MM_SHUFFLE(2, 3, 0, 1)));
      min = _mm_min_ps(min, _mm_shuffle_ps(min, min, _MM_SHUFFLE(1, 0, 3, 2)));
  
      // Deal with the last 0-3 elements the slow way
      float result = FLT_MAX;
      if (count > 3) _mm_store_ss(&result, min);
      for (; i < count; i++)
      {
          result = min(values[i], result);
      }
  
      return result;
    }
  

  为获得最佳SSE性能，您应确保在对齐的地址处发生负载 . 如有必要，您可以像上面代码中的最后几个一样处理前几个未对齐的元素 .

  另外需要注意的是内存带宽 . 如果在该循环期间没有使用ClusterCentroid结构的几个成员，那么您将从内存中读取比实际需要更多的数据，因为在高速缓存行大小的块中读取内存，每个块为64字节 .

  回复于 2024-05-18T23:42:23+08:00
• 21

  这可能是双向的，但我会尝试以下结构：

  std::vector<float> centDists(num_centroids); //<-- one for each thread. 
  for (size_t p=0; p<num_points; ++p) {
      Point& pt = points[p];
      for (size_t c=0; c<num_centroids; ++c) {
          const float dist = ...;
          centDists[c]=dist;
      }
      pt.min_idx it= min_element(centDists.begin(),centDists.end())-centDists.begin();    
  }
  

  显然，你现在必须在内存上迭代两次，这可能会损害缓存命中率（你也可以将它分成子范围）但另一方面，每个内部循环都应该易于矢量化和展开 - 所以你只需要衡量它是否值得 .

  即使你坚持你的版本，我也会尝试使用 local variables 来跟踪最小索引和距离，并将结果应用到最后 .
  理性的是，对 pt.min_dist 的每次读取或写入都是通过指针实现的，这取决于编译器的优化，可能会也可能不会降低您的性能 .

  对矢量化很重要的另一件事是将 array of Structs （在这种情况下为蜈蚣）变成 struct of arrays （例如，对于点的每个坐标，例如一个数组），因为这样您就不需要额外的收集指令来加载数据以供SIMD指令使用 . 有关该主题的更多信息，请参见Eric Brumer's talk .

  编辑：我的系统的一些数字（haswell，clang 3.5）：
  我用你的基准测试和我的系统做了一个简短的测试，上面的代码使算法减慢了大约10％ - 基本上，没有任何东西可以被矢量化 .

  然而，当将AoS应用于质心的SoA变换时，距离计算被矢量化，与使用AoS到SoA转换的原始结构相比，导致整体运行时间减少约40％ .

  回复于 2024-05-18T23:42:23+08:00
• 6

  一种可能的微优化：将min_dist和min_index存储在局部变量中 . 编译器可能必须以编写方式更频繁地写入内存;在某些体系结构上，这会对性能产生很大影响 . 有关其他示例，请参阅my answer here .

  亚当斯建议同时做4比较也是一个很好的建议 .

  但是，您最好的加速来自减少必须检查的质心数量 . 理想情况下，在质心周围构建一个kd树（或类似物），然后查询它以找到最近的点 .

  如果您没有任何树构建代码，这里是我最喜欢的“穷人”最近点搜索：

  Sort the points by one coordinate, e.g. cent.pos[0]
  Pick a starting index for the query point (pt)
  Iterate forwards through the candidate points until you reach the end, OR when abs(pt.pos[0] - cent.pos[0]) > min_dist
  Repeat the previous step going the opposite direction.
  

  搜索的额外停止条件意味着您应该跳过相当多的积分;你也保证不会跳过任何比你已经找到的最好的点 .

  所以对于你的代码，这看起来像

  // sort centroid by x coordinate.
  min_index = -1;
  min_dist = numeric_limits<float>::max();
  
  // pick the start index. This works well if the points are evenly distributed.
  float min_x = centroids[0].pos[0];
  float max_x = centroids[num_centroids-1].pos[0];
  float cur_x = pt.pos[0];
  float t = (max_x - cur_x) / (max_x - min_x);
  // TODO clamp t between 0 and 1
  int start_index = int(t * float(num_centroids))
  
  // Forward search
  for (int i=start_index ; i < num_centroids; ++i)
  {
      const ClusterCentroid& cent = centroids[i];
      if (fabs(cent.pos[0] - pt.pos[0]) > min_i)
          // Everything to the right of this must be further min_dist, so break.
          // This is where the savings comes from!
          break; 
      const float dist = ...;
      if (dist < min_dist)
      {
          min_dist = dist;
          min_index = i;
      }
  }
  
  // Backwards search
  for (int i=start_index ; i >= 0; --i)
  {
      // same as above
  }
  pt.min_dist = min_dist
  pt.min_index = min_index
  

  （请注意，这假设您计算点之间的距离，但是您的程序集表明它是距离的平方 . 相应地调整中断条件） .

  构建树或对质心进行排序会产生轻微的开销，但这应该通过在更大的循环中（在点数上）更快地进行计算来抵消 .

  回复于 2024-05-18T23:42:23+08:00