Выполнение матричной операции над двумя большими матрицами

Question

У меня есть две большие матрицы (40000 * 4096), и я хотел бы сравнить и сопоставить каждую строку первой матрицы со всеми строками для второй матрицы, и в результате выход будет иметь размер (40000 *40000).Тем не менее, поскольку мне нужно делать это несколько тысяч раз, это невероятное время, которое занимает 26 тысяч секунд на каждую итерацию, то есть 5000 раз ... Я был бы рад, если бы вы могли дать мне несколько умных предложений.Спасибо.PS это то, что я сделал до сих пор только для одной итерации (1 из 5000)

def matcher(Antigens, Antibodies,ind):
    temp = np.zeros((Antibodies.shape[0],Antibodies.shape[1]))
    output = np.zeros((Antibodies.shape[0],1))
    for i in range(len(Antibodies)):
        temp[i] = np.int32(np.equal(Antigens[ind],Antibodies[i]))
        output[i] = np.sum(temp[i])
    return output
output = [matcher(gens,Antibodies) for gens in Antigens]

Answer 1

Хорошо, я думаю, что понимаю, какова ваша цель:

Подсчет количества совпадений строк (матрица антиген против антитела).Каждая строка результирующего вектора (40000 x 1) представляет количество точных совпадений между 1 строкой антигена и всей строкой антител (поэтому значения от 0 до 40_000).

Я сделал несколько поддельных данных:

import numpy as np
import numba as nb

num_mat = 5       # number of matrices
num_row = 10_000  # number of rows per matrix
num_elm = 4_096   # number of elements per row
dim = (num_mat,num_row,num_elm)

Antigens = np.random.randint(0,256,dim,dtype=np.uint8)
Antibodies = np.random.randint(0,256,dim,dtype=np.uint8)

Здесь есть один важный момент: я уменьшил матрицы до наименьшего типа данных, который может представлять данные, чтобы уменьшить объем их памяти.Я не уверен, как выглядят ваши данные, но, надеюсь, вы тоже можете это сделать.

Кроме того, следующий код предполагает, что ваши измерения выглядят фальшивыми данными:

(числоматриц, строк, элементов)

@nb.njit
def match_arr(arr1, arr2):
    for i in range(arr1.shape[0]): #4096 vs 4096
        if arr1[i] != arr2[i]:
            return False
    return True

@nb.njit
def match_mat_sum(ag, ab):
    out = np.zeros((ag.shape[0])) # 40000
    for i in range(ag.shape[0]):
        tmp = 0
        for j in range(ab.shape[0]):
            tmp += match_arr(ag[i], ab[j])
        out[i] = tmp
    return out

@nb.njit(parallel=True)
def match_sets(Antigens, Antibodies):
    out = np.empty((Antigens.shape[0] * Antibodies.shape[0], Antigens.shape[1])) # 5000 x 40000
    # multiprocessing per antigen matrix, may want to move this as suits your data
    for i in nb.prange(Antigens.shape[0]):
        for j in range(Antibodies.shape[0]):
            out[j+(5*i)] = match_mat_sum(Antigens[i], Antibodies[j]) # need to figure out the index to avoid race conditions
    return out

Я сильно опираюсь на Нумбу.Одна из ключевых оптимизаций заключается не в проверке эквивалентности целых строк с помощью np.equal(), а в написании пользовательской функции match_arr(), которая прерывается, как только находит несоответствующий элемент.Надеюсь, это позволит нам пропустить кучу сравнений.

Сравнение времени:

%timeit match_arr(arr1, arr2)
314 ns ± 0.361 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

%timeit np.equal(arr1, arr2)
1.07 µs ± 5.35 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

match_mat_sum

Эта функция просто вычисляет средний шаг (40 000 x 1вектор), который представляет сумму точных совпадений между двумя матрицами.Этот шаг уменьшает две матрицы, такие как: (mxn), (oxn) -> (m)

match_sets()

Последняя функция распараллеливает эту операцию с явными параллельными цикламидо nb.prange.Возможно, вы захотите переместить эту функцию в другой цикл в зависимости от того, как выглядят ваши данные (например, если у вас одна матрица антигена, но 5000 матриц антител, вам следует переместиться на prange во внутреннюю петлю, или вы не будете использовать параллелизацию).Поддельные данные предполагают наличие антигена и матрицы антител.

Еще одна важная вещь, которую следует здесь отметить, это индексирование массива out.Чтобы избежать условий гонки, каждый явный цикл должен записывать в уникальное пространство.Опять же, в зависимости от ваших данных, вам нужно будет индексировать правильное «место», чтобы поместить результат.

На Ryzen 1600 (6-ядерном) с 16 гигабайтами оперативной памяти, используя эти поддельные данные, ясгенерировал результат за 10,2 секунды.

Ваши данные примерно в 3200 раз больше.Предполагая линейное масштабирование, полный набор займет около 9 часов, при условии, что у вас достаточно памяти.

Вы также можете написать некоторый пакетный загрузчик, вместо того, чтобы загружать 5000 гигантских матриц непосредственно в память.

Answer 2

Эту проблему можно решить с помощью смешанного широковещательного вещания и модуля Numberxpr, который выполняет операции быстро, минимизируя хранение промежуточных значений.

import numexpr as ne

# expand arrays dimensions to support broadcasting when doing comparison
Antigens, Antibodies = Antigens[None, :, :], Antibodies[:, None, :]
output = ne.evaluate('sum((Antigens==Antibodies)*1, axis=2)')
# *1 is a hack because numexpr does not currently support sum on bool

Это может быть быстрее, чем ваше текущее решение, нодля таких больших массивов это займет некоторое время.

Производительность Numberxpr для этих операций немного слабовата, но вы можете по крайней мере использовать трансляцию внутри цикла:

output = np.zeros((Antibodies.shape[0],)*2, dtype=np.int32)
for row, out_row in zip(Antibodies, output):
    (row[None,:]==Antigens).sum(1, out=out_row)