用于非常特殊情况的快速点积答案

【问题标题】：Fast dot product for a very special case用于非常特殊情况的快速点积
【发布时间】：2010-05-13 10:00:02
【问题描述】：

给定一个大小为 L 的向量 X，其中 X 的每个标量元素都来自二进制集合 {0,1}，如果大小为 L 的向量 Y 由的整数值元素。我建议，必须有一种非常快速的方法来做到这一点。

假设我们有L=4; X[L]={1, 0, 0, 1}; Y[L]={-4, 2, 1, 0}，我们必须找到z=X[0]*Y[0] + X[1]*Y[1] + X[2]*Y[2] + X[3]*Y[3]（在这种情况下会给我们-4）。

很明显，X 可以用二进制数字表示，例如L=32 的整数类型 int32。然后，我们要做的就是找到这个整数与 32 个整数数组的点积。您对如何快速完成有任何想法或建议？

【问题讨论】：

您的数据有多大和有多稀疏？
“很明显，X 可以用二进制数字表示”——是的，但这并不明显会带来任何性能改进。当然，除非内存大小很重要。是这样吗？
您的特殊情况是否足够特殊以至于您可以依赖扩展？例如，对于任何类型的 SIMD，使用向量运算可能比将 X 打包到位域中获得更多的速度优势。您必须进行比较。
@Konrad：对于一般算法，内存大小无关紧要
L 是固定常数吗？如果是这样，那么它可能会对最佳解决方案产生一些影响......您可以提供的任何信息/约束通常都有助于找到更快的解决方案。

标签： c++ c algorithm math

【解决方案1】：

这确实需要分析，但您可能需要考虑另一种方法：

int result=0;
int mask=1;
for ( int i = 0; i < L; i++ ){
    if ( X & mask ){
        result+=Y[i];
    }
    mask <<= 1;
}

通常位移和按位运算比乘法快，但是，if 语句可能比乘法慢，尽管使用分支预测和大 L 我猜它可能会更快。但是，您确实必须对其进行分析，以确定它是否会导致任何加速。

正如在下面的 cmets 中所指出的，手动或通过编译器标志（例如 GCC 上的“-funroll-loops”）展开循环也可以加快速度（省略循环条件）。

编辑
在下面的 cmets 中，提出了以下很好的调整：

int result=0;
for ( int i = 0; i < L; i++ ){
    if ( X & 1 ){
        result+=Y[i];
    }
    X >>= 1;
}

【讨论】：

