Как рассчитать очень большую корреляционную матрицу

Question

У меня есть np.array наблюдений z, где z.shape равно (100000, 60).Я хочу эффективно рассчитать корреляционную матрицу 100000x100000, а затем записать на диск координаты и значения только этих элементов> 0,95 (это очень малая часть от общего числа).

Моя версия грубой силы выглядит следующим образом, но, что неудивительно, очень медленно:

for i1 in range(z.shape[0]):
    for i2 in range(i1+1):
        r = np.corrcoef(z[i1,:],z[i2,:])[0,1]
        if r > 0.95:
            file.write("%6d %6d %.3f\n" % (i1,i2,r))

Я понимаю, что сама матрица корреляции может быть вычислена намного более эффективно в одномоперация с использованием np.corrcoef (z), но тогда требования к памяти огромны.Я также знаю, что можно разбить набор данных на блоки и вычислить подразделы размера укуса матрицы корреляции за один раз, но программирование этого и отслеживание индексов кажется излишне сложным.

Есть лидругой способ (например, использование memmap или pytables), который прост в кодировании и не предъявляет чрезмерных требований к физической памяти?

Answer 1

Вы можете использовать scipy.spatial.distance.pdist с metric = correlation, чтобы получить все корреляции без симметричных членов.К сожалению, это все равно оставит вам около 5e10 терминов, которые, вероятно, переполнят вашу память.

Вы можете попробовать переформулировать KDTree (который теоретически может обрабатывать косинусное расстояние и, следовательно, корреляционное расстояние)фильтровать для более высоких корреляций, но с 60 измерениями вряд ли вы бы сильно ускорились. Проклятие размерности - отстой.

Лучше всего, вероятно, использовать грубое форсирование блоков данных с использованием scipy.spatial.distance.cdist(..., metric = correlation), а затем сохранять только высокие корреляции в каждом блоке.Как только вы узнаете, насколько большой блок может обрабатывать ваша память без замедления из-за архитектуры памяти вашего компьютера, он должен быть намного быстрее, чем делать по одному за раз.

Answer 2

Поэкспериментировав с решением memmap, предложенным другими, я обнаружил, что, хотя это было быстрее, чем мой первоначальный подход (который занимал около 4 дней на моем Macbook), он все еще занимал очень много времени (по крайней мере, один день) -предположительно из-за неэффективной поэлементной записи в выходной файл.Это было неприемлемо, учитывая то, что мне нужно было выполнять вычисления несколько раз.

В конце концов, лучшее решение (для меня) - войти на портал Amazon Web Services EC2, создать экземпляр виртуальной машины (запускс изображением, оснащенным Anaconda Python) с 120+ ГБ ОЗУ, загрузите файл входных данных и выполните вычисления (используя метод умножения матриц) полностью в памяти ядра.Это завершилось примерно за две минуты!

Для справки, код, который я использовал, был в основном такой:

import numpy as np
import pickle
import h5py

# read nparray, dimensions (102000, 60)

infile = open(r'file.dat', 'rb')
x = pickle.load(infile)
infile.close()     

# z-normalize the data -- first compute means and standard deviations
xave = np.average(x,axis=1)
xstd = np.std(x,axis=1)

# transpose for the sake of broadcasting (doesn't seem to work otherwise!)
ztrans = x.T - xave
ztrans /= xstd

# transpose back
z = ztrans.T

# compute correlation matrix - shape = (102000, 102000)
arr = np.matmul(z, z.T)   
arr /= z.shape[0]

# output to HDF5 file
with h5py.File('correlation_matrix.h5', 'w') as hf:
        hf.create_dataset("correlation",  data=arr)

Answer 3

Исходя из моих грубых расчетов, вам нужна корреляционная матрица, содержащая 100 000 ^ 2 элементов.Это занимает около 40 ГБ памяти, при условии плавающих значений.
Это, вероятно, не поместится в памяти компьютера, в противном случае вы можете просто использовать corrcoef.Есть причудливый подход, основанный на собственных векторах, который я не могу найти прямо сейчас, и который попадает в (обязательно) сложную категорию ... Вместо этого полагайтесь на тот факт, что для данных с нулевым средним значение ковариация может быть найдена с использованием точечного произведения.

z0 = z - mean(z, 1)[:, None]
cov = dot(z0, z0.T)
cov /= z.shape[-1]

И это можно превратить в корреляцию путем нормализации по отклонениям

sigma = std(z, 1)
corr = cov
corr /= sigma
corr /= sigma[:, None]

Конечно, использование памяти все еще остается проблемой.Вы можете обойти это с массивами, отображаемыми в память (убедитесь, что он открыт для чтения и записи) и параметром out, равным dot (другой пример см. Оптимизация кода больших данных с небольшим объемом оперативной памяти )

N = z.shape[0]
arr = np.memmap('corr_memmap.dat', dtype='float32', mode='w+', shape=(N,N)) 
dot(z0, z0.T, out=arr)
arr /= sigma
arr /= sigma[:, None]

Затем вы можете перебрать полученный массив и найти индексы с большим коэффициентом корреляции.(Вы можете найти их непосредственно с помощью where(arr > 0.95), но сравнение создаст очень большой логический массив, который может или не может поместиться в памяти).