logo标签列表

非常特殊情况下的快速点积

c++calgorithmmath
15
给定一个大小为L的向量X，其中X的每个标量元素都来自二进制集合{0,1}，要找到一个点积z=dot(X,Y)，如果向量Y的大小为L，则其由整数值元素组成。我建议有一种非常快速的方法可以做到这一点。
假设我们有，并且我们必须找到z=X[0]*Y[0] + X[1]*Y[1] + X[2]*Y[2] + X[3]*Y[3](在这种情况下将给我们-4)。
很显然，X可以用二进制数字表示，例如，对于L=32，使用int32类型的整数。然后，我们所要做的就是找到该整数与包含32个整数的数组的点积。你有什么想法或建议如何非常快地做到这一点吗？
1
你的数据有多大，有多稀疏？ - Chris H
1
“很明显，X可以用二进制数字表示” - 是的，但这并不明显会带来任何性能提升。除非，当然，内存大小很重要。是这种情况吗？ - Konrad Rudolph
2
你的特殊情况是否足够特殊，可以依赖于扩展？例如，使用任何类型的SIMD，您可能会通过使用向量操作获得比将X打包到位域中获得更多的速度优势。您需要进行比较。 - Steve Jessop
@Konrad：对于一般算法，内存大小并不重要。 - psihodelia
L是一个固定常数吗？如果是的话，它可能对最优解产生一些影响...您能提供的任何信息/约束通常都有助于找到更快的解决方案。 - mikera
14个回答
7

这确实需要进行分析，但你可能需要考虑另一种选择：

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

通常位移和按位操作比乘法更快，然而if语句可能比乘法慢，尽管使用分支预测和大型L，我猜它可能会更快。但是你真的需要对其进行分析，以确定是否会加速。

正如下面的评论所指出的那样，展开循环可以通过手动或通过编译器标志（例如在GCC上的“-funroll-loops”）来加速此过程（省略循环条件）。

编辑
在下面的评论中，提出了以下良好的调整建议：

int result=0;
for ( int i = 0; i < L; i++ ){
    if ( X & 1 ){
        result+=Y[i];
    }
    X >>= 1;
}
如果你的L值为32左右，可以考虑循环展开。这样完全消除循环和条件语句中的条件指令，以最大限度地利用并行性。至少在Intel/AMD CPU上，这段代码最昂贵的部分似乎是分支操作。 - Kris
我认为，在一些 MUL 操作只需要一个 CPU 时钟周期的架构上，result+=Y[i]*X&1; 会更快，因为不需要分支。 - psihodelia
消除条件语句！result += Y[i] & -(X&1); - R.. GitHub STOP HELPING ICE
把 if 给去掉，然后这样写：result+=Y[i] * (X & 1); - The Nomadic Coder
@SemicolonWarrier 无论哪种方式，都需要进行性能分析...没有性能分析的情况下，无法确定分支预测或乘法哪个更优。使用显式的分支预测提示（例如GCC的__builtin_expect），也是值得尝试的。在所有情况下，这都需要进行性能分析。 - Michael Aaron Safyan
显示剩余8条评论
4

试试这个：

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

这样做避免了使用条件语句，而是使用位运算符将标量值与零或一进行屏蔽。
另一种同样有效的位运算符替代方法是：result+=Y[i] & -((X>>i)&1);不确定哪个更快.... - mikera
从逻辑上讲，你的第二种方法包含了一个更少的操作，因此“应该”更快。 - Goz
是的，我同意... 不过多年来我学到了不要轻易假设编译器优化和微处理器特性！ - mikera
4
根据我的测试，它们实际上以相同的速度运行。 - kennytm
非常酷，谢谢KennyTM（我手头没有C编译器！）。我已经很长一段时间以来一直在基于标准假设下开展位运算，在现代处理器上基本上是免费的 :-)。 - mikera
4

研究一下SSE2是否有帮助？它已经具备点积运算，此外，您还可以轻松地并行执行4（或者可能是8，我忘记了寄存器大小）个简单迭代的原始循环。

SSE还具有一些简单的逻辑操作，因此它可能能够进行加法而不使用任何条件操作... 同样，您需要查看可用的操作。

不，谢谢。SSE是一种特定于机器的实现。我正在寻找一种通用算法。 - psihodelia
1
@psihodelia：你一定要考虑使用SSE；这个问题会从流处理中获益匪浅。 - Konrad Rudolph
2
@psihodelia：你说的话毫无意义。通用算法只能从渐近性能方面有意义地进行比较。如果你想比较性能，你需要针对具体的机器，并在现实环境中进行基准测试（代码和数据缓存压力等）。@John：SSE4.1会更有帮助，它具有扩展指令来扩展位掩码。我快速编写了一个示例，在48个指令中完成L=32，没有分支，额外的一行缓存数据：http://gist.github.com/400685 - Ants Aasma
这很酷。尽管我认为SSE4现在还不够主流。 - Mr. Boy
3

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

