1

次のような三重対角行列を作成する 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 ループを置き換えて、三重対角行列の作成をベクトル化しようとしています。

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 ループの代わりにベクトル化を使用してそのような行列を作成することは可能ですか? 最終的には、はるかに大きく複雑な三重対角行列を作成する必要があるため、ベクトル化を使用してプロセスを高速化したいと考えています。

4

2 に答える 2

2

Luis が指摘するように、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

問題は、この構成で多くの加算を行ったことであり、ほとんどの加算はゼロです。このソリューションを好まないもう 1 つの理由は、完全な行列を生成することです。一般的にはお勧めしませんが、視覚的に直感的な優れたソリューションです。(ちなみに、diag のヘルプが、三重対角行列を作成するためのこのソリューションを正確に示していることに気付きました。)

はるかに良い選択は、疎行列の使い方を学ぶことです。三重対角行列は SPARSE です。その機能を使用します。そのためには、spdiags の使用方法を学習するか、少なくともスパースについて学習してください。

A をスパース三重対角行列として構築しましょう。節約を確認できるように、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 ファイル交換で見つけた独自のコードblktridiagを使用してこれを行う方法も示す必要があると思います。実際にはブロックの三重対角配列を生成するように設計されていますが、スカラー ブロックしかないので、この問題でも機能します。

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 がそれを実行します。ただし、多くの場合、スパース形式は大きな利点になります。疎行列の使い方を学びます。あなたはあなたがしたことを幸せに思うでしょう。

于 2013-09-06T00:33:29.513 に答える
1

for以下を使用してループをベクトル化できますsub2ind。

m = 5;
Fo = 0.35;

A = zeros(m); % initialize
A(sub2ind([m m],2:m, 1:m-1)) = 1;
A(sub2ind([m m],1:m, 1:m)) = 2;
A(sub2ind([m m],1:m-1, 2:m)) = 3;

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

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

ループを置き換えるアプローチの問題は、2 つのインデックスがブロック (2 つのインデックスのすべての組み合わせ) を定義すると見なされ、対角線 (最初のインデックスの各値と 2 番目の各対応する値) ではないと見なされることです。

于 2013-09-05T23:20:42.833 に答える