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

一个问题是我们在这个组合中做了很多加法，而且大部分加法都是零。不喜欢这个解决方案的另一个原因是它会生成一个完整的矩阵。我一般不建议这样做，但它是一个很好的、视觉上直观的解决方案。（顺便说一句，我刚刚注意到 diag 的帮助准确地显示了这个创建三对角矩阵的解决方案。）

更好的选择是学习使用稀疏矩阵。三对角矩阵是稀疏的。使用该功能。为此，请学习如何使用 spdiags，或者至少了解 sparse。

让我们将 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.

我想我还应该展示如何使用我自己的代码blktridiag来执行此操作，该代码位于 MATLAB Central 文件交换中。它实际上是为生成块三对角数组而设计的，但它也可以解决这个问题，因为我们只有标量块。

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;

您替换循环的方法的问题在于，这两个索引被认为是定义一个块（两个索引的所有组合），而不是对角线（第一个索引的每个值与第二个索引中的每个对应值）。

于 2013-09-05T23:20:42.833 回答