python - 快速生成numpy数组

Question

我想通过使用每个元素的行向量、列向量和值向量来生成一个稀疏的 numpy ndarray。

例如，如果我有

row_index=np.array([0,1,2])
column_index=np.array([2,1,0])
value=np.array([4,5,6])

然后我想要一个矩阵

[0,0,4
0,5,0
6,0,0]

scipy.sparsenumpy 中是否有一个函数可以通过 using执行类似的操作scipy.sparse.csc_matrix((data, (row_ind, col_ind)), [shape=(M, N)])？如果没有，有没有办法在没有for循环的情况下生成矩阵？我想加快代码但 scipy.sparse在计算过程中很慢，我想要的矩阵不是那么大。

标签： pythonarraysnumpy

Answer 1

如果您想要的矩阵不是很大，那么创建一个常规（非稀疏）ndarray 可能会更快。例如，您可以使用以下代码仅使用以下代码生成密集矩阵numpy：

row_index = np.array([0, 1, 2])
column_index = np.array([2, 1, 0])
values = np.array([4, 5, 6])

# numpy dense
M = np.zeros((np.max(row_index) + 1, np.max(column_index) + 1))
M[row_index, column_index] = values

在我的机器上，创建矩阵（最后两行）大约需要 6.3 μs 才能运行。我将它与以下使用的代码进行了比较scipy.sparse：

# scipy sparse
M = scipy.sparse.csc_matrix((values, (row_index, column_index)),
                            shape=(np.max(row_index) + 1, np.max(column_index) + 1))

运行大约需要 80 μs。因为你要的是创建稀疏数组的方法，所以我把第一个实现改成下面的代码，这样创建的ndarray就变成了稀疏数组：

# numpy sparse
M = np.zeros((np.max(row_index) + 1, np.max(column_index) + 1))
M[row_index, column_index] = values
M = scipy.sparse.csc_matrix(M)

运行大约需要 82 μs。这段代码的瓶颈显然是创建稀疏矩阵的操作。

请注意，该scipy.sparse方法可以很好地作为矩阵大小的函数进行缩放，最终成为较大矩阵的最快方法（在我的机器上，从大约 360×360 开始）。请参阅下图，了解每种方法的速度与矩阵大小的函数关系，从 10×10 矩阵到 1000×1000 矩阵。图中的一些异常值很可能是由于我机器上的其他程序干扰造成的。此外，我不确定在 ~360×360 和 ~510×510 的 numpy 密集方法中“跳跃”背后的技术细节。我还添加了用于运行此比较的代码，以便您可以在自己的机器上运行它。

import timeit
import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse

def generate_indices(num_values):
    row_index = np.arange(num_values)
    column_index = np.arange(num_values)[::-1]
    values = np.arange(num_values)
    
    return row_index, column_index, values

def numpy_dense(N, row_index, column_index, values):
    start = timeit.default_timer()
    for _ in range(N):
        M = np.zeros((np.max(row_index) + 1, np.max(column_index) + 1))
        M[row_index, column_index] = values
    end = timeit.default_timer()
    
    return (end - start) / N

def numpy_sparse(N, row_index, column_index, values):
    start = timeit.default_timer()
    for _ in range(N):
        M = np.zeros((np.max(row_index) + 1, np.max(column_index) + 1))
        M[row_index, column_index] = values
        M = scipy.sparse.csc_matrix(M)
    end = timeit.default_timer()
    
    return (end - start) / N

def scipy_sparse(N, row_index, column_index, values):
    start = timeit.default_timer()
    for _ in range(N):
        M = scipy.sparse.csc_matrix((values, (row_index, column_index)),
                                    shape=(np.max(row_index) + 1, np.max(column_index) + 1))
    end = timeit.default_timer()
    
    return (end - start) / N

ns = np.arange(10, 1001, 10)  # matrix size to try with

runtimes_numpy_dense, runtimes_numpy_sparse, runtimes_scipy_sparse = [], [], []

for n in ns:
    print(n)
    indices = generate_indices(n)
    # number of iterations for timing
    # ideally, you want this to be as high as possible,
    # but I didn't want to wait very long for this plot
    N = 1000 if n < 500 else 100
    runtimes_numpy_dense.append(numpy_dense(N, *indices))
    runtimes_numpy_sparse.append(numpy_sparse(N, *indices))
    runtimes_scipy_sparse.append(scipy_sparse(N, *indices))

fig, ax = plt.subplots()
ax.plot(ns, runtimes_numpy_dense, 'x-', markersize=4, label='numpy dense')
ax.plot(ns, runtimes_numpy_sparse, 'x-', markersize=4, label='numpy sparse')
ax.plot(ns, runtimes_scipy_sparse, 'x-', markersize=4, label='scipy sparse')

ax.set_yscale('log')
ax.set_xlabel('Matrix size')
ax.set_ylabel('Runtime (s)')
ax.legend()

plt.show()