感谢有趣的解决方案，但 X 不能是数组 - 它应该是整数。
@psihodelia，不明白这个问题。我已经更新了。
也可以用-bit替换掩码[bit]
好的，现在你的解决方案好多了！然而，有趣的是，它是否可以进一步优化，例如if(X&(1
@psihodelia，是的，可以稍微优化一下，让它只移动一个。我刚刚进行了更改。

【解决方案2】：

研究 SSE2 的建议是否有用？它已经具有点积类型的操作，而且您可以轻松地并行执行 4 次（或者可能是 8 次，我忘记了寄存器大小）简单循环的简单迭代。 SSE 也有一些简单的逻辑类型的操作，因此它可以在不使用任何条件操作的情况下进行加法而不是乘法......再次，您必须查看可用的操作。

【讨论】：

不，谢谢。 SSE 是特定于机器的实现。我正在寻找一种通用算法。
@psihodelia：你绝对应该考虑 SSE；这个问题将从流处理中受益匪浅。
@psihodelia：你没有任何意义。一般算法只能在渐近性能方面进行有意义的比较。如果您想比较性能，您需要在实际设置（代码和数据缓存压力等）中进行特定于机器的基准测试。@John：SSE4.1 会更有帮助，它有扩展指令来扩展位掩码。我整理了一个快速示例，它在 48 条指令中执行 L=32，没有分支，一个额外数据的缓存行：gist.github.com/400685
这很酷。尽管即使我认为 SSE4 还不足以成为目前的主流目标。

【解决方案3】：

试试这个：

int result=0;
for ( int i = 0; i < L; i++ ){
    result+=Y[i] & (~(((X>>i)&1)-1));
}

这避免了条件语句，并使用按位运算符用零或一屏蔽标量值。

【讨论】：

也适用的按位运算符的替代方案是：result+=Y[i] & -((X>>i)&1);不确定哪个更快....
从逻辑上讲，您的第二种方法包含较少的操作，因此“应该”更快
是的，我同意......尽管多年来我已经学会了在编译器优化和微处理器怪癖方面不要假设任何事情！
实际上它们的运行速度与我测试过的相同。
酷，感谢 KennyTM（我手边没有 C 编译器！）。我一直在研究标准假设，即按位运算在现代处理器上基本上是免费的：-)

【解决方案4】：

由于大小明确并不重要，我认为以下可能是最有效的通用代码：

int result = 0;
for (size_t i = 0; i < 32; ++i)
    result += Y[i] & -X[i];

位编码X 不会带来任何好处（即使循环可能会提前终止，正如@Mathieu 正确指出的那样）。但是在循环中省略if确实。

当然，正如其他人所指出的，循环展开可以大大加快这一速度。

【讨论】：

我认为，在某些架构上，Y[i]*X[i] 比 AND 和否定需要更少的指令。
@psihodelia：更少的指令并不一定意味着更少的运行时间（尤其是在 CISC 架构上）。
@psihodelia：这无关紧要。你可以很确定否定 + 按位，并且将总是比乘法快。
对于一般算法而言，涉及多少操作并不重要。在 RISC 架构上，一个操作应该比两个快。

【解决方案5】：

此解决方案与 Micheal Aaron 的解决方案相同，但速度稍快（根据我的测试）：

long Lev=1;
long Result=0
for (int i=0;i<L;i++) {
  if (X & Lev)
     Result+=Y[i];
  Lev*=2;
}

我认为有一种数字方法可以快速建立单词中的下一个设置位，如果您的 X 数据非常稀疏但目前无法找到所述数字公式，则应该可以提高性能。

【讨论】：

MA 已修改他的解决方案以包含此优化。他有我的赞成票

【解决方案6】：

我已经看到了一些带有技巧的响应（以避免分支），但没有一个得到正确的循环恕我直言：/

优化@Goz答案：

int result=0;
for (int i = 0, x = X; x > 0; ++i, x>>= 1 )
{
   result += Y[i] & -(int)(x & 1);
}

优点：

不需要每次都做i的位移操作(X>>i)
如果X 的高位包含 0，则循环会更快停止

现在，我确实想知道它是否运行得更快，特别是因为 for 循环的过早停止可能不那么容易展开循环（与编译时常量相比）。

【讨论】：

不错。如何在循环开始时增加 i 以跳过 X 中的零？这是额外的代码，但它避免了内存访问，所以它应该是一个净赢.....
一直在考虑这个问题：因为我们已经努力通过按位操作来消除分支，所以用可变长度循环将它们重新添加回来可能没有意义，至少在以下情况下L 是一个小的固定数。如果 1 和 0 是随机的 50%，那么您平均只在向量的每一侧节省一次内存访问。
+1 很早就出局了。但是，如果最高位有 1，则它不会提供胜利，因为您仍然可以进行同样多的右移。但是，编译器的流水线可能更容易。
你能向int解释一下演员阵容吗？这不是很没必要吗？（x & 1 是 int。）
@Rudolph：我不知道，我只是从Goz 复制了这一点。 @mikera/Goz：我也不确定终止条件，因为我担心它可能会随着循环展开而搞砸，但是在每个循环中移位 1 确实看起来更快:)

【解决方案7】：

如何将移位循环与小型查找表结合起来？

    int result=0;

    for ( int x=X; x!=0; x>>=4 ){
        switch (x&15) {
            case 0: break;
            case 1: result+=Y[0]; break;
            case 2: result+=Y[1]; break;
            case 3: result+=Y[0]+Y[1]; break;
            case 4: result+=Y[2]; break;
            case 5: result+=Y[0]+Y[2]; break;
            case 6: result+=Y[1]+Y[2]; break;
            case 7: result+=Y[0]+Y[1]+Y[2]; break;
            case 8: result+=Y[3]; break;
            case 9: result+=Y[0]+Y[3]; break;
            case 10: result+=Y[1]+Y[3]; break;
            case 11: result+=Y[0]+Y[1]+Y[3]; break;
            case 12: result+=Y[2]+Y[3]; break;
            case 13: result+=Y[0]+Y[2]+Y[3]; break;
            case 14: result+=Y[1]+Y[2]+Y[3]; break;
            case 15: result+=Y[0]+Y[1]+Y[2]+Y[3]; break;
        }
        Y+=4;
    }

这将取决于编译器在优化 switch 语句方面的表现，但根据我的经验，他们现在在这方面做得很好......

【讨论】：

我比我更喜欢你的解决方案。

【解决方案8】：

这个问题可能没有一般的答案。您需要在所有不同的目标下分析您的代码。性能将取决于编译器优化，例如大多数现代 CPU（x86、PPC、ARM 都有自己的实现）上可用的循环展开和 SIMD 指令。

【讨论】：

【解决方案9】：

对于 small L，您可以使用 switch 语句代替循环。例如，如果 L = 8，您可以：

