7

(2N)x(2N) numpy 配列の主対角線に沿ってあるサイズ 2x2 の対角ブロックを抽出するきちんとした方法を探しています (つまり、そのようなブロックが N 個存在します)。これは numpy.diag を一般化し、主対角線に沿って要素を返します。これは 1x1 ブロックと見なすことができます (もちろん、numpy はこのようには表現しません)。

これをもう少し広く表現すると、(MN)x(MN) 配列から N MxM ブロックを抽出したい場合があります。これは scipy.linalg.block_diag の補数のようです。ブロックをブロックをブロック対角行列に変換する方法 (NumPy)で適切に説明されており、block_diag が配置した場所からブロックを引き出します。

その質問の解決策からコードを借りて、これを設定する方法は次のとおりです。

>>> a1 = np.array([[1,1,1],[1,1,1],[1,1,1]])
>>> a2 = np.array([[2,2,2],[2,2,2],[2,2,2]])
>>> a3 = np.array([[3,3,3],[3,3,3],[3,3,3]])
>>> import scipy.linalg
>>> scipy.linalg.block_diag(a1, a2, a3)
array([[1, 1, 1, 0, 0, 0, 0, 0, 0],
       [1, 1, 1, 0, 0, 0, 0, 0, 0],
       [1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 2, 2, 2, 0, 0, 0],
       [0, 0, 0, 2, 2, 2, 0, 0, 0],
       [0, 0, 0, 2, 2, 2, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 3, 3, 3],
       [0, 0, 0, 0, 0, 0, 3, 3, 3],
       [0, 0, 0, 0, 0, 0, 3, 3, 3]])

次に、次のような関数が必要です

>>> A = scipy.linalg.block_diag(a1, a2, a3)
>>> extract_block_diag(A, M=3)
array([[[1, 1, 1],
        [1, 1, 1],
        [1, 1, 1]],
       [[2, 2, 2],
        [2, 2, 2],
        [2, 2, 2]],
       [[3, 3, 3],
        [3, 3, 3],
        [3, 3, 3]]])

numpy.diag との類推を続けるために、非対角ブロックを抽出することもできます: N - k 番目のブロック対角線上のそれらの k。(そして、ブロックを主対角線から離して配置できるようにする block_diag の拡張は確かに便利ですが、それはこの質問の範囲ではありません。) 上記の配列の場合、次の結果が得られます。

>>> extract_block_diag(A, M=3, k=1)
array([[[0, 0, 0],
        [0, 0, 0],
        [0, 0, 0]],
       [[0, 0, 0],
        [0, 0, 0],
        [0, 0, 0]]])

この質問でカバーされている stride_tricks の使用は、この機能に似たものを生成することを目的としていますが、ストライドがバイトレベルで動作することは理解していますが、これはあまり適切ではないようです。

動機として、これは、要素自体がスカラーではなく 2x2 行列である共分散行列 (つまり、分散) の対角要素を抽出したい状況から生じます。

編集:クリスからの提案に基づいて、次の試みを行いました:

def extract_block_diag(A,M,k=0):
    """Extracts blocks of size M from the kth diagonal
    of square matrix A, whose size must be a multiple of M."""

    # Check that the matrix can be block divided
    if A.shape[0] != A.shape[1] or A.shape[0] % M != 0:
        raise StandardError('Matrix must be square and a multiple of block size')

    # Assign indices for offset from main diagonal
    if abs(k) > M - 1:
        raise StandardError('kth diagonal does not exist in matrix')
    elif k > 0:
        ro = 0
        co = abs(k)*M 
    elif k < 0:
        ro = abs(k)*M
        co = 0
    else:
        ro = 0
        co = 0

    blocks = np.array([A[i+ro:i+ro+M,i+co:i+co+M] 
                       for i in range(0,len(A)-abs(k)*M,M)])
    return blocks

上記のデータに対して次の結果が返されます。

D = extract_block_diag(A,3)
[[[1 1 1]
  [1 1 1]
  [1 1 1]]

 [[2 2 2]
  [2 2 2]
  [2 2 2]]

 [[3 3 3]
  [3 3 3]
  [3 3 3]]]

D = extract_block_diag(A,3,-1)
[[[0 0 0]
  [0 0 0]
  [0 0 0]]

 [[0 0 0]
  [0 0 0]
  [0 0 0]]]
4

5 に答える 5

4

出発点として、次のようなものを使用できます

>>> a
array([[1, 1, 1, 0, 0, 0, 0, 0, 0],
   [1, 1, 1, 0, 0, 0, 0, 0, 0],
   [1, 1, 1, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 2, 2, 2, 0, 0, 0],
   [0, 0, 0, 2, 2, 2, 0, 0, 0],
   [0, 0, 0, 2, 2, 2, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 3, 3, 3],
   [0, 0, 0, 0, 0, 0, 3, 3, 3],
   [0, 0, 0, 0, 0, 0, 3, 3, 3]]) 

