単一のループで行列とベクトルを乗算するルーチンをベクトル化して高速化するためのいくつかの質問/回答を見つけましたが、もう少し一般的なことをしようとしています。つまり、任意の数の行列を乗算してから実行します。任意の回数操作します。
任意の層数と光周波数から薄膜反射を計算するための一般的なルーチンを作成しています。各光周波数Wに対して、各層には屈折率と、屈折率と層の厚さに依存Nする関連する 2x2 転送マトリックスLと 2x2 インターフェイスマトリックスがあります。が層の数で、 が周波数の数であるI場合、インデックスを nxm 行列にベクトル化できますが、各周波数での反射を計算するには、ネストされたループを実行する必要があります。私は最終的にこれをフィッティングルーチンの一部として使用しているので、スピードアップするためにできることは何でも大歓迎です.nm
これは最小限の作業例を提供する必要があります。
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 周波数の 3 層 (空気/物質/基板) 問題の 2 つの異なる偏波に対してこれを 2 回実行すると、0.952 秒かかりますが、3 層システムの明示的な式 (ベクトル化) を使用してまったく同じ問題を解くには 0.0265 秒かかります。問題は、3 層を超えると、明示的な式が急速に扱いにくくなり、上記は完全に一般的ですが、層の数ごとに異なるサブルーチンが必要になることです。
このコードをベクトル化するか、それ以外の方法で高速化する希望はありますか?
（コードを短くするためにいくつかのことをコードから除外したことを編集して追加したため、実際に反射率を計算するためにこれを使用しないでください）
編集：明確にするために、IとLはレイヤーごと、周波数ごとに異なるため、ループごとに変化します。単に指数を取るだけではうまくいきません。現実世界の例として、空気中のシャボン玉の最も単純なケースを取り上げます。3 つの層 (空気/石鹸/空気) と 2 つのインターフェイスがあります。特定の周波数の場合、完全な伝達行列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 を計算し、前のループの結果と乗算し、スタックの一番下に到達するまで続けます。次に、1 番目と 3 番目の値を取得し、比率を取得します。これがその周波数での反射率です。それから次の周波数に移り、すべてをやり直します。
以下で説明するように、答えには、各レイヤーのブロック対角行列が何らかの形で含まれると思われます。