MATLAB Answers

0

処理の高速化

Asked by yuuji yamada on 30 Oct 2019
Latest activity Edited by Yoshio
on 6 Nov 2019
以前、処理の高速化ができないものかと思い、下記URLの投稿をしたものです。
https://jp.mathworks.com/matlabcentral/answers/482062-arrayfun
上記の投稿でGが単位行列と回答をしてしまいましたが私が勘違いをしておりGは
単位行列ではないことが判明しました。その場合、どのような記述になりますでしょうか。
御回答して頂いたコードにGを足して下記のように書き換えるところまではできたのですが
その先で難儀しています。どう書きかえればよいでしょうか。
%Xp = [F1(:, i), F2(:, i)];
%T(i, :, :) = ((diag(Pn(:, i)) * inv(Xp' * G * Xp) * Xp' * G) * BB)'
                  ↓
%Tnew(i,:) = ((diag(Pn(:, i)) * inv(Xp' * G * Xp) * Xp' * G) * BB)'
%Tnew(i,:) = ((diag(Pn(:, i)) * ((Xp' * G * Xp)\(Xp' * G)))* BB)';

  3 Comments

Yoshio
on 31 Oct 2019
Gに関して、対角行列であるとか、対称行列であるなど構造上の特徴はありますか?
yuuji yamada on 31 Oct 2019
御回答ありがとうございます。
Gの構造は64x64の行列になります。
Yoshio
on 31 Oct 2019
実対称行列であれば、GをLL'と分解して計算が簡略化できる可能性があるのですが、Gの例はあるでしょうか?
サンプルが無いと検証もできませんので、よろしくお願いします。

Sign in to comment.

1 Answer

Yoshio
Answer by Yoshio
on 4 Nov 2019
Edited by Yoshio
on 6 Nov 2019

ベクトル化前のコードも入れていますので、参考とされてください。Gが対角行列の場合は、ベクトル化しても誤差は無視できそうですが、そうで無い場合は、誤差が大きくなる場合がありそうですので、あくまで参考にお考えください。
逆行列の公式をそのまま使うのは、数値計算上は望ましくないので、誤差評価は十分行ってください。
clear
rng
num=250000; m = 64;
%G = diag(rand(m,1));
G = rand(m);
F1=rand(m, num);
F2=rand(m, num);
Pn=rand(2, num);
BB = rand(m, 1);
T = zeros(num, 1, 2);
T0 = zeros(num,2);
Tnew1 = T0;
Tnew2 = T0;
Tnew3 = T0;
tic
for i=1:num
Xp = [F1(:, i), F2(:, i)];
%T(i, :, :) = ((diag(Pn(:, i)) * inv(Xp' * G * Xp) * Xp' * G) * BB)';
% T0(i,:) = ((diag(Pn(:, i)) * inv(Xp' * G * Xp) * Xp' * G) * BB)';
% T1(i,:) = ((diag(Pn(:, i)) * inv(Xp'*G*Xp)* Xp'*G) * BB)' ;
% T2(i,:) = ((diag(Pn(:, i)) * ((Xp'*G*Xp)\(Xp'*G)))* BB)';
% T3(i,:) = BB'*(G'*Xp*inv(Xp'*G*Xp)'*diag(Pn(:, i)));
f1 = F1(:, i);
f2 = F2(:, i);
g1 = G'*f1;
g2 = G'*f2;
A = g1'*f1;
B = g1'*f2;
C = g2'*f1;
D = g2'*f2;
Det = A*D-B*C;
a = Pn(1,i);
b = Pn(2,i);
Tnew1(i,:) = BB'*[g1 g2]*[a*[D;-B] b*[-C;A]]/Det;
Tnew2(i,:) = BB'*[a*(D*g1-B*g2) b*(-C*g1+A*g2)]/Det;
end
toc
% max(max(abs(T0-T1)))
% max(max(abs(T0-T2)))
% max(max(abs(T1-T3)))
% max(max(abs(T0-Tnew1)))
% max(max(abs(T0-Tnew2)))
max(max(abs(Tnew1-Tnew2)))
tic
G1 = G'*F1;
G2 = G'*F2;
A = sum(G1.*F1);
B = sum(G1.*F2);
C = sum(G2.*F1);
D = sum(G2.*F2);
Det = A.*D-B.*C;
Tnew3 = [Pn(1,:).*(BB'*(D.*G1-B.*G2)./Det); Pn(2,:).*(BB'*(-C.*G1+A.*G2)./Det)]';
max(max(abs(Tnew2-Tnew3)))
toc

  0 Comments

Sign in to comment.