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

Question

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

import numpy as np 

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

В настоящее время он занимает у меня почти 10 секунд с помощью numpy.linalg.norm:

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:

• why isn't numpy.mean multithreaded?

• Why does multiprocessing use only a single core after I import numpy?

• multithreaded blas in python/numpy (не помогло)

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

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

Answer 1

копания в код np.linalg.norm, я вывел, что для этих параметров, то найти максимум матрицы сингулярных значений по размерности N

Сначала приготовьте небольшой массив выборки.Сделать N первый размер, чтобы устранить rollaxis операцию:

In [268]: N=10; A1 = np.random.rand(N,8,8)+1j*np.random.rand(N,8,8) 

In [269]: np.linalg.norm(A1,ord=2,axis=(1,2)) 
Out[269]: 
array([ 5.87718306, 5.54662999, 6.15018125, 5.869058 , 5.80882818, 
     5.86060462, 6.04997992, 5.85681085, 5.71243196, 5.58533323]) 

эквивалент операции:

In [270]: np.amax(np.linalg.svd(A1,compute_uv=0),axis=-1) 
Out[270]: 
array([ 5.87718306, 5.54662999, 6.15018125, 5.869058 , 5.80882818, 
     5.86060462, 6.04997992, 5.85681085, 5.71243196, 5.58533323]) 

одинаковые значения, и то же самое время:

In [271]: timeit np.linalg.norm(A1,ord=2,axis=(1,2)) 
1000 loops, best of 3: 398 µs per loop 
In [272]: timeit np.amax(np.linalg.svd(A1,compute_uv=0),axis=-1) 
1000 loops, best of 3: 389 µs per loop 

И большую часть времени, проведенного в svd, который производит массив (N, 8):

In [273]: timeit np.linalg.svd(A1,compute_uv=0) 
1000 loops, best of 3: 366 µs per loop 

Так что, если вы хотите ускорить работу norm, вы можете более подробно изучить это значение svd. svd использует np.linalg._umath_linalg функции - это файл .so - скомпилирован.

c код в https://github.com/numpy/numpy/blob/97c35365beda55c6dead8c50df785eb857f843f0/numpy/linalg/umath_linalg.c.src

Он уверен, похоже, это самый быстрый вы получите. Нет петли уровня Python. Любой цикл в этом коде c или функция lapack, которую он вызывает.

Answer 2

np.linalg.norm(A, ord=2) вычисляет спектральную норму путем нахождения наибольшего сингулярного значения с использованием SVD. Однако, так как ваши 8х8 подматрицы Эрмитова, их наибольшие сингулярные значения будут равны максимуму их абсолютных собственных (see here):

import numpy as np 

def random_symmetric(N, k): 
    A = np.random.randn(N, k, k) 
    A += A.transpose(0, 2, 1) 
    return A 

N = 100000 
k = 8 
A = random_symmetric(N, k) 

norm1 = np.abs(np.linalg.eigvalsh(A)).max(1) 
norm2 = np.linalg.norm(A, ord=2, axis=(1, 2)) 

print(np.allclose(norm1, norm2)) 
# True 

Eigendecomposition на эрмитовой матрицы совсем немного быстрее, чем СВД:

In [1]: %%timeit A = random_symmetric(N, k) 
np.linalg.norm(A, ord=2, axis=(1, 2)) 
    ....: 
1 loops, best of 3: 1.54 s per loop 

In [2]: %%timeit A = random_symmetric(N, k) 
np.abs(np.linalg.eigvalsh(A)).max(1) 
    ....: 
1 loops, best of 3: 757 ms per loop