计算 ~1m Hermitian 矩阵的谱范数：“numpy.linalg.norm”太慢

Question

我想计算N个 8x8 Hermitian 矩阵的谱范数，其中N接近 1E6。以这 100 万个随机复数 8x8 矩阵为例：

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 无法提高速度的程度？或者我的代码效率非常低并且我正在做一些根本错误的事情？

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

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

• 为什么导入 numpy 后多处理仅使用单个核心？

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

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

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

Answer 1

深入研究代码np.linalg.norm，我推断出，对于这些参数，它正在寻找 N 维上矩阵奇异值的最大值

\n\n

首先生成一个小样本数组。制作N第一个维度来消除rollaxis操作：

\n\n

In [268]: N=10; A1 = np.random.rand(N,8,8)+1j*np.random.rand(N,8,8)\n\nIn [269]: np.linalg.norm(A1,ord=2,axis=(1,2))\nOut[269]: \narray([ 5.87718306,  5.54662999,  6.15018125,  5.869058  ,  5.80882818,\n        5.86060462,  6.04997992,  5.85681085,  5.71243196,  5.58533323])\n

\n\n

等效操作：

\n\n

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

\n\n

相同的值和相同的时间：

\n\n

In [271]: timeit np.linalg.norm(A1,ord=2,axis=(1,2))\n1000 loops, best of 3: 398 \xc2\xb5s per loop\nIn [272]: timeit np.amax(np.linalg.svd(A1,compute_uv=0),axis=-1)\n1000 loops, best of 3: 389 \xc2\xb5s per loop\n

\n\n

而大部分时间都花在svd，它生成一个 (N,8) 数组：

\n\n

In [273]: timeit np.linalg.svd(A1,compute_uv=0)\n1000 loops, best of 3: 366 \xc2\xb5s per loop\n

\n\n

因此，如果您想加快速度norm，则必须进一步研究加快速度svd。 svd使用 np.linalg._umath_linalg函数 - 即.so文件）。

\n\n

代码c在https://github.com/numpy/numpy/blob/97c35365beda55c6dead8c50df785eb857f843f0/numpy/linalg/umath_linalg.c.src

\n\n

看起来这确实是您能获得的最快速度。没有 Python 级别的循环。任何循环都在该c代码中，或者lapack调用的函数中。

\n