各対角番号に対する対角和からなる配列

写真のような正方行列Aを持っている際に、自身より右下の第k対角成分の和をとった行列Bを作成したいのですが、何か簡単に書く方法はないでしょうか。

現在は縮小行列を作り、各縮小行列の第j番目の対角成分を取得する形で計算していますが、かなり遅いし効率も悪いなと感じています。

A = reshape(1:9,3,3);
B = zeros(3,3);

for i = 1:3
    A_partial = A(i:end,i:end);
  
    for j = 1:3-i+1
        B(i,i+j-1) = sum(diag(A_partial,j-1));
        B(i+j-1,i) = sum(diag(A_partial,-1*(j-1)));
    end

end
A
B

[Atsushi Ueno]

メモ：同じ主旨の問題がCodyに登録されました。逆の対角方向を選択するオプションも求められます。これを解くと他のCody参加者の回答を閲覧する事が出来ます。突飛なアイデアを持つ回答者が沢山いますよ。

Sub-Diagonal Sum - MATLAB Cody - MATLAB Central (mathworks.com)

[Yamada Taro]

@Atsushi Uenoさん、有難うございます！！こんなサイトがあるんですね、見てみます。

Answer 1

[NO MIYA]

もっと簡単な方法があると思いますが、とりあえずこんなのでどうでしょうか。

A⇒AA、B⇒BBと考えてください。

一応n×n（3<n）の正方行列に対応できるようにもしておきました。

padarrayを使うにはImage Processing Toolboxが必要なのですが、やっていることは配列の0埋めなので自作関数でいくらでもやりようがあると思います（結果、計算量が増えるかもしれませんが）

for文を回す回数が投稿者様のが二重ループで3(n)回以上なのに対し、下記プログラムは2(n-1)回で面倒な対角要素を省いたので最低限の改良はできたかと思います。

AA = reshape(1:9,3,3); % 基の正方行列
BB = AA;               % 出力したい配列
n  = size(AA,1);       % 正方行列の1辺の長さ

%%% 処理
for i = 2:1:n
    AA_plus = padarray(AA(i:end,i:end),[i-1,i-1],0,'post');
    BB = BB + AA_plus;
end
BB

BB = 3×3

15 12 7 8 14 8 3 6 9

[NO MIYA]

追記です。

コメントが不親切なプログラムだけ貼ってしまったのでアルゴリズムを説明させていただきます。

字が汚くて申し訳ないのですが、Bは下図のように分解できるので（丁度投稿者様がプログラム中の 'A_partial' でやっていらっしゃる感じです）

1. 'A_partial'に当たる正方行列を抽出（Bの事前割り当て初期化をzerosではなく、B = A; にすることで繰り返しを一回減らせます）
2. 0で埋めてサイズをn×nに合わせる
3. 足す

[Atsushi Ueno]

同じ考えです。

• 任意サイズの2次元行列で目的の出力Bが得られます
• ゼロ埋めで行列サイズを合わせなくても同じ事が出来ます

A = ones(3,4);
B = A;
for i = 1:min(size(A))-1
    B(1:end-i,1:end-i) = B(1:end-i,1:end-i) + A(1+i:end,1+i:end);
end
B

B = 3×4

3 3 2 1 2 2 2 1 1 1 1 1

[Yamada Taro]

@NO MIYAさん、@Atsushi Uenoさん、ご丁寧にありがとうございます。

おかげさまで無事理解でき、実装もできました。 一点疑問なのですが、padarrayを用いる理由としてはスライスに代入する処理があまり早くないからしょうか？

とりあえず現状ではpadarrayをいれていないので、スライスにスライスを代入する形にしています。

[NO MIYA]

padarrayを使わずとも、

@Atsushi Uenoさんのようにインデックス指定すれば良い話でした。

勉強になりました。

[Atsushi Ueno]

@Yamada Taroさんの疑問に対しては「両者とも処理速度の面は気にせず回答した」事になります。

自分は処理速度ではなく単にシンプルな記述方法の提案でした(計算量も減らせますが)。スライスの代入よりメモリ上で連続した行列全体をコピーする方が速いと思いますが、0埋めにも計算量が必要なのでどっこいどっこいといったところでしょうか。(めんどいから検証する気が無い)

「行数分のストライドを持つ線形インデックスでアクセス」すれば処理速度が上がるのでは？と思いました(BLASの実装を参照)。しかし、下記の理由で処理速度が上がる事はないだろうと思いました。

• 「自身より右下の第k対角成分」を判定しないと折り返し反対側まで加算してしまう
• これはfortranやCで実装する場合の話で、MATLABなら行列操作に任せるべき
• MATLABの行列操作もBLAS相当かそれ以上に高効率な実装のはずです(読めない)
• MATLABはまた静的型付けでないので、行列が連続したメモリ配置かどうか不明

そもそもMATLABはインタプリタなので、速度重視ならコンパイルすべきです。

A = ones(3,4);
h = size(A,1) + 1;
B = A;
for i = 1:numel(A)
    B(i) = sum(A(i:h:end)); % 添字の大小が逆転する項は加算してはならないが実装していない
                            %(実装すると結局結構な計算量になるので処理速度は期待できない)
end
B % これでは左下の三角形が期待通りに計算されない

B = 3×4

3 3 2 1 3 2 2 1 3 2 1 1

[Yamada Taro]

@Atsushi Uenoさん

ご返信ありがとうございます、こちら返事が遅くなってしまい大変申し訳ございません。

fortranは触ったことはないのですが、matlabであればできる限り既にかなり高速化がされている行列演算にできる限り頼る形にした方が良いとのこと、勉強になりました。

また速度重視ならコンパイルすべき、とのご指摘もありがとうございます。こちら試したことがないので、高速化が必要になった際に挑戦してみようと思います。ありがとうございました。

Answer 2

[Hernia Baby]

せっかくなので@Yamada Taroさんのイメージどおりに動くよう、対角成分をうまく使う関数を作成しました

clc,clear;
A = reshape(1:9,3,3);
B = zeros(size(A));

キモとなる関数です

1. diagを使って対角成分のみをとる
2. flipして反転
3. cumsumで累積和
4. flipして順序を戻す
5. diagで対角行列をつくる

f = @(X,x) diag(flip(cumsum(flip(diag(X,x)))),x);

上記の関数を - (行番号 - 1) ~ 列番号 - 1 分行います

-1としているのは角っこは計算する必要がないからです。

時間も計測してみましょうかね

tic
for ii = -height(A)+1:width(A)-1
    B = B + f(A,ii);
end
toc

Elapsed time is 0.007458 seconds.

確認してみましょう

B

B = 3×3

15 12 7 8 14 8 3 6 9

[Yamada Taro]

有難うございます。おかげさまで無事実装できました。