向量化回归图计算答案

【问题标题】：Vectorize a regression map calculation向量化回归图计算
【发布时间】：2018-12-30 04:51:56
【问题描述】：

我在字段B(x,y,t) 上计算时间序列A(t) 的回归图，方法如下：

A=1:10; %time
B=rand(100,100,10); %x,y,time

rc=nan(size(B,1),size(B,2));
for ii=size(B,1)
  for jj=1:size(B,2)
     tmp = cov(A,squeeze(B(ii,jj,:))); %covariance matrix
     rc(ii,jj) = tmp(1,2); %covariance A and B
  end
end
rc = rc/var(A); %regression coefficient

有没有办法对代码进行矢量化/加速？或者也许一些我不知道的内置函数可以达到相同的结果？

【问题讨论】：

标签： matlab performance regression vectorization covariance

【解决方案1】：

为了将此算法向量化，您必须“亲自动手”并自己计算协方差。如果您查看cov 内部，您会发现它有很多行输入检查和很少的实际计算行，总结一下关键步骤：

y = varargin{1};
x = x(:);
y = y(:);
x = [x y];
[m,~] = size(x);
denom = m - 1;
xc = x - sum(x,1)./m;  % Remove mean
c = (xc' * xc) ./ denom;

为了稍微简化一下上面的内容：

x = [x(:) y(:)];
m = size(x,1);
xc = x - sum(x,1)./m;
c = (xc' * xc) ./ (m - 1);

现在这是相当简单的矢量化...

function q51466884
A = 1:10; %time
B = rand(200,200,10); %x,y,time
%% Test Equivalence:
assert( norm(sol1-sol2) < 1E-10);
%% Benchmark:
disp([timeit(@sol1), timeit(@sol2)]);

%%
function rc = sol1()
rc=nan(size(B,1),size(B,2));
for ii=1:size(B,1)
  for jj=1:size(B,2)
     tmp = cov(A,squeeze(B(ii,jj,:))); %covariance matrix
     rc(ii,jj) = tmp(1,2); %covariance A and B
  end
end
rc = rc/var(A); %regression coefficient
end

function rC = sol2()  
m = numel(A);
rB = reshape(B,[],10).'; % reshape
% Center:
cA = A(:) - sum(A)./m;
cB = rB - sum(rB,1)./m;
% Multiply:
rC = reshape( (cA.' * cB) ./ (m-1), size(B(:,:,1)) ) ./ var(A);
end

end

我得到了这些时间：[0.5381 0.0025]，这意味着我们在运行时节省了两个数量级：)

请注意，优化算法的很大一部分是假设您的数据中没有任何“奇怪”，例如 NaN 值等。查看 cov.m 内部以查看我们跳过的所有检查。

【讨论】：