如何分配在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 

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.]])