matlab - 将任意数量的矩阵乘以任意次数

Question

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

我正在编写一个通用例程来计算任意层数的薄膜反射与光频率的关系。对于每个光频率 W，每一层都有一个折射率 N 和一个关联的 2x2 传输矩阵 L 和 2x2 接口(interface)矩阵 I 这取决于折射率和层的厚度。如果 n 是层数，而 m 是频率数，那么我可以将索引矢量化为 n x m 矩阵，但是为了计算反射每个频率，我都必须做嵌套循环。由于我最终将其用作试衣程序的一部分，因此我将不胜感激我能做些什么来加快速度。

这应该提供一个最低限度的工作示例:

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，将它们与前一个循环的结果相乘，然后继续直到我到达堆栈的底部。然后我获取第一个和第三个值，取比率，这就是该频率下的反射率。然后我转到下一个频率并重新做一遍。

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

Answer 1

不在 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 block 矩阵，其形式如下:

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

e- 为 exp(-1i*e)，e+ 为 exp(1i*e)'

请参阅 kron 了解如何获得 block 矩阵形式，以矢量化传播 C=C*I*L 只需取 M^n