将矩阵转换为 LAPACK 的最快实现

Question

我有一个函数可以将矩阵转换为 SciPy 中 LAPACK 例程所需的函数。我编写了简短的代码来说明我的工作：

import numpy as np
import time

N = 3

def convert_matrix_lapack(Amat,out):
    for ii in range(N):
        out[(2*N-2-ii):(3*N-2-ii), ii] = Amat[:, ii]

A = np.arange(N**2).reshape(N, N)
A_lapack = np.zeros((3*N-2,N), dtype=np.float)

tt = time.time()
convert_matrix_lapack(A, A_lapack)
print(time.time() - tt)

在实践中通用矩阵（N=3），

>>> A
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

变成

>>> A_lapack
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 2.],
       [0., 1., 5.],
       [0., 4., 8.],
       [3., 7., 0.],
       [6., 0., 0.]])

如何使用 NumPy 中的内置函数来固定任意 N 的代码（我的目标是 N 小于 50）？

Answer 1

这应该可行：

from scipy.sparse import diags
import numpy as np

N = 3
A = np.arange(N**2).reshape(N,N)


offsets = np.arange(-N+1, 1)
A_lapack = np.flipud(diags(A, offsets, shape=(3*N-2, N)).toarray())

它是如何工作的：

1. 它使用scipy.sparse.diags 创建一个稀疏对角矩阵
2. 使用numpy.flipud 翻转它

如果我们包括使用np.zeros 的数组分配步骤，速度是相当的，但对于N=100 的数组来说更慢：

• 0.00021839141845703125秒（原码）
• 0.0017552375793457031 秒（我的解决方案）