计算~1m埃尔米特矩阵的谱范数:“numpy.linalg.normam”太慢

Question

我想计算N8x8埃尔米特矩阵的谱范数，其中N接近1E6。以这100万个随机复数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无法提高速度的程度？还是我的代码效率很低，而且我做了一些根本错误的事情？

我注意到，在进行计算时，只有两个CPU内核被100%利用。考虑到这一点，我查看了之前的StackOverflow问题：

• 为什么numpy.mean不是多线程的？

• 为什么在我导入numpy之后，多处理只使用单个内核？

• python/numpy中的多线程blas（没有帮助）

和其他几个，但不幸的是，我仍然没有解决方案。

我考虑将我的数组分割成更小的块，并并行处理这些块（可能在GPU上使用CUDA）。在numpy/Python中有办法做到这一点吗？我还不知道我的代码中的瓶颈在哪里，也就是说，它是CPU还是内存绑定，或者可能是不同的东西。

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找到最大奇异值来计算谱范数。然而，由于8x8子矩阵是埃尔米特矩阵，它们的最大奇异值将等于它们的绝对本征值的最大值（见此处）：

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

Hermitian矩阵上的特征分解比SVD:快很多

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