文档建议反馈控制台
首页
学习
活动
专区
工具
TVP
最新优惠活动
文章/答案/技术大牛
发布
社区首页 >问答首页 >如何在python中加速多个内部产品

问如何在python中加速多个内部产品
EN

Stack Overflow用户
提问于 2014-11-17 02:53:53
回答 5查看 923关注 0票数 17

我有一些执行以下操作的简单代码。

它遍历所有可能的长度为n的列表F和+-1个条目。对于每个条目，它都会迭代所有可能长度的2n列表S和+-1条目，其中$S$的前一半只是后半部分的副本。代码计算F与长度为n的S的每个子列表的内积。对于每个F，S，它计算直到第一个非零的内积为零的内积。

下面是代码。

代码语言：javascript
复制
#!/usr/bin/python

from __future__ import division
import itertools
import operator
import math

n=14
m=n+1
def innerproduct(A, B):
    assert (len(A) == len(B))
    s = 0 
    for k in xrange(0,n):
        s+=A[k]*B[k]
    return s

leadingzerocounts = [0]*m
for S in itertools.product([-1,1], repeat = n):
    S1 = S + S
    for F in itertools.product([-1,1], repeat = n):
        i = 0
        while (i<m):
            ip = innerproduct(F, S1[i:i+n])
            if (ip == 0):
                leadingzerocounts[i] +=1
                i+=1
            else:
                break

print leadingzerocounts

n=14的正确输出为

代码语言：javascript
复制
[56229888, 23557248, 9903104, 4160640, 1758240, 755392, 344800, 172320, 101312, 75776, 65696, 61216, 59200, 59200, 59200]

使用pypy，当n= 14时，这需要1分18秒。不幸的是，我真的想运行16,18,20,22,24,26。我不介意使用numba或cython，但如果可能的话，我更愿意使用python。

任何帮助加快这一进程的人都非常感谢。

我将在这里记录下最快的解决方案。(如果我错过了更新的答案，请告诉我。)

  • n = 22 at 9m35.081s by Eisenstat (C)
  • n = 18 at 1m16.344s by Eisenstat (pypy)
  • n = 18 at 2m54.998s by Tupteq (pypy)
  • n = 14 at 26s by Neil (numpy)
  • n - 14 at 11m59.192s by kslote1
performance
algorithm
cython
numba
python
EN

回答 5

Stack Overflow用户

回答已采纳

发布于 2014-11-20 03:24:32

这个新的代码通过利用问题的循环对称性获得了另一个数量级的加速。此Python版本使用Duval算法枚举项链；C版本使用蛮力。两者都包含了下面描述的加速。在我的机器上，C版本在100秒内解决了n= 20！粗略的计算表明，如果你让它在单个内核上运行一周，它可以解决n= 26，并且，如下所述，它符合并行性。

代码语言：javascript
复制
import itertools


def necklaces_with_multiplicity(n):
    assert isinstance(n, int)
    assert n > 0
    w = [1] * n
    i = 1
    while True:
        if n % i == 0:
            s = sum(w)
            if s > 0:
                yield (tuple(w), i * 2)
            elif s == 0:
                yield (tuple(w), i)
        i = n - 1
        while w[i] == -1:
            if i == 0:
                return
            i -= 1
        w[i] = -1
        i += 1
        for j in range(n - i):
            w[i + j] = w[j]