int result = 0;
for (size_t i = 0; i < 32; ++i)
    result += Y[i] & -X[i];
< p > 二进制编码 < code > X 没有任何实际作用（即使循环可能根据 @Mathieu 的正确说明提前终止）。但是在循环内省略 < code > if 却能有所帮助。< /p> < p > 当然，像其他人指出的那样，循环展开可以极大地加速此过程。< /p>
我认为，在某些体系结构中，Y[i]*X[i]所需的指令比AND和否定少。 - psihodelia
@psihodelia：指令越少并不一定意味着运行时间更短（特别是在CISC架构上）。 - Drew Hall
@psihodelia：那不相关。你可以非常确定否定+按位与运算将始终比乘法更快。 - Konrad Rudolph
对于一般算法来说，所涉及的操作次数确实很重要。在RISC架构中，一次操作应该比两次更快。 - psihodelia
2

这个解决方案与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已经修改了他的解决方案，包括这个优化。他得到了我的支持。 - Elemental
2

我看到很多人用位运算技巧（避免分支）来回答问题，但在我看来，没有一个人得出了正确的循环。

优化@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开头的零，这样可以避免内存访问，虽然会增加一些额外的代码，但总体上应该是有益的。 - mikera
我一直在思考这个问题：既然我们已经通过位运算消除了分支，那么使用可变长度循环来添加分支似乎就不是很明智，至少对于L为小固定数字的情况而言。如果1和0随机分布达到50%，那么平均每个向量的两侧只能节省一个内存访问。 - mikera
+1 早期退出很好。但是，如果最高位为1，则不会提供胜利，因为您仍然需要执行同样多的右移操作。但是，这可能更容易为编译器进行流水线处理。 - Goz
请问能否解释一下将类型转换为 int 的原因？这样做不是有些多余吗？(x & 1 已经是一个 int。) - Konrad Rudolph
@Rudolph：我不知道，我只是从Goz中复制了这一部分。@mikera/Goz：我也不确定终止条件，因为我担心它可能会影响循环展开，但在每个循环中进行位移似乎更快 :) - Matthieu M.
2
如何将移位循环与小型查找表结合使用？
    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语句的能力，但在我的经验中，现代编译器在这方面做得相当不错...
我喜欢你的解决方案比我的更好。 - Joseph Quinsey
1
result = 0;
for(int i = 0; i < L ; i++)
    if(X[i]!=0)
      result += Y[i];
这种方法比 if (X[i] == 1) 更好在哪里？为什么它不是一个好主意？ - Konrad Rudolph
它实际上比if(X[i] == 1)差劣。 这不是一个好主意，因为在条件语句中除了bool以外的任何东西都被认为是不好的。 :) - Pratik Deoghare
@TheMachineCharmer：那么为什么你的答案中不写X[i] == 1呢？从效率的角度来看，这是相同的（！），尽管其他答案完全避免了条件语句，但这将是一个好的答案。 - Konrad Rudolph
2
通常与零进行比较比与其他值进行比较更快，因为大多数机器都有基于零的跳转指令，例如跳转如果为零、跳转如果不是零等。与其他值进行比较通常意味着进行减法运算，然后再比较是否为零。 - Michael Aaron Safyan
@Michael Aaron Safyan：这是一个很棒的主意，回答了Konrad Rudalf的第一个问题并澄清了我的疑惑。谢谢！ :) - Pratik Deoghare
1

对于小型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()内自动内联，并带有八个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个快速指令。一个聪明的编译器可以想象展开其他建议的答案以得到相同的结果。

但我仍在仔细检查。

1
我喜欢这个！我的唯一问题是指令缓存成本是否超过了操作数量的明显优势？ - mikera
1
@mikera: 计时。这是唯一的方法来判断。不幸的是，由于处理器难以在其ICache中保存所有代码并执行其他有用的工作，因此大规模展开版本可能不会更快。 - Donal Fellows
2
更大的问题是这个程序有8个间接跳转。仅仅是惩罚的周期就超过了100个，这还不包括实际添加代码的时间。 - Ants Aasma
1

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

另一个相关的问题是您的编译器是否会为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已经在缓存中，那么这个答案几乎无关紧要。）
请注意，我并不质疑较少的内存传输可能会提高性能；我质疑这是否“显而易见”。当然，缓存是一个非常重要的因素。 - Konrad Rudolph
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接