>>> M = 3
>>> [a[i*M:(i+1)*M,i*M:(i+1)*M] for i in range(a.shape[0]/M)]
[array([[1, 1, 1],
   [1, 1, 1],
   [1, 1, 1]]), array([[2, 2, 2],
   [2, 2, 2],
   [2, 2, 2]]), array([[3, 3, 3],
   [3, 3, 3],
   [3, 3, 3]])]
于 2012-05-31T12:12:28.083 に答える
4

ビューで行うこともできます。これはおそらくインデックス作成アプローチよりも高速です。

import numpy as np
import scipy.linalg

a1 = np.array([[1,1,1],[1,1,1],[1,1,1]])
a2 = np.array([[2,2,2],[2,2,2],[2,2,2]])
a3 = np.array([[3,3,3],[3,3,3],[3,3,3]])

b = scipy.linalg.block_diag(a1, a2, a3)
b[1,4] = 4
b[1,7] = 5
b[4,1] = 6
b[4,7] = 7
b[7,1] = 8
b[7,4] = 9
print b

def extract_block_diag(a, n, k=0):
    a = np.asarray(a)
    if a.ndim != 2:
        raise ValueError("Only 2-D arrays handled")
    if not (n > 0):
        raise ValueError("Must have n >= 0")

    if k > 0:
        a = a[:,n*k:] 
    else:
        a = a[-n*k:]

    n_blocks = min(a.shape[0]//n, a.shape[1]//n)

    new_shape = (n_blocks, n, n)
    new_strides = (n*a.strides[0] + n*a.strides[1],
                   a.strides[0], a.strides[1])

    return np.lib.stride_tricks.as_strided(a, new_shape, new_strides)

print "-- Diagonal blocks:"
c = extract_block_diag(b, 3)
print c

print "-- They're views!"
c[0,1,2] = 9
print b[1,2]

print "-- 1st superdiagonal"
c = extract_block_diag(b, 3, 1)
print c

print "-- 2nd superdiagonal"
c = extract_block_diag(b, 3, 2)
print c

print "-- 3rd superdiagonal (empty)"
c = extract_block_diag(b, 3, 3)
print c

print "-- 1st subdiagonal"
c = extract_block_diag(b, 3, -1)
print c

print "-- 2nd subdiagonal"
c = extract_block_diag(b, 3, -2)
print c
于 2012-06-02T13:23:04.000 に答える
0

https://stackoverflow.com/a/8070716/219229での unutbu の回答に基づいて、次のようになります:(ところで、バイトレベルでの操作の何が問題なのですか?)

from numpy.lib.stride_tricks import as_strided

...

def extract_block_diag(A, M=3, k=0):
   ny, nx = A.shape
   ndiags = min(map(lambda x: x//M, A.shape))
   offsets = (nx*M+M, nx, 1)
   strides = map(lambda x:x*A.itemsize, offsets)
   if k > 0:
       B = A[:,k*M]
       ndiags = min(nx//M-k, ny//M)
   else:
       k = -k
       B = A[k*M]
       ndiags = min(nx//M, ny//M-k)
   return as_strided(B, shape=(ndiags,M,M),
                     strides=((nx*M+M)*A.itemsize, nx*A.itemsize, A.itemsize))

使用例:

# a1, a2, a3 from your example
>>> bigA = scipy.linalg.block_diag(a1, a2, a3)
>>> extract_block_diag ( bigA, 3 )
array([[[1, 1, 1],
        [1, 1, 1],
        [1, 1, 1]],

       [[2, 2, 2],
        [2, 2, 2],
        [2, 2, 2]],

       [[3, 3, 3],
        [3, 3, 3],
        [3, 3, 3]]])
>>> extract_block_diag ( bigA, 3 )[2]
array([[3, 3, 3],
       [3, 3, 3],
       [3, 3, 3]])
>>> extract_block_diag(np.arange(1,9*9+1).reshape(9,9),3,1)
[[[ 4  5  6]
 [13 14 15]
 [22 23 24]]

[[34 35 36]
 [43 44 45]
 [52 53 54]]]

上記の関数はビューを返すことに注意してください。返された配列の配列内で何かを変更すると、オリジナルも影響を受けます。必要に応じてコピーを作成します。

于 2013-11-16T02:12:39.320 に答える
0

単純なアプローチを使用したくない特別な理由はありますか?

>>> A=np.array([[1,1,1,1],[2,2,2,2],[3,3,3,3],[4,4,4,4]])
>>> M1s,M2s=0,2 # start from A[M1s,M2s]
>>> M=2  # extract an MxM block
>>> for a in A[M1s:M1s+M]:
...   print a[M2s:M2s+M]
... 
[1 1]
[2 2]
>>>
于 2012-05-31T12:02:56.933 に答える