3

我找到了几个问题/答案，用于向量化和加速在单个循环中将矩阵和向量相乘的例程，但我正在尝试做一些更通用的事情，即将任意数量的矩阵相乘，然后执行操作任意次数。

我正在编写一个通用例程，用于计算任意层数与光学频率的薄膜反射。对于每个光学频率W，每一层都有一个折射率N和一个相关的 2x2 传递矩阵L和 2x2 界面矩阵I，这取决于折射率和层的厚度。如果n是层数，并且m是频率数，那么我可以将索引向量化为一个 nxm 矩阵，但是为了计算每个频率的反射，我必须进行嵌套循环。由于我最终将其用作拟合例程的一部分，因此我将不胜感激我能做的任何事情来加快它的速度。

这应该提供一个最小的工作示例：

W = 1260:0.1:1400; %frequency in cm^-1
N = rand(4,numel(W))+1i*rand(4,numel(W)); %dummy complex index of refraction
D = [0 0.1 0.2 0]/1e4; %thicknesses in cm
[n,m] = size(N);
r = zeros(size(W));

for x = 1:m %loop over frequencies
  C = eye(2); % first medium is air

  for y = 2:n %loop over layers

    na = N(y-1,x);
    nb = N(y,x);

    %I = InterfaceMatrix(na,nb); % calculate the 2x2 interface matrix
    I = [1 na*nb;na*nb 1]; % dummy matrix

    %L = TransferMatrix(nb) % calculate the 2x2 transfer matrix
    L = [exp(-1i*nb*W(x)*D(y)) 0; 0 exp(+1i*nb*W(x)*D(y))]; % dummy matrix

    C = C*I*L;
  end

  a = C(1,1);
  c = C(2,1);

  r(x) = c/a; % reflectivity, the answer I want.
end

对于具有 2562 个频率的三层（空气/材料/基板）问题的两个不同极化运行两次，需要 0.952 秒，而使用三层系统的显式公式（矢量化）解决完全相同的问题需要 0.0265 秒。问题是，超过 3 层，显式公式很快变得难以处理，我必须为每个层数有一个不同的子例程，而以上是完全通用的。

是否有希望将此代码矢量化或以其他方式加速它？

（编辑补充说我在代码中留下了一些东西来缩短它，所以请不要尝试使用它来实际计算反射率）

编辑：为了澄清，每个层I和L每个频率都不同，所以它们在每个循环中都会改变。简单地取指数是行不通的。对于一个真实世界的例子，以空气中的肥皂泡为例。共有三层（空气/肥皂/空气）和两个接口。对于给定的频率，完整的传输矩阵C为：

C = L_air * I_air2soap * L_soap * I_soap2air * L_air;

和I_air2soap ~= I_soap2air。因此，我从L_air = eye(2)连续的层开始然后向下，计算 I_(y-1,y) 和 L_y，将它们与前一个循环的结果相乘，一直持续到我到达堆栈的底部。然后我抓住第一个和第三个值，取比率，这就是那个频率的反射率。然后我进入下一个频率并重新做一遍。

我怀疑答案将以某种方式涉及每一层的块对角矩阵，如下所述。

4

2 回答 2

2

不是在 matlab 旁边，所以这只是一个启动器，而不是双循环，您可以写na*nb为Nab=N(1:end-1,:).*N(2:end,:); 指数中的项nb*W(x)*D(y)可以写为e=N(2:end,:)*W'*D; 结果I*L是一个 2x2 块矩阵，具有以下形式：

M = [1, Nab; Nab, 1]*[e-, 0;0, e+] = [e- , Nab*e+ ; Nab*e- , e+]

使用e-as exp(-1i*e) 和e+as exp(1i*e)'

查看kron如何获得块矩阵形式，以矢量化传播C=C*I*L只需M^n

于 2013-03-20T20:55:55.967 回答
2

@Lama 通过建议块矩阵让我走上了正确的道路，但最终的答案最终变得更加复杂，所以我把它放在这里是为了后代。由于每一层的传输和接口矩阵都不同，我在层上的循环中保留，但构造一个大的稀疏块矩阵，其中每个块代表一个频率。

W = 1260:0.1:1400; %frequency in cm^-1
N = rand(4,numel(W))+1i*rand(4,numel(W)); %dummy complex index of refraction
D = [0 0.1 0.2 0]/1e4; %thicknesses in cm
[n,m] = size(N);
r = zeros(size(W));

C = speye(2*m); % first medium is air
even = 2:2:2*m;
odd = 1:2:2*m-1;

for y = 2:n %loop over layers

  na = N(y-1,:);
  nb = N(y,:);

  % get the reflection and transmission coefficients from subroutines as a vector
  % of length m, one value for each frequency
  %t = Tab(na, nb);
  %r = Rab(na, nb);
  t = rand(size(W)); % dummy vector for MWE
  r = rand(size(W)); % dummy vector for MWE

  % create diagonal and off-diagonal elements.  each block is [1 r;r 1]/t
  Id(even) = 1./t;
  Id(odd) = Id(even);
  Io(even) = 0;
  Io(odd) = r./t;

  It = [Io;Id/2].';
  I = spdiags(It,[-1 0],2*m,2*m);

  I = I + I.';

  b = 1i.*(2*pi*D(n).*nb).*W;
  B(even) = -b;
  B(odd)  =  b;
  L = spdiags(exp(B).',0,2*m,2*m);

  C = C*I*L;
end

a = spdiags(C,0);
a = a(odd).';

c = spdiags(C,-1);
c = c(odd).';

r = c./a; % reflectivity, the answer I want.

使用上面提到的 3 层系统，它并没有显式公式那么快，但它很接近，并且在一些分析后可能会变得更快一点。原始代码的完整版本为 0.97 秒，公式为 0.012 秒，稀疏对角线版本为 0.065 秒。

于 2013-03-21T01:16:10.357 回答