int dot8(unsigned int X, const int Y[])
{
    switch (X)
    {
       case 0: return 0;
       case 1: return Y[0];
       case 2: return Y[1];
       case 3: return Y[0]+Y[1];
       // ...
       case 255: return Y[0]+Y[1]+Y[2]+Y[3]+Y[4]+Y[5]+Y[6]+Y[7];
    }
    assert(0 && "X too big");
}

如果 L = 32，您可以编写一个 dot32() 函数，该函数调用 dot8() 四次次，尽可能内联。（如果您的编译器拒绝内联 dot8()，您可以将 dot8() 重写为宏以强制内联。）添加：

int dot32(unsigned int X, const int Y[])
{
    return dot8(X >> 0  & 255, Y + 0)  +
           dot8(X >> 8  & 255, Y + 8)  +
           dot8(X >> 16 & 255, Y + 16) +
           dot8(X >> 24 & 255, Y + 24);
}

正如 mikera 指出的那样，这种解决方案可能会产生指令缓存成本；如果是这样，使用 dot4() 函数可能会有所帮助。

进一步更新：这可以结合 mikera 的解决方案：

static int dot4(unsigned int X, const int Y[])
{
    switch (X)
    {
        case 0: return 0;
        case 1: return Y[0];
        case 2: return Y[1];
        case 3: return Y[0]+Y[1];
        //...
        case 15: return Y[0]+Y[1]+Y[2]+Y[3];
    }
}

在 CYGWIN 上使用 gcc 4.3.4 使用 -S -O3 选项查看生成的汇编代码，我有点惊讶地看到它在 dot32() 中自动内联，带有 8 个 16 项跳转表。

但添加 __attribute__((__noinline__)) 似乎会产生更好看的汇编程序。

另一种变体是在 switch 语句中使用 fall-throughs，但 gcc 添加了 jmp 指令，看起来并没有更快。

编辑--全新的答案：在考虑了 Ants Aasma 提到的 100 个循环惩罚和其他答案之后，以上可能不是最优的。相反，您可以手动展开循环，如下所示：

int dot(unsigned int X, const int Y[])
{
    return (Y[0] & -!!(X & 1<<0)) +
           (Y[1] & -!!(X & 1<<1)) +
           (Y[2] & -!!(X & 1<<2)) +
           (Y[3] & -!!(X & 1<<3)) +
           //...
           (Y[31] & -!!(X & 1<<31));
}

这在我的机器上生成 32 x 5 = 160 条快速指令。可以想象，智能编译器可以展开其他建议的答案以给出相同的结果。

但我仍在仔细检查。

【讨论】：

我喜欢这个！我唯一的问题是指令缓存成本是否超过了操作数量方面的明显优势？
@mikera：计时。这是唯一的判断方法。唉，大规模展开的版本很可能不会更快，因为处理器更难将所有代码保存在其 ICache 中（并且还做其他有用的工作）。
更大的问题是这个有8个间接跳转。那是 100 多个循环的惩罚，不包括实际添加东西的时间。

【解决方案10】：

result = 0;
for(int i = 0; i < L ; i++)
    if(X[i]!=0)
      result += Y[i];

【讨论】：

这比if (X[i] == 1) 有什么优势？为什么这不是一个好主意？
其实不如if(X[i] == 1)。这不是一个好主意，因为在条件语句中包含除bool 之外的任何内容都被认为是不好的。 :)
@TheMachineCharmer：那为什么不在你的答案中写上X[i] == 1？从效率的角度来看，它是相同的（！），这将是一个很好的答案，即使它使用了一个已被其他答案完全避免的条件。
通常与零的比较比与其他值的比较快，因为大多数机器都有一个零时跳转、非零跳转和其他基于零的跳转语句。与其他值比较通常意味着减法，然后与零进行比较。
@Michael Aaron Safyan：这是个好主意，回答了 Konrad Rudalf 的第一个问题并澄清了我的疑问。谢谢！ :)

【解决方案11】：

从主内存加载X 和Y 所花费的时间很可能占主导地位。如果您的 CPU 架构就是这种情况，那么加载较少时算法会更快。这意味着将X 存储为位掩码并将其扩展至一级缓存将加速整个算法。

另一个相关问题是您的编译器是否会为Y 生成最佳负载。这高度依赖于 CPU 和编译器。但总的来说，如果编译器能够准确地看到什么时候需要哪些值，它会有所帮助。您可以手动展开循环。但是，如果 L 是常量，则将其留给编译器：