def leading_zero_counts(n):
    assert isinstance(n, int)
    assert n > 0
    assert n % 2 == 0
    counts = [0] * n
    necklaces = list(necklaces_with_multiplicity(n))
    for combo in itertools.combinations(range(n - 1), n // 2):
        for v, multiplicity in necklaces:
            w = list(v)
            for j in combo:
                w[j] *= -1
            for i in range(n):
                counts[i] += multiplicity * 2
                product = 0
                for j in range(n):
                    product += v[j - (i + 1)] * w[j]
                if product != 0:
                    break
    return counts


if __name__ == '__main__':
    print(leading_zero_counts(12))

C版本：

代码语言：javascript
复制
#include <stdio.h>

enum {
  N = 14
};

struct Necklace {
  unsigned int v;
  int multiplicity;
};

static struct Necklace g_necklace[1 << (N - 1)];
static int g_necklace_count;

static void initialize_necklace(void) {
  g_necklace_count = 0;
  for (unsigned int v = 0; v < (1U << (N - 1)); v++) {
    int multiplicity;
    unsigned int w = v;
    for (multiplicity = 2; multiplicity < 2 * N; multiplicity += 2) {
      w = ((w & 1) << (N - 1)) | (w >> 1);
      unsigned int x = w ^ ((1U << N) - 1);
      if (w < v || x < v) goto nope;
      if (w == v || x == v) break;
    }
    g_necklace[g_necklace_count].v = v;
    g_necklace[g_necklace_count].multiplicity = multiplicity;
    g_necklace_count++;
   nope:
    ;
  }
}

int main(void) {
  initialize_necklace();
  long long leading_zero_count[N + 1];
  for (int i = 0; i < N + 1; i++) leading_zero_count[i] = 0;
  for (unsigned int v_xor_w = 0; v_xor_w < (1U << (N - 1)); v_xor_w++) {
    if (__builtin_popcount(v_xor_w) != N / 2) continue;
    for (int k = 0; k < g_necklace_count; k++) {
      unsigned int v = g_necklace[k].v;
      unsigned int w = v ^ v_xor_w;
      for (int i = 0; i < N + 1; i++) {
        leading_zero_count[i] += g_necklace[k].multiplicity;
        w = ((w & 1) << (N - 1)) | (w >> 1);
        if (__builtin_popcount(v ^ w) != N / 2) break;
      }
    }
  }
  for (int i = 0; i < N + 1; i++) {
    printf(" %lld", 2 * leading_zero_count[i]);
  }
  putchar('\n');
  return 0;
}

您可以通过利用符号对称性(4x)并仅迭代通过第一个内积测试的向量(渐近地，O(sqrt(n))x)来获得一些加速。

代码语言：javascript
复制
import itertools


n = 10
m = n + 1


def innerproduct(A, B):
    s = 0
    for k in range(n):
        s += A[k] * B[k]
    return s


leadingzerocounts = [0] * m
for S in itertools.product([-1, 1], repeat=n - 1):
    S1 = S + (1,)
    S1S1 = S1 * 2
    for C in itertools.combinations(range(n - 1), n // 2):
        F = list(S1)
        for i in C:
            F[i] *= -1
        leadingzerocounts[0] += 4
        for i in range(1, m):
            if innerproduct(F, S1S1[i:i + n]):
                break
            leadingzerocounts[i] += 4
print(leadingzerocounts)

C版本，以了解我们的性能在多大程度上输给了PyPy ( PyPy的16大致相当于C的18 )：

代码语言：javascript
复制
#include <stdio.h>

enum {
  HALFN = 9,
  N = 2 * HALFN
};

int main(void) {
  long long lzc[N + 1];
  for (int i = 0; i < N + 1; i++) lzc[i] = 0;
  unsigned int xor = 1 << (N - 1);
  while (xor-- > 0) {
    if (__builtin_popcount(xor) != HALFN) continue;
    unsigned int s = 1 << (N - 1);
    while (s-- > 0) {
      lzc[0]++;
      unsigned int f = xor ^ s;
      for (int i = 1; i < N + 1; i++) {
        f = ((f & 1) << (N - 1)) | (f >> 1);
        if (__builtin_popcount(f ^ s) != HALFN) break;
        lzc[i]++;
      }
    }
  }
  for (int i = 0; i < N + 1; i++) printf(" %lld", 4 * lzc[i]);
  putchar('\n');
  return 0;
}

这个算法的并行性令人尴尬，因为它只是累加xor的所有值。对于C版本，粗略的计算表明，几千个小时的CPU时间足以计算n = 26，按照EC2上的当前汇率，这相当于几百美元。毫无疑问，有一些优化需要进行(例如，矢量化)，但对于这样的一次性优化，我不确定程序员的更多努力是值得的。

票数 22
EN

Stack Overflow用户

发布于 2014-11-17 03:14:46

提高n倍的一个非常简单的方法是更改以下代码：

代码语言：javascript
复制
def innerproduct(A, B):
    assert (len(A) == len(B))
    for j in xrange(len(A)):
        s = 0 
        for k in xrange(0,n):
            s+=A[k]*B[k]
    return s

至

代码语言：javascript
复制
def innerproduct(A, B):
    assert (len(A) == len(B))
    s = 0 
    for k in xrange(0,n):
        s+=A[k]*B[k]
    return s

(我不知道为什么在j上有循环，但它每次都进行相同的计算，所以没有必要。)

票数 7
EN

Stack Overflow用户

发布于 2014-11-20 02:47:15

我已经尝试将它转换到NumPy数组中，并借用了这个问题：itertools product speed up

这就是我所得到的(这里可能有更多的加速)：

代码语言：javascript
复制
def find_leading_zeros(n):
    if n % 2:
        return numpy.zeros(n)
    m = n+1
    leading_zero_counts = numpy.zeros(m)
    product_list = [-1, 1]
    repeat = n
    s = (numpy.array(product_list)[numpy.rollaxis(numpy.indices((len(product_list),) * repeat),
                                                  0, repeat + 1).reshape(-1, repeat)]).astype('int8')
    i = 0
    size = s.shape[0] / 2
    products = numpy.zeros((size, size), dtype=bool)
    while i < m:
        products += (numpy.tensordot(s[0:size, 0:size],
                                     numpy.roll(s, i, axis=1)[0:size, 0:size],
                                     axes=(-1,-1))).astype('bool')
        leading_zero_counts[i] = (products.size - numpy.sum(products)) * 4
        i += 1

    return leading_zero_counts

为n=14运行时，我得到：

代码语言：javascript
复制
>>> find_leading_zeros(14)
array([ 56229888.,  23557248.,   9903104.,   4160640.,   1758240.,
        755392.,    344800.,    172320.,    101312.,     75776.,
        65696.,     61216.,     59200.,     59200.,     59200.])

所以一切看起来都很好。至于速度：

代码语言：javascript
复制
>>> timeit.timeit("find_leading_zeros_old(10)", number=10)
28.775046825408936
>>> timeit.timeit("find_leading_zeros(10)", number=10)
2.236745834350586

看看你怎么想的。

编辑：

最初的版本使用了2074MB的内存用于N=14，所以我删除了连接数组，转而使用numpy.roll。此外，将数据类型更改为使用布尔数组，将使n=14的内存降至277MB。

在时间上，编辑速度又快了一点：

代码语言：javascript
复制
>>> timeit.timeit("find_leading_zeros(10)", number=10)
1.3816070556640625

EDIT2：

加上David指出的对称性，我再减少一次。它现在使用213MB。与之前的编辑相比的比较时间：

代码语言：javascript
复制
>>> timeit.timeit("find_leading_zeros(10)", number=10)
0.35357093811035156

我现在可以在我的mac上在14秒内完成n=14，我认为这对“纯python”来说是不错的。

票数 2
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：

https://stackoverflow.com/questions/26960749

复制
相关文章

相似问题

领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2024 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 深公网安备号 44030502008569

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号 | 京公网安备号11010802020287

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档