Эффективное вычисление коэффициента корреляции по столбцам

Оригинальный вопрос

Я сопоставляю строку P размера n с каждым столбцом матрицы O размера n×m. Я создал следующий код:

import numpy as np

def ColumnWiseCorrcoef(O, P):
    n = P.size
    DO = O - (np.sum(O, 0) / np.double(n))
    DP = P - (np.sum(P) / np.double(n))
    return np.dot(DP, DO) / np.sqrt(np.sum(DO ** 2, 0) * np.sum(DP ** 2))

Это более эффективно, чем наивный подход:

def ColumnWiseCorrcoefNaive(O, P):
    return np.corrcoef(P,O.T)[0,1:O[0].size+1]

Вот время, которое я получаю с numpy-1.7.1-MKL на ядре Intel:

O = np.reshape(np.random.rand(100000), (1000,100))
P = np.random.rand(1000)

%timeit -n 1000 A = ColumnWiseCorrcoef(O, P)
1000 loops, best of 3: 787 us per loop
%timeit -n 1000 B = ColumnWiseCorrcoefNaive(O, P)
1000 loops, best of 3: 2.96 ms per loop

Теперь вопрос: можете ли вы предложить еще более быструю версию кода для этой проблемы? Выдавливание дополнительных 20% было бы здорово.

ОБНОВЛЕНИЕ май 2017

Через некоторое время я вернулся к этой проблеме, перезапустил и расширил задачу и тесты.

  1. Используя einsum, я расширил код до случая, когда P - это не строка, а матрица. Таким образом, задача состоит в том, чтобы соотнести все столбцы O со всеми столбцами P.

  2. Мне было любопытно, как можно решить одну и ту же проблему на разных языках, обычно используемых для научных вычислений. Я реализовал ее (с помощью других людей) в MATLAB, Julia, и R. MATLAB и Julia - самые быстрые, и у них есть специальная программа. вычислить столбцовую корреляцию. У R также есть специальная рутина, но она самая медленная.

  3. В текущей версии numpy (1.12.1 от Anaconda) einsum по-прежнему побеждает выделенные функции, которые я использовал.

Все сценарии и сроки доступны по адресу https://github.com/ikizhvatov/efficient-columnwise-correlation.

Источник

2 ответа

Решение

Мы можем представить np.einsum к этому; однако, ваш пробег может варьироваться в зависимости от вашей компиляции и от того, использует ли он SSE2 или нет. Экстра einsum запросы на замену операций суммирования могут показаться посторонними, но numpy ufuncs не использует SSE2, пока numpy 1.8 пока einsum делает, и мы можем избежать нескольких if заявления.

Запустив это на ядре Opteron с Intel MKL Blas я получаю странный результат, как я ожидал dot позвоните, чтобы занять большую часть времени.

def newColumnWiseCorrcoef(O, P):
    n = P.size
    DO = O - (np.einsum('ij->j',O) / np.double(n))
    P -= (np.einsum('i->',P) / np.double(n))
    tmp = np.einsum('ij,ij->j',DO,DO)
    tmp *= np.einsum('i,i->',P,P)          #Dot or vdot doesnt really change much.
    return np.dot(P, DO) / np.sqrt(tmp)

O = np.reshape(np.random.rand(100000), (1000,100))
P = np.random.rand(1000)

old = ColumnWiseCorrcoef(O,P)
new = newColumnWiseCorrcoef(O,P)

np.allclose(old,new)
True

%timeit ColumnWiseCorrcoef(O,P)
100 loops, best of 3: 1.52ms per loop

%timeit newColumnWiseCorrcoef(O,P)
1000 loops, best of 3: 518us per loop

Запустив это снова с системой intel с intel mkl, я получаю нечто более разумное / ожидаемое:

%timeit ColumnWiseCorrcoef(O,P)
1000 loops, best of 3: 524 µs per loop

%timeit newColumnWiseCorrcoef(O,P)
1000 loops, best of 3: 354 µs per loop

Снова на машине Intel с чем-то немного большим:

O = np.random.rand(1E5,1E3)
P = np.random.rand(1E5)

%timeit ColumnWiseCorrcoef(O,P)
1 loops, best of 3: 1.33 s per loop

%timeit newColumnWiseCorrcoef(O,P)
1 loops, best of 3: 791 ms per loop
Источник

Это менее эффективное решение, но оно легко читаемо и состоит всего из одной строки кода:

np.einsum("ij,i->j", spst.zscore(A), spst.zscore(B)) / A.shape[0]

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

np.einsum("ij,ij->j", spst.zscore(A), spst.zscore(B)) / A.shape[0]

Согласно моему тесту синхронизации на моем iMac с четырехъядерным процессором Intel Core i5 с тактовой частотой 2,8 ГГц, он примерно в 3 раза медленнее:

      import scipy.stats as spst

def ColumnWiseCorrcoef(O, P):
    n = P.size
    DO = O - (np.sum(O, 0) / np.double(n))
    DP = P - (np.sum(P) / np.double(n))
    return np.dot(DP, DO) / np.sqrt(np.sum(DO ** 2, 0) * np.sum(DP ** 2))

def newColumnWiseCorrcoef(O, P):
    n = P.size
    DO = O - (np.einsum('ij->j',O) / np.double(n))
    P -= (np.einsum('i->',P) / np.double(n))
    tmp = np.einsum('ij,ij->j',DO,DO)
    tmp *= np.einsum('i,i->',P,P)          #Dot or vdot doesnt really change much.
    return np.dot(P, DO) / np.sqrt(tmp)

def OneLineColumnWiseCorrcoef(A, B):
    return np.einsum("ij,i->j", spst.zscore(A), spst.zscore(B)) / A.shape[0]

O = np.reshape(np.random.rand(100000), (1000,100))
P = np.random.rand(1000)

old = ColumnWiseCorrcoef(O,P)
new = OneLineColumnWiseCorrcoef(O, P)
np.allclose(old,new)
True

%timeit ColumnWiseCorrcoef(O,P)
# 325 µs ± 19.8 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

%timeit newColumnWiseCorrcoef(O,P)
# 274 µs ± 713 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

%timeit OneLineColumnWiseCorrcoef(O,P)
# 1.02 ms ± 4.4 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
Источник
Другие вопросы по тегам