计算几何级数的上三角矩阵最快方法（Python）

Question

计算几何级数的上三角矩阵最快方法（Python）

3

谢谢您的帮助。使用Python（主要是numpy），我正在尝试计算一个上三角矩阵，其中每一行“j”都是几何级数的前j项，所有行都使用相同的参数。

例如，如果我的参数是B（其中abs（B）=＜1，即B在[-1,1]范围内），那么第一行将是[1 B B ^ 2 B ^ 3 ... B ^（N-1）]，第二行将是[0 1 B B ^ 2 ... B ^（N-2）] ...第N行将是[0 0 0 ... 1]。

这个计算对于贝叶斯Metropolis-Gibbs采样器非常关键，因此需要为“B”的新值重复执行数千次。

目前我已经尝试了两种方法：

方法1 - 大多数向量化：

B_Matrix = np.triu(np.dot(np.reshape(B**(-1*np.array(range(N))),(N,1)),np.reshape(B**(np.array(range(N))),(1,N))))

本质上, 这是一个Nx1矩阵和1xN矩阵相乘的上三角部分:

上三角 ([1 B^(-1) B^(-2) ... B^(-(N-1))]' * [1 B B^2 B^3 ... B^(N-1)])

这个方法在小N时很好用 (代数上是正确的), 但是对于大的N它会出现错误。而且当B=0时也会出现错误 (但实际上应该是允许的)。我相信这是因为对于小的B和大的N，B^(-N) ~ inf导致的。

方法二:

B_Matrix = np.zeros((N,N))
B_Row_1 = B**(np.array(range(N)))
for n in range(N):
    B_Matrix[n,n:] = B_Row_1[0:N-n]

这样的方法按行填充矩阵，但使用了循环，这会减慢计算速度。

我想知道是否有人遇到过这种情况，或者有更好的方法来更快地计算这个矩阵。

我以前从未在stackoverflow上发过帖子，但没有看到这个问题，所以想问一下。

如果有更好的地方可以询问此事，请告诉我，并且如果需要提供更多细节，请告诉我。

- Alexander Monte Calvo

我认为你在这里首先要问自己的问题是：我需要多快？如果你想要带有循环的版本（最终可能更快），你可能想要检查一下完全静态类型的 cython 函数能够达到多快。如果还不够快，你仍然可以选择完全用 C 编写此函数，并将其集成到你的 Python 代码中。 - cel

修改了行，谢谢！我想我需要去了解一下Cython，虽然我没有经验，但我听说它可以让程序运行更快。 - Alexander Monte Calvo

这个数组将是“只读”的吗，还是你会在原地进行更改？ - Warren Weckesser

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Warren Weckesser · Accepted Answer

你可以使用 scipy.linalg.toeplitz：

In [12]: n = 5

In [13]: b = 0.5

In [14]: toeplitz(b**np.arange(n), np.zeros(n)).T
Out[14]: 
array([[ 1.    ,  0.5   ,  0.25  ,  0.125 ,  0.0625],
       [ 0.    ,  1.    ,  0.5   ,  0.25  ,  0.125 ],
       [ 0.    ,  0.    ,  1.    ,  0.5   ,  0.25  ],
       [ 0.    ,  0.    ,  0.    ,  1.    ,  0.5   ],
       [ 0.    ,  0.    ,  0.    ,  0.    ,  1.    ]])

如果你只是“只读”地使用数组，你可以通过 numpy 的步幅技巧来快速创建一个仅使用 2*n-1 个元素（而不是 n^2）的数组：

In [55]: from numpy.lib.stride_tricks import as_strided

In [56]: def make_array(b, n):
   ....:     vals = np.zeros(2*n - 1)
   ....:     vals[n-1:] = b**np.arange(n)
   ....:     a = as_strided(vals[n-1:], shape=(n, n), strides=(-vals.strides[0], vals.strides[0]))
   ....:     return a
   ....: 

In [57]: make_array(0.5, 4)
Out[57]: 
array([[ 1.   ,  0.5  ,  0.25 ,  0.125],
       [ 0.   ,  1.   ,  0.5  ,  0.25 ],
       [ 0.   ,  0.   ,  1.   ,  0.5  ],
       [ 0.   ,  0.   ,  0.   ,  1.   ]])

如果您要就地修改数组，请复制make_array(b, n)返回的结果。也就是说，arr = make_array(b, n).copy()。

make_array2函数包含了评论中@Jaime提出的建议：

In [30]: def make_array2(b, n):
   ....:     vals = np.zeros(2*n-1)
   ....:     vals[n-1] = 1
   ....:     vals[n:] = b
   ....:     np.cumproduct(vals[n:], out=vals[n:])
   ....:     a = as_strided(vals[n-1:], shape=(n, n), strides=(-vals.strides[0], vals.strides[0]))
   ....:     return a
   ....: 

In [31]: make_array2(0.5, 4)
Out[31]: 
array([[ 1.   ,  0.5  ,  0.25 ,  0.125],
       [ 0.   ,  1.   ,  0.5  ,  0.25 ],
       [ 0.   ,  0.   ,  1.   ,  0.5  ],
       [ 0.   ,  0.   ,  0.   ,  1.   ]])

make_array2的速度比make_array快两倍以上：

In [35]: %timeit make_array(0.99, 600)
10000 loops, best of 3: 23.4 µs per loop

In [36]: %timeit make_array2(0.99, 600)
100000 loops, best of 3: 10.7 µs per loop