処理の高速化
2 views (last 30 days)
Show older comments
以前、処理の高速化ができないものかと思い、下記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をLL'と分解して計算が簡略化できる可能性があるのですが、Gの例はあるでしょうか?
サンプルが無いと検証もできませんので、よろしくお願いします。
Answers (1)
Yoshio
on 4 Nov 2019
Edited: 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
See Also
Categories
Find more on Resizing and Reshaping Matrices 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!