template<int I> inline void calcZ(int (&X)[L], int(&Y)[L], int &Z) {
  Z += X[I] * Y[I]; // Essentially free, as it operates in parallel with loads.
  calcZ<I-1>(X,Y,Z);
}
template< > inline void calcZ<0>(int (&X)[L], int(&Y)[L], int &Z) {
  Z += X[0] * Y[0];
}
inline int calcZ(int (&X)[L], int(&Y)[L]) {
    int Z = 0;
    calcZ<L-1>(X,Y,Z);
    return Z;
}

（Konrad Rudolph 在评论中对此提出质疑，想知道内存使用。这不是现代计算机架构中真正的瓶颈，内存和 CPU 之间的带宽才是。如果 Y 以某种方式存在，这个答案几乎无关紧要已经在缓存中了。）

【讨论】：

注：我不怀疑更少的内存传输可能会提高性能。我质疑这是否“显而易见”。当然，缓存是一个非常重要的因素。

【解决方案12】：

您可以将位向量存储为 int 序列，其中每个 int 将几个系数打包为位。然后，逐分量乘法相当于位与。有了这个，您只需要计算可以这样完成的设置位数：

inline int count(uint32_t x) {
    // see link
}

int dot(uint32_t a, uint32_t b) {
    return count(a & b);
}

有关计算设置位的一些技巧，请参阅http://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel

编辑：抱歉，我刚刚意识到只有一个向量包含 {0,1} 的元素，而另一个不包含。此答案仅适用于两个向量都限制为 {0,1} 集合中的系数的情况。

【讨论】：

【解决方案13】：

使用x[i] = 1 所在位置的链表表示X。要找到所需的总和，您需要 O(N) 操作，其中 N 是您列表的大小。

【讨论】：

如果 X 是稀疏的并且 X 是相对静态的，这是一个好主意。显然，最快的实现将是某种连续数组，而不是此处建议的链表。
@Elemental 你的意思是稀疏的、相对静态的和大的。
是的，提问者暗示 X 的大小小于 32，这绝对不算大。

【解决方案14】：

好吧，如果它是 1，你希望所有位都过去，如果它是 0，则没有。所以你想以某种方式将 1 变成 -1（即 0xffffffff）并且 0 保持不变。那只是-X ....所以你做...

Y & (-X)

对于每个元素...工作完成了吗？

Edit2：举一个代码示例，您可以这样做并避免分支：

int result=0;
for ( int i = 0; i < L; i++ )
{
   result+=Y[i] & -(int)((X >> i) & 1);
}

当然，最好将 1 和 0 保留在整数数组中，从而避免移位。

编辑：还值得注意的是，如果 Y 中的值大小为 16 位，那么您可以执行其中的 2 个操作，并且每个操作都可以进行操作（如果您有 64 位寄存器，则可以执行 4 个操作）。不过，这确实意味着将 X 值 1 乘 1 取反为更大的整数。

ie YVals = -4, 3 in 16-bit = 0xFFFC, 0x3 ... 放入 1 32-bit，你得到 0xFFFC0003。如果您有 1, 0 作为 X val，那么您将形成 0xFFFF0000 和 2 的位掩码，并且您在 1 位和运算中得到 2 个结果。

另一个编辑：

如果您想要关于如何执行第二种方法的代码类似，这应该可以工作（尽管它利用了未指定的行为，因此它可能不适用于每个编译器.. 适用于我的每个编译器'虽然遇到过）。

union int1632
{
     int32_t i32;
     int16_t i16[2];
};

int result=0;
for ( int i = 0; i < (L & ~0x1); i += 2 )
{
    int3264 y3264;
    y3264.i16[0] = Y[i + 0];
    y3264.i16[1] = Y[i + 1];

    int3264 x3264;
    x3264.i16[0] = -(int16_t)((X >> (i + 0)) & 1);
    x3264.i16[1] = -(int16_t)((X >> (i + 1)) & 1);

    int3264 res3264;
    res3264.i32  = y3264.i32 & x3264.i32;

    result += res3264.i16[0] + res3264.i16[1];    
}

if ( i < L )
    result+=Y[i] & -(int)((X >> i) & 1);

希望编译器会优化分配（我不确定，但这个想法可以重新工作，这样它们肯定是）并给你一个小的加速，因为你现在只需要做 1 按位而不是 2。虽然速度提升会很小......

【讨论】：

我已经检查了您的解决方案，但它给了我完全错误的答案。必须找到二进制字与整数数组的点积。
嗯，你必须早点把结果加起来。您还需要将 Y 向量中的值提升为整数才能正常工作。
代码示例有帮助吗？您仍然必须对两种方式进行分析，因为分支更便宜是合理的……但是，我的方法是恒定时间。
@Goz - 查看了文本中的工作示例；你的算法无法做到这一点。
@Goz - 是的，谢谢！你的算法确实很好。我已经检查过了，效果很好。