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

Question

在此先感谢您的帮助。

使用 Python（主要是 numpy），我试图计算一个上三角矩阵，其中每一行 "j" 是几何级数的第一个 j 项，所有行都使用相同的参数.

例如，如果我的参数是 B（其中 abs(B)=<1，即 [-1,1] 中的 B），则第 1 行将是 [1 B B^2 B^3 ... B^(N-1)]，第 2 行将是 [0 1 B B^2...B^(N-2)] ... 第 N 行将是 [0 0 0 ... 1].

此计算是贝叶斯大都会-吉布斯采样器的关键，因此需要对 "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]

所以这只是逐行填充矩阵，但是使用了一个循环来减慢速度。

我想知道是否有人运行以前了解过这个，或者对如何更快地计算这个矩阵有更好的想法。

我以前从未在 Whosebug 上发过帖子，但在任何地方都没有看到这个问题，所以我想问一下。

让我知道是否有更好的地方可以问这个问题，以及我是否应该提供更多详细信息。

Answer 1

您可以使用 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.    ]])

如果你对数组的使用是严格的"read only"，你可以玩一些numpy strides的技巧来快速创建一个只使用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

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

Fastest way to compute upper-triangular matrix of geometric series (Python)

python

numpy

linear-algebra

vectorization