matlab中二次型的快速计算

Question

我有以下情况：

    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

目前我用两个 'for' 循环分别计算了 A(i) 和 H 的逆：

    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

我还没有找到优化代码的方法，我在论坛上看到的答案是针对 A 不变且不依赖于索引的情况。有没有更有效的计算方法？

Answer 1

我不知道涉及什么大小，但以下方法应该几乎总是更快和更准确，因为它不需要矩阵求逆（这是我们应该始终尽量避免的事情）：

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(:))

在我的笔记本电脑上，原始方法的时间为 10.3 秒，建议的方法为 0.44 秒。

matlab中二次型的快速计算

Fast calculation of quadratic form in matlab

performance

matlab

quadratic