schnelle Berechnung der quadratischen Form in Matlab

Question

Ich habe den folgenden Fall:

Y(i) - m x 1 vecotr , i = 1,...,N 
    A(i) - m x m symmetric matrix , i = 1,...,N 
    H(i,j) = 0.5*(Y(i)-Y(j))'(A(i)^-1+A(j)^-1)(Y(i)-Y(j)) |i,j = 1,...,N 

Derzeit I das Inverse von A (i) separat berechnet und H mit zwei 'für' Schleifen:

for i= 1:N 
     A_inv(:,:,i) = inv(A(:,:,i)); 
    end 

    H= zeros(N,N); 
    for j=1:N 
     for i=(j+1):N 
      x = Y(:,1,j)-Y(:,1,i); 
      H(j,i) = 0.5*(x'*(A_inv(:,:,i)+A_inv(:,:,j))*x); 
      H(i,j) = H(j,i); 
     end 
    end 

Ich habe keinen Weg gefunden, den Code zu optimieren, die Antworten, die ich in den Foren gesehen habe, betreffen Fälle, in denen A konstant ist und nicht von den Indizes abhängt. Gibt es eine effizientere Möglichkeit, es zu berechnen?

Answer 1

Ich weiß nicht, welche Größen beteiligt sind, aber die folgende Vorgehensweise sollte fast immer beide schneller und genauer, da es sich nicht um eine Matrix Inversion erfordert (was etwas ist, sollte man immer zu vermeiden versuchen):

clearvars 

% Build sample inputs 

N = 100; 
m = 500; 

Y = randn(m, 1, N); 
for i = N:-1:1 
    t = randn(m); 
    A(:,:,i) = t' * t; 
end 

% original method 

tic 
for i= 1:N 
    A_inv(:,:,i) = inv(A(:,:,i)); 
end 

H= zeros(N,N); 
for j=1:N 
    for i=(j+1):N 
     x = Y(:,1,j)-Y(:,1,i); 
     H(j,i) = 0.5*(x'*(A_inv(:,:,i)+A_inv(:,:,j))*x); 
     H(i,j) = H(j,i); 
    end 
end 
toc 

H_old = H; 

clear A_inv H x 


% proposed method 

tic 

for i = N:-1:1 
    Ai_invY = A(:,:,i) \ Y(:,:); 
    H(i,:) = Y(:,i)' * Ai_invY(:,i) ... % Yi' * Ai_inv * Yi 
      - 2 * Y(:,i)' * Ai_invY ... % - Yi' * Ai_inv * Yj - Yj' * Ai_inv * Yi 
      + sum(Y(:,:) .* Ai_invY, 1); % + Yj' * Ai_inv * Yj 
end 
H = (H + H')/2; 

toc 

% check difference 

plot(H_old(:) - H(:)) 

In meinem Laptop bekomme ich 10,3 Sekunden für die ursprüngliche Methode und 0,44 Sekunden für die vorgeschlagene Methode.