処理の高速化

以前、処理の高速化ができないものかと思い、下記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)';

[Yoshio]

Gに関して、対角行列であるとか、対称行列であるなど構造上の特徴はありますか？

[yuuji yamada]

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

Gの構造は64ｘ64の行列になります。

[Yoshio]

実対称行列であれば、ＧをLL'と分解して計算が簡略化できる可能性があるのですが、Gの例はあるでしょうか？

サンプルが無いと検証もできませんので、よろしくお願いします。

Answer 1

[Yoshio]

ベクトル化前のコードも入れていますので、参考とされてください。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