matlab: ブロック対角成分の抽出

matlab で，ある行列からブロック対角成分を抽出したい場合， 効率的な実装はあるだろうか．

少し調べてみよう．

対角行列，とくにブロック対角行列の効率的作成法について，以前の記事で述べた．
今回は，逆に，すでに与えられた行列から，ブロック対角成分を抜き出すことを考える．

ブロック対角行列とは，厳密には，正方行列であって，その対角成分も正方行列であるものをさす（Wolfram MathWorld の Block Diagonal Matrix）．ここではもう少し緩やかに，正方とは限らないブロック行列を行・列の重複なく左上から右下に並べ，それ以外の成分が0であるような行列全般を扱うこととする．そのため，本記事でいうブロック対角成分は，必ずしも正方とは限らない．というか，対角成分を含まないブロックがあったりしてもよい．

ブロックサイズが揃っている（同一の）場合

以前の記事が応用できる．成分の抽出のために，マスク行列を作成するのだ．

まず，抽出元となる行列を A とする．抽出するのは [ m by n ] のブロック k 個だとすると，A は [ m*k by n*k ] の行列のはずである．

抽出元の行列(A)が手元になければ，試しにを作ってみよう．

m = 2;
n = 3;
k = 2;
A = rand(m*k, n*k);

さて，ここからブロック対角成分を抽出するためのマスクは，以下のように作る．ここが以前の記事の応用である．

C1 = sparse(true(m, n*k));
C1(m*k,n*k) = 0;
C2 = reshape(C1,[m*k*n,k]);
C2(m*(k*n+1),k) = 0;
C = reshape(C2, [m*k,n*k+1]);
C(:,end) = [];

これでマスクができた．つまり，この logical 行列 C を logical indexing (前に記事を書いた)に使うわけだ．C の様子は，
spy(C)
とすると視覚的に確認できる．こんな感じ：

spy(C) で表示されるフィギュア．所望のブロック対角の箇所(2x3のブロック2つ)が非零になっていることがわかる．非零要素の個数は nz=12.

あとは，成分を抽出すればいい．単にリニアインデクス(これも前に記事を書いた)順にベクトル d として抽出してよいなら

d = A(C);

でOK. ブロック毎に分けて格納したいなら，この下の「ブロックサイズが揃っていない（同一でない）場合」を参照．もとの行列と同形に再構築したい，つまり「元の行列からブロック対角成分以外を 0 にした行列(R)を作りたい」なら

R = sparse(m*k, n*k);
R(C) = A(C);

とすればよい．中身を確認すると，

>> A
A =
    0.9029    0.5036    0.5244    0.1588    0.6901    0.4772
    0.9123    0.6278    0.1008    0.9195    0.3914    0.9934
    0.8822    0.4148    0.5407    0.8950    0.4128    0.9987
    0.7367    0.4762    0.4750    0.4849    0.2321    0.0981

>> full(R)
ans =
    0.9029    0.5036    0.5244         0         0         0
    0.9123    0.6278    0.1008         0         0         0
         0         0         0    0.8950    0.4128    0.9987
         0         0         0    0.4849    0.2321    0.0981

となっており，ばっちり所望の動作を実現できている．

ブロックサイズが揃っていない（同一でない）場合

ブロックサイズ k11 x k12, k21 x k22, k31 x k32, ..., kn1 x kn2 のブロック成分を抽出する．
ここは素直にループを回してもいいだろう．

% ブロック対角成分のサイズを K に格納してあるとしよう．
K = [ ...
    k11, k12; ...
    k21, k22; ...
    k31, k32; ...
    ...
    kn1, kn2; ...
    ];

% 例えばこんな感じ．抽出元の行列を A とする．
A=rand(6,7);
K = [ 3,2; 1,2; 2,3];

% ブロック対角の個数は n 個．
n = size(K, 1);

% ブロック対角成分をセルに格納していく．
D = cell(n, 1);
kz1 = 0;
kz2 = 0;
for m=1:n
    kz1q = kz1 + K(m, 1);
    kz2q = kz2 + K(m, 2);
    D{m} = A( kz1+1:kz1q, kz2+1:kz2q );
    kz1 = kz1q;
    kz2 = kz2q;
end

% これで抽出は完了．

% ブロック対角行列を再構成してみよう．
% 入れ物となる sparse 行列を作成しておく．非零要素数 sum(prod(K,2)) を先に割り当て．
R = sparse([], [], [], sum(K(:, 1)), sum(K(:, 2)), sum(prod(K,2)));
kz1 = 0;
kz2 = 0;
for m=1:n
    kz1q = kz1 + K(m, 1);
    kz2q = kz2 + K(m, 2);
    R( kz1+1:kz1q, kz2+1:kz2q ) = D{m};
    kz1 = kz1q;
    kz2 = kz2q;
end

ブロックサイズが揃っていない場合にも，ループを回さず超高速で実行する方法はあるかもしれない．ご存知の方は情報頂けると嬉しいです．