[解決済み] Matlabで3角形の行列を作るループをベクトル化する

質問

3角形の行列を作るforループがあるのですが、こんな感じです。

m = 5;
Fo = 0.35;

A(1,1) = 1+2*Fo;
A(1,2) = 1+Fo;

for i = 2:m-1
    A(i,i-1) = 1;
    A(i,i) = 2;
    A(i,i+1) = 3;
end

A(m,m-1) = 4;
A(m,m) = 5;

を出力するところ。

A =
    1.7000    1.3500         0         0         0
    1.0000    2.0000    3.0000         0         0
         0    1.0000    2.0000    3.0000         0
         0         0    1.0000    2.0000    3.0000
         0         0         0    4.0000    5.0000

私は、for-loopを置き換えるために以下を使用して、3角形の行列の作成をベクトル化しようとしています。

i = 2:m-1;
A(i,i-1) = 1;
A(i,i) = 2;
A(i,i+1) = 3;

残念ながら、この出力は正しくありません。

A =
    1.7000    1.3500         0         0         0
    1.0000    2.0000    3.0000    3.0000    3.0000
    1.0000    2.0000    3.0000    3.0000    3.0000
    1.0000    2.0000    3.0000    3.0000    3.0000
         0         0         0    4.0000    5.0000

このような行列を、for-loopの代わりにベクトル化を使って作成することは可能でしょうか？私は最終的にもっと大きくて複雑な3角形の行列を作る必要があるので、ベクトル化を使って処理を高速化したいと思っています。

どのように解決するのですか？

ルイスさんのご指摘の通り、sub2indは動作します。diagも使えると思います。例えば、こんな感じです。

m = 5;
Fo = 0.35;

d = [1+2*Fo; repmat(2,m - 2,1); 5];
ud = [1+Fo;repmat(3,m-2,1)];
sd = [ones(m-2,1);4];

A = diag(d) + diag(ud,1) + diag(sd,-1)
A =
          1.7         1.35            0            0            0
            1            2            3            0            0
            0            1            2            3            0
            0            0            1            2            3
            0            0            0            4            5

問題は、この構成で多くの足し算を行っており、足し算のほとんどが0であることです。この解決策を嫌うもう一つの理由は、完全な行列を生成してしまうからです。一般的にはお勧めしませんが、視覚的に直感的でいい解決法です。(ところで、今気づいたのですが、diagのヘルプには、まさにこの解法で3角形の行列が生成されています)。

それよりも、疎な行列の使い方を学ぶ方がはるかに良い選択です。三角行列はスパース（sparse）です。その能力を利用しましょう。そのためには、spdiagsの使い方を学ぶか、最低でもsparseについて学ぶ必要があります。

Aをスパースな3角形行列として構築してみましょう。mの値を大きくして、節約できることがわかるようにします。

m = 500;
Fo = 0.35;

d = [1+2*Fo; repmat(2,m - 2,1); 5];
ud = [1+Fo;repmat(3,m-2,1)];
sd = [ones(m-2,1);4];

A = diag(d) + diag(ud,1) + diag(sd,-1);
As = spdiags([[sd;0],d,[0;ud]],[-1 0 1],m,m);

whos A*
  Name        Size               Bytes  Class     Attributes
  A         500x500            2000000  double              
  As        500x500              27976  double    sparse

つまり、スパースフォームは28k、フルバージョンは2メガバイトで保存できたのです。

これらの配列をスパースにした状態で使い始めると、本当の意味でメリットが出てきます。例えば、バックスラッシュを使う。

y = rand(m,1);

tic,x = A\y;toc
Elapsed time is 0.002847 seconds.

tic,xs = As\y;toc
Elapsed time is 0.000290 seconds.

私自身のコードでこの方法を示すことも必要でしょう。 ブラックトライディグ MATLAB Centralのファイル交換で見つけることができます。これはブロック3角形配列を生成するために設計されたものですが、スカラーブロックを持っているだけなので、この問題でも動作します。

Ab = blktridiag(reshape(d,1,1,m),reshape(sd,1,1,m-1),reshape(ud,1,1,m-1));
As - Ab
ans =
   All zero sparse: 500-by-500

最後に，natan が指摘するように，完全な行列が目的であれば，関数 full は，疎な行列の入力に対してそれを実現します．しかし、多くの場合、疎な形式は大きな利点となる。疎な行列の使い方を学んでください。そうすれば、きっと満足できるはずです。