処理の高速化

22 views (last 30 days)
yuuji yamada
yuuji yamada on 30 Oct 2019
Edited: 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
yuuji yamada
yuuji yamada on 31 Oct 2019
御回答ありがとうございます。
Gの構造は64x64の行列になります。
Yoshio
Yoshio on 31 Oct 2019
実対称行列であれば、GをLL'と分解して計算が簡略化できる可能性があるのですが、Gの例はあるでしょうか?
サンプルが無いと検証もできませんので、よろしくお願いします。

Sign in to comment.

Answers (1)

Yoshio
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

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!