Как ускорить несколько внутренних продуктов в python

Ответ 1

Этот новый код получает ускорение по порядку, используя циклическую симметрию проблемы. Эта версия Python перечисляет ожерелья с алгоритмом Дюваля; версия C использует грубую силу. Оба включают ускорения, описанные ниже. На моей машине версия C решает n = 20 за 100 секунд! Расчет по сравнению с конвертом предполагает, что если вы позволите ему работать в течение недели на одном ядре, мог бы сделать n = 26, и, как отмечено ниже, он поддается значению parallelism.

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:

#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).

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 (16 для PyPy примерно эквивалентно 18 для C):

#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 вычисление задней части конверта позволяет предположить, что для вычисления n = 26 потребуется несколько тысяч часов процессорного времени, что составляет несколько сотен долларов при текущих ставках на EC2. Несомненно, есть некоторые оптимизации, которые нужно сделать (например, векторизация), но для одноразового использования я не уверен, сколько стоит больше усилий программиста.

Ответ 2

Одной из очень простых скоростей с коэффициентом n является изменение этого кода:

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

к

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, но он просто делает один и тот же расчет каждый раз, поэтому не нужен.)

Ответ 3

У меня есть передача, передающая это в массивы NumPy и заимствованные из этого вопроса: скорость продукта itertools

Это то, что у меня есть (здесь может быть больше ускорений):

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 Я получаю:

>>> find_leading_zeros(14)
array([ 56229888.,  23557248.,   9903104.,   4160640.,   1758240.,
        755392.,    344800.,    172320.,    101312.,     75776.,
        65696.,     61216.,     59200.,     59200.,     59200.])

Итак, все выглядит хорошо. Что касается скорости:

>>> timeit.timeit("find_leading_zeros_old(10)", number=10)
28.775046825408936
>>> timeit.timeit("find_leading_zeros(10)", number=10)
2.236745834350586

Посмотрите, что вы думаете.

EDIT:

В исходной версии используется 2074 МБ памяти для N = 14, поэтому я удалил конкатенированный массив и вместо этого использовал numpy.roll. Также, изменяя типы данных для использования логического массива, память уменьшается до 277 МБ для n = 14.

Время, когда редактирование будет немного быстрее:

>>> timeit.timeit("find_leading_zeros(10)", number=10)
1.3816070556640625

EDIT2:

Итак, добавив в симметрию, как указал Дэвид, я уменьшаю это снова. Теперь он использует 213 МБ. Сравнение времени сравнения по сравнению с предыдущими изменениями:

>>> timeit.timeit("find_leading_zeros(10)", number=10)
0.35357093811035156 

Теперь я могу сделать случай n = 14 за 14 секунд в моей книге mac, что неплохо для "чистого python", я думаю.

Ответ 4

Я пытался ускорить это, и я потерпел неудачу плохо:( Но я отправляю код, как-то быстрее, но не достаточно быстро для таких значений, как n=24.

Мои предположения

Ваши списки состоят из значений, поэтому я решил использовать числа вместо списков - каждый бит представляет одно из возможных значений: если бит установлен, то это означает 1, если он обнулен, это означает -1. Единственным возможным результатом умножения {-1, 1} является 1 или -1, поэтому вместо умножения я использовал побитовое XOR. Я также заметил, что есть симметрия, поэтому вам нужно только проверить подмножество (одну четверть) возможных списков и умножить результат на 4 (Дэвид объяснил это в своем ответе).

Наконец, я поместил результаты возможных операций в таблицы, чтобы устранить необходимость вычислений. Это занимает много памяти, но кому это нужно (для n=24 было около 150 МБ)?

И тогда @David Эйзенстат ответил на вопрос:) Итак, я взял его код и изменил его на основе бит. Это примерно в 2-3 раза быстрее (для n=16 потребовалось ~ 30 секунд, по сравнению с ~ 90 решений David), но я думаю, что этого еще недостаточно, чтобы получить результаты для n=26 или около того.

import itertools

n = 16
m = n + 1
mask = (2 ** n) - 1

# Create table of sum results (replaces innerproduct())
tab = []
for a in range(2 ** n):
    s = 0
    for k in range(n):
        s += -1 if a & 1 else 1
        a >>= 1
    tab.append(s)

# Create combination bit masks for combinations
comb = []
for C in itertools.combinations(range(n - 1), n // 2):
    xor = 0
    for i in C:
       xor |= (1 << i)
    comb.append(xor)

leadingzerocounts = [0] * m
for S in xrange(2 ** (n-1)):
    S1 = S + (1 << (n-1))
    S1S1 = S1 + (S1 << n)

    for xor in comb:
        F = S1 ^ xor

        leadingzerocounts[0] += 4
        for i in range(1, m):
            if tab[F ^ ((S1S1 >> i) & mask)]:
                break
            leadingzerocounts[i] += 4

print(leadingzerocounts)

Выводы

Я думал, что придумал что-то блестящее и надеялся, что все это беспорядок с битами даст большой прирост скорости, но толчок был разочаровывающим: (

Я думаю, что причина в том, как Python использует операторы - он вызывает функцию для каждой арифметической (или логической) операции, даже если это может быть сделано с помощью одной команды ассемблера (я надеялся, что pypy сможет упростить операции для этого уровень, но это не так). Таким образом, возможно, если бы C (или ASM) использовалось с этим решением для работы с битами, оно бы отлично работало (возможно, вы могли бы добраться до n=24).

Ответ 5

По моему мнению, хороший способ повысить производительность - использовать встроенные функции python.

Сначала используется карта для расчета произведения записей:

>>> a =[1,2,3]
>>> b = [4,5,6]
>>>map(lambda x,y : x*y, a , b)
[4, 10, 18]

Затем используйте reduce для вычисления сумм:

>>> reduce(lambda v,w: v+w, map(lambda x,y :x*y, a, b))
32

Итак, ваша функция станет

def innerproduct(A, B):
    assert (len(A) == len(B))
    return reduce(lambda v,w: v+w, map(lambda x,y :x*y, A, B))

Затем мы можем взять все эти "для циклов" и заменить их генераторами и поймать StopIteration.

#!/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))
    return reduce(lambda v,w: v+w, map(lambda x,y :x*y, A, B))


leadingzerocounts = [0]*m

S_gen = itertools.product([-1,1], repeat = n)

try:
    while(True):
       S = S_gen.next()
       S1 = S + S
       F_gen = itertools.product([-1,1], repeat = n)
       try:
           while(True):
               F = F_gen.next()
               for i in xrange(m):
                   ip = innerproduct(F, S1[i:i+n])
                   if (ip == 0):
                       leadingzerocounts[i] +=1
                       i+=1
                   else:
                      break
       except StopIteration:
           pass

except StopIteration as e:
    print e

print leadingzerocounts

Я наблюдал ускорение для меньшего n, но у моей jalopy не хватало лошадиных сил для вычисления моей версии или исходного кода для n = 14. Еще одним способом ускорить это было бы замещение строки:

    F_gen = itertools.product([-1,1], repeat = n)