arrayfunの使い方について

現在、Ｍａｔｌａｂ2018ｂを使用し、アプリを作成しているのですが、処理を高速化させたい部分があります。

それはforループを使用して約25万回ループしている箇所です。

arrayfun関数を使用するとforループを使用しないで書くことができるようなのですが下記のような処理の場合でも

arrayfun関数でforループを使用せず書くことができるものでしょうか。

改善したい箇所は大体、下記のような感じになっており、最終的に25万行のデータＳを作成したいと思っています。

S=zeros(250000, 1, 2);

for i=1:250000

Xp = [F1(:, i), F2(:, i)];

S(i, :, :) = Fn(:, i) * Xp' * Xp * Xp' * C;

end

[Yoshio]

ダミーのデータで良いので、

F1, F2 ,Fn,C をn=25として定義し、for i = 1:nで動くコードを乗せて頂けますか？

その後、高速化を考えましょう。

[yuuji yamada]

下記にダミーデータで動作するコードを書いて見ました。よろしくお願いいたします。

num=25;

G=eye(64);

F1=rand(64, num);

F2=rand(64, num);

Pn=rand(2, num);

B = rand(64, 1);

T = zeros(num, 1, 2);

for i=1:num

Xp = [F1(:, i), F2(:, i)];

T(i, :, :) = ((diag(Pn(:, i)) * inv(Xp' * G * Xp) * Xp' * G) * B)';

end

[Yoshio]

ダミーデータ作成ありがとうございます。

細かい話で恐縮ですが、Gが単位行列となっていますが、Xp' * G * Xpを計算しているのは、Gを単位行列以外の場合にも対応させたいということでしょうか？

[yuuji yamada]

御回答ありがとうございます。

当方、該当処理部分の数式の意味を理解していなくて恐縮なのですが

Ｇの内容は固定で最初にeye(64)で定義した状態から変更されることはないです。

これで回答になっていますでしょうか。

Answer 1

[Yoshio]

ダミーのコードを参照して、以下のようにしてみました。Tは3次元配列ですが、計算結果だけからだと2次元でもよさそうなので、Tnewとしてこれを求めています。

行列を求めるところを陽に公式を適用してベクトル化していますが、今の所arrayfunが使えるるコードが浮かびません。

clear 
rng
num=250000; m = 64;
%G=eye(m);
F1=rand(m, num);
F2=rand(m, num);
Pn=rand(2, num);
BB = rand(m, 1);
T = zeros(num, 1,  2);
Tnew = zeros(num,2);
tic 
for i=1:num
    %Xp = [F1(:, i), F2(:, i)];
    %T(i, :, :) = ((diag(Pn(:, i)) * inv(Xp' * G * Xp) * Xp' * G) * BB)'
    %Tnew(i,:) = ((diag(Pn(:, i)) * inv(Xp' * Xp) * Xp') * BB)' 
    %Tnew(i,:) =  ((diag(Pn(:, i)) * ((Xp'*Xp)\Xp'))* BB)';
    %Tnew(i,:) = BB*Xp*((diag(Pn(:, i)) * inv(Xp' * Xp))'
    f1 = F1(:, i);
    f2 = F2(:, i);
    A = f1'*f1;
    B = f1'*f2;
    C = f1'*f2;
    D = f2'*f2;
    E = [D -B;-C A]/(A*D-B*C);
    a = Pn(1,i);
    b = Pn(2,i);
    Tnew(i,:) = BB'*[f1 f2]*([a*E(1,:);b*E(2,:)])';
end
toc

[Yoshio]

上記プログラムでループの中を見直したところ、Eの逆行列をベクトルで明示的に表すことで、ループ無しで書けました。25万行計算、2.3 GHz Intel Core i5のMac Miniでと0.4秒程度です。

tic
A = sum(F1.*F1);
B = sum(F1.*F2);
C = B;
D = sum(F2.*F2);
% E = [D -B;-C A]/(A*D-B*C);
Det =  A.*D-B.*C;
Tnew3 = [Pn(1,:).*(BB'*(D.*F1-C.*F2)./Det); Pn(2,:).*(BB'*(-B.*F1+A.*F2)./Det)]';
max(max(abs(Tnew-Tnew3)))
toc

[yuuji yamada]

御回答ありがとうございます。

さっそく試してみます。ありがとうございました。