Python加速奇异值分解

Question

我想计算3D矩阵的每个切片的奇异值分解。

我使用numpy和scipy来计算SVD，但它们都比MATLAB实现慢得多。numpy和scipy版本大约需要7秒，而MATLAB版本只需要0.7秒。

有没有一种方法可以加速Python中的SVD计算？

Python

import time
import scipy.linalg
import numpy.linalg
A = np.random.rand(100, 100, 1000) + 1j * np.random.rand(100, 100, 1000)
S = np.empty((A.shape[2], min(A.shape[0:1])))
t1 = time.time()
for i in range(A.shape[2]):
S[i, :] = numpy.linalg.svd(A[:, :, i], compute_uv=False)
print("[numpy] Elapsed time: {:.3f} s".format(time.time() - t1))
t1 = time.time()
for i in range(A.shape[2]):
S[i, :] = scipy.linalg.svdvals(A[:, :, i])
print("[scipy] Elapsed time: {:.3f} s".format(time.time() - t1))
# [numpy] Elapsed time: 7.137 s
# [scipy] Elapsed time: 7.435 s

MATLAB

A = randn(100, 100, 1000) + 1j * randn(100, 100, 1000);
S = nan(size(A,3), min(size(A, [1 2])));
tic;
for i = 1:size(A, 3)
S(i, :) = svd(A(:,:,i));
end
toc;
% Elapsed time is 0.702556 seconds.

这是np.show_config():的输出

blas_mkl_info:
NOT AVAILABLE
blis_info:
NOT AVAILABLE
openblas_info:
library_dirs = ['D:\a\1\s\numpy\build\openblas_info']
libraries = ['openblas_info']
language = f77
define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
library_dirs = ['D:\a\1\s\numpy\build\openblas_info']
libraries = ['openblas_info']
language = f77
define_macros = [('HAVE_CBLAS', None)]
lapack_mkl_info:
NOT AVAILABLE
openblas_lapack_info:
library_dirs = ['D:\a\1\s\numpy\build\openblas_lapack_info']
libraries = ['openblas_lapack_info']
language = f77
define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
library_dirs = ['D:\a\1\s\numpy\build\openblas_lapack_info']
libraries = ['openblas_lapack_info']
language = f77
define_macros = [('HAVE_CBLAS', None)]
Supported SIMD extensions in this NumPy install:
baseline = SSE,SSE2,SSE3
found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2,AVX512F,AVX512CD,AVX512_SKX
not found = AVX512_CLX,AVX512_CNL
None

Answer 1

我不确定这是否能解决您的问题，但您不需要循环任何内容，因为np.linalg.svd()已经处理了n维数组。(即使没有，你基本上也不需要NumPy中的循环。(

这就是我做你正在做的事情的方式：

import time
import numpy as np
import numpy.linalg as la
shape = (100, 100, 1000)
rng = np.random.default_rng(42)
A = rng.random(shape) + 1j * rng.random(shape)
t1 = time.perf_counter()
S = la.svd(A.T, compute_uv=False)
t2 = time.perf_counter()
print(f"[numpy] Elapsed time: {t2 - t1:.3f} s")

在我的电脑上，这需要0.83秒(使用"英特尔MKL"(。不过我没有尝试你的MATLAB代码，所以我不确定这是否会加快你的速度。

Answer 2

对于支持"英特尔数学内核库"(MKL(的计算机，安装使用MKL的NumPy/SciPy版本可以显著缩短计算时间。感谢@joni和@kwinkunks提供的信息

在我的例子中，计算时间从OpenBLAS的7秒减少到"英特尔MKL "的0.68秒

这可以通过根据以下教程通过源代码构建NumPy或SciPy来完成：从源构建NumPy/SciPy

或者，安装带有预构建MKL支持的NumPy和SciPy版本的Anaconda平台可以简化这一过程。