如何在NumPy中不使用循环为大矩阵分配正方形子矩阵

Question

我应该如何使用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

- Joseph Hastings

Answer 1

您的测试不正确：对于i，j=0,0，您的M2 [] =分配只是覆盖了同一矩阵块。

当使用M1时获得对称矩阵的原因是因为您在单个循环中分配了M1值。

如果将循环拆分为两个：

for x in range(n):
      M1[i:i+n,j+x] = m[x]
for x in range(n): 
      M1[j+x,i:i+n] = m[x]

M1显然与M2相同。

综上所述，以下代码可行（相当于您的M2计算），但是！仅当对角线上方和下方的块之间没有重叠时才有效。 如果存在重叠，则必须决定如何处理。

xs=np.arange(4).reshape(2,2)
ys=np.zeros((7,7))
ys[i:i+n,j:j+n]=xs
ys[j:j+n,i:i+n]=xs.T
print ys
>> array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  2.,  3.,  0.,  0.],
       [ 0.,  0.,  2.,  0.,  0.,  0.,  0.],
       [ 0.,  1.,  3.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]])

- sega_sai

---

啊，说得好，代码写成这样很难进行向量化，因为循环中的两个赋值可能会相互影响。我以为我碰到了一个微妙的索引特性，但实际上解释要简单得多。谢谢！ - Joseph Hastings