我如何向量化这个循环,它使用 numpy 数组填充更大矩阵的两个方形子矩阵(也保持更大的矩阵对称):
for x in range(n):
assert m[x].shape == (n,)
M[i:i+n,j+x] = m[x]
M[j+x,i:i+n] = m[x]
这很诱人,但不同意上面的循环(参见下面的示例分歧):
assert m.shape == (n,n)
M[i:i+n,j:j+n] = m
M[j:j+n,i:i+n] = m
这是一个小例子(n>1 时崩溃):
from numpy import arange,empty,NAN
from numpy.testing import assert_almost_equal
for n in (1,2,3,4):
# make the submatrix
m = (10 * arange(1, 1 + n * n)).reshape(n, n)
N = n # example below, submatrix is the whole thing
# M1 using loops, M2 "vectorized"
M1 = empty((N, N))
M2 = empty((N, N))
M1.fill(NAN)
M2.fill(NAN)
i,j = 0,0 # not really used when (n == N)
# this results in symmetric matrix
for x in range(n):
assert m[x].shape == (n,)
M1[i:i+n,j+x] = m[x]
M1[j+x,i:i+n] = m[x]
# this does not work as expected
M2[i:i+n,j:j+n] = m
M2[j:j+n,i:i+n] = m
assert_almost_equal(M1,M1.transpose(),err_msg="M not symmetric?")
print "M1\n",M1,"\nM2",M2
assert_almost_equal(M1,M2,err_msg="M1 (loop) disagrees with M2 (vectorized)")
我们最终得到:
M1 = [10 30
30 40] # symmetric
M2 = [10 20
30 40] # i.e. m