将任意数目的矩阵乘以任意次数

Question

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

我正在编写一个计算薄膜反射的普通例程，它由任意层数与光学频率成正比。对于每个光学频率W，每个层都有一个折射指数N和一个相关的2x2传输矩阵L和2x2接口矩阵I，这取决于折射率和该层的厚度。如果n是层数，m是频率数，那么我可以将索引向量化为n矩阵，但是为了计算每个频率上的反射，我必须执行嵌套循环。因为我最终是使用它作为一个适当的例行程序的一部分，我能做的任何事情来加快它将是非常感谢的。

这应提供一个最低限度的工作实例：

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对于每个层和每个频率是不同的，因此它们在每个循环中都会发生变化。简单地取指数是行不通的。举个现实世界的例子，以空气中的肥皂泡沫为例。有三个层(air/soap/air)和两个接口。对于给定的频率，完整的传输矩阵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，将它们与前一个循环的结果相乘，然后一直进行到堆栈的底部。然后我得到第一个和第三个值，取比率，也就是那个频率上的反射率。然后我转到下一个频率，然后再做一次。

我怀疑答案将涉及到每一层的块对角矩阵，如下所述。

Answer 1

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

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秒。

Answer 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-作为exp(-1i*e)，e+作为exp(1i*e)‘

请参阅kron关于如何获取块矩阵形式，将传播C=C*I*L矢量化只需取M^n