arrayfunの使い方について

28 views (last 30 days)
yuuji yamada
yuuji yamada on 25 Sep 2019
Commented: yuuji yamada on 14 Oct 2019
現在、Matlab2018bを使用し、アプリを作成しているのですが、処理を高速化させたい部分があります。
それはforループを使用して約25万回ループしている箇所です。
arrayfun関数を使用するとforループを使用しないで書くことができるようなのですが下記のような処理の場合でも
arrayfun関数でforループを使用せず書くことができるものでしょうか。
改善したい箇所は大体、下記のような感じになっており、最終的に25万行のデータSを作成したいと思っています。
S=zeros(250000, 1, 2);
for i=1:250000
Xp = [F1(:, i), F2(:, i)];
S(i, :, :) = Fn(:, i) * Xp' * Xp * Xp' * C;
end
  4 Comments
Yoshio
Yoshio on 3 Oct 2019
ダミーデータ作成ありがとうございます。
細かい話で恐縮ですが、Gが単位行列となっていますが、Xp' * G * Xpを計算しているのは、Gを単位行列以外の場合にも対応させたいということでしょうか?
yuuji yamada
yuuji yamada on 3 Oct 2019
御回答ありがとうございます。
当方、該当処理部分の数式の意味を理解していなくて恐縮なのですが
Gの内容は固定で最初にeye(64)で定義した状態から変更されることはないです。
これで回答になっていますでしょうか。

Sign in to comment.

Accepted Answer

Yoshio
Yoshio on 8 Oct 2019
Edited: Yoshio on 8 Oct 2019
ダミーのコードを参照して、以下のようにしてみました。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
  2 Comments
Yoshio
Yoshio on 14 Oct 2019
Edited: Yoshio on 14 Oct 2019
上記プログラムでループの中を見直したところ、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
yuuji yamada on 14 Oct 2019
御回答ありがとうございます。
さっそく試してみます。ありがとうございました。

Sign in to comment.

More Answers (0)

Categories

Find more on ループと条件付きステートメント in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!