Многопоточное целочисленное матричное умножение в NumPy/SciPy

Выполнение чего-то типа

import numpy as np
a = np.random.rand(10**4, 10**4)
b = np.dot(a, a)

использует несколько ядер, и он работает хорошо.

Элементы в a, однако, являются 64-битными поплавками (или 32-разрядными в 32-разрядных платформах?), и я хотел бы умножить 8-разрядные целые массивы. Однако попробуйте следующее:

a = np.random.randint(2, size=(n, n)).astype(np.int8)

приводит к точечному продукту, не использующему несколько ядер, и, таким образом, на моем ПК медленнее на 1000 раз.

array: np.random.randint(2, size=shape).astype(dtype)

dtype    shape          %time (average)

float32 (2000, 2000)    62.5 ms
float32 (3000, 3000)    219 ms
float32 (4000, 4000)    328 ms
float32 (10000, 10000)  4.09 s

int8    (2000, 2000)    13 seconds
int8    (3000, 3000)    3min 26s
int8    (4000, 4000)    12min 20s
int8    (10000, 10000)  It didn't finish in 6 hours

float16 (2000, 2000)    2min 25s
float16 (3000, 3000)    Not tested
float16 (4000, 4000)    Not tested
float16 (10000, 10000)  Not tested

Я понимаю, что NumPy использует BLAS, который не поддерживает целые числа, но если я использую обертки SciPy BLAS, т.е.

import scipy.linalg.blas as blas
a = np.random.randint(2, size=(n, n)).astype(np.int8)
b = blas.sgemm(alpha=1.0, a=a, b=a)

вычисление многопоточное. Теперь blas.sgemm работает с точно таким же временем, как np.dot для float32, но для не-float он преобразует все в float32 и выводит поплавки, чего нет np.dot. (Кроме того, b теперь находится в F_CONTIGUOUS порядке, что является меньшей проблемой).

Итак, если я хочу сделать умножение целых матриц, я должен выполнить одно из следующих действий:

Используйте NumPy больно медленно np.dot и радуйся, что я могу сохранить 8-битные целые числа.
Используйте SciPy sgemm и используйте 4-кратную память.
Используйте Numpy np.float16 и используйте только память 2x, с предостережением, что np.dot намного медленнее на массивах float16, чем на массивах float32, больше, чем int8.
Найти оптимизированную библиотеку для многопоточного умножения целых матриц (на самом деле, Mathematica делает это, но я бы предпочел решение Python), идеально поддерживая 1-битные массивы, хотя 8-битные массивы тоже отлично... (я на самом деле намерен сделать умножение матриц над конечным полем Z/2Z, и я знаю, что могу сделать это с помощью Sage, что довольно Pythonic, но, опять же, есть ли что-то строго Python?)

Могу ли я выполнить опцию 4? Существует ли такая библиотека?

Отказ от ответственности: я фактически запускаю NumPy + MKL, но я пробовал аналогичный тест на vanilly NumPy с аналогичными результатами.

Ответ 1

Вариант 5 - Сбросьте пользовательское решение: Разделите матричный продукт в нескольких субпродуктах и выполните их параллельно. Это можно относительно легко реализовать со стандартными модулями Python. Подпрограммы вычисляются с помощью numpy.dot, который освобождает глобальную блокировку интерпретатора. Таким образом, можно использовать потоки, которые относительно легки и могут обращаться к массивам из основного потока для эффективности памяти.

Реализация:

import numpy as np
from numpy.testing import assert_array_equal
import threading
from time import time


def blockshaped(arr, nrows, ncols):
    """
    Return an array of shape (nrows, ncols, n, m) where
    n * nrows, m * ncols = arr.shape.
    This should be a view of the original array.
    """
    h, w = arr.shape
    n, m = h // nrows, w // ncols
    return arr.reshape(nrows, n, ncols, m).swapaxes(1, 2)


def do_dot(a, b, out):
    #np.dot(a, b, out)  # does not work. maybe because out is not C-contiguous?
    out[:] = np.dot(a, b)  # less efficient because the output is stored in a temporary array?


def pardot(a, b, nblocks, mblocks, dot_func=do_dot):
    """
    Return the matrix product a * b.
    The product is split into nblocks * mblocks partitions that are performed
    in parallel threads.
    """
    n_jobs = nblocks * mblocks
    print('running {} jobs in parallel'.format(n_jobs))

    out = np.empty((a.shape[0], b.shape[1]), dtype=a.dtype)

    out_blocks = blockshaped(out, nblocks, mblocks)
    a_blocks = blockshaped(a, nblocks, 1)
    b_blocks = blockshaped(b, 1, mblocks)

    threads = []
    for i in range(nblocks):
        for j in range(mblocks):
            th = threading.Thread(target=dot_func, 
                                  args=(a_blocks[i, 0, :, :], 
                                        b_blocks[0, j, :, :], 
                                        out_blocks[i, j, :, :]))
            th.start()
            threads.append(th)

    for th in threads:
        th.join()

    return out


if __name__ == '__main__':
    a = np.ones((4, 3), dtype=int)
    b = np.arange(18, dtype=int).reshape(3, 6)
    assert_array_equal(pardot(a, b, 2, 2), np.dot(a, b))

    a = np.random.randn(1500, 1500).astype(int)

    start = time()
    pardot(a, a, 2, 4)
    time_par = time() - start
    print('pardot: {:.2f} seconds taken'.format(time_par))

    start = time()
    np.dot(a, a)
    time_dot = time() - start
    print('np.dot: {:.2f} seconds taken'.format(time_dot))

В этой реализации я получаю ускорение примерно x4, что является физическим числом ядер в моей машине:

running 8 jobs in parallel
pardot: 5.45 seconds taken
np.dot: 22.30 seconds taken