Вычисление спектральных норм эрмитовых матриц ~1m: `numpy.linalg.norm` выполняется слишком медленно. ⇐ Python

[Anonymous]

Я хотел бы вычислить спектральные нормы эрмитовых матриц N 8x8, при этом N близок к 1E6. В качестве примера возьмем 1 миллион случайных комплексных матриц размером 8x8:

Код: Выделить всё

import numpy as np

array = np.random.rand(8,8,1e6)  + 1j*np.random.rand(8,8,1e6)

Сейчас использование numpy.linalg.norm занимает у меня почти 10 секунд:

Код: Выделить всё

np.linalg.norm(array, ord=2, axis=(0,1))

Я попробовал использовать приведенный ниже код Cython, но это дало лишь незначительное улучшение производительности:

Код: Выделить всё

import numpy as np
cimport numpy as np
cimport cython

np.import_array()

DTYPE = np.complex64

@cython.boundscheck(False)
@cython.wraparound(False)
def function(np.ndarray[np.complex64_t, ndim=3] Array):
assert Array.dtype == DTYPE
cdef int shape0 = Array.shape[2]
cdef np.ndarray[np.float32_t, ndim=1] normarray = np.zeros(shape0, dtype=np.float32)
normarray = np.linalg.norm(Array, ord=2, axis=(0, 1))
return normarray

Я также попробовал numba и некоторые другие функции scipy (например, scipy.linalg.svdvals) для вычисления сингулярных значений этих матриц. Все по-прежнему слишком медленно.

Нельзя ли сделать это быстрее? numpy уже оптимизирован до такой степени, что невозможно получить прирост скорости при использовании Cython или numba? Или мой код крайне неэффективен и я делаю что-то в корне неправильно?

Я заметил, что только два ядра моего процессора используются на 100% при выполнении вычислений. Имея это в виду, я рассмотрел предыдущие вопросы StackOverflow:

• почему numpy.mean не является многопоточным?
• Почему многопроцессорная обработка использует только одно ядро ​​после импорта numpy?
• многопоточный blas в Python /numpy (не помогло)

и несколько других, но, к сожалению, у меня до сих пор нет решения.

Я подумал о том, чтобы разделить массив на более мелкие фрагменты и обработать их параллельно (возможно, на графическом процессоре с использованием CUDA). Есть ли способ сделать это в numpy/Python? Я пока не знаю, где находится узкое место в моем коде, т. е. связано ли оно с процессором, с памятью или, возможно, с чем-то другим.

Подробнее здесь: https://stackoverflow.com/questions/336 ... -norm-is-t