嗯，这就是我的意思，但@Yvon，你是对的，它比将NaNs 分配给0 要慢。此外，使用 TIMEIT 功能进行更好的测试（Elkan 的答案似乎不正确）。我确信有一些小的改进和其他方法，但我想我会添加这个。
Test Results:
func4 %# Jon potential improves on func3
runTime: 0.00374437,
error:5.32907e-15
func3 %# Elkan's comment on Yvon's answer
runTime: 0.00375083,
error:5.32907e-15
func2 %# Expanded Norm Yvon
runTime: 0.00386379,
error:5.32907e-15
func1 %# Original
runTime: 0.00515395,
error:0
func5 %# Jon's no zero-assign version of func3
runTime: 0.00884188,
error:5.32907e-15
func6 %# Jon's no zero-assign version of Original (norm-mean directly)
runTime: 0.0105508,
error:2.66454e-15
代码：
function test()
a = randn(80,3000);
[r,c] = size(a);
b = randperm(r*c);
nanix = b(1:round(numel(b)*0.3)); % randomly select 30% of values to be NaN
a(nanix) = NaN;
% nnx = sum(~isnan(a));
f = {...
@() func1(a); %# Original
@() func2(a); %# Expanded Yvon
@() func3(a); %# Elkan's comment
@() func4(a); %# Jon potential improves on 3
@() func5(a); %# Jon's no zero-assign Elkan's
@() func6(a); %# Jon's no zero-assign on Original
};
%Time each function/method
timings = cellfun(@timeit, f);
validity = cellfun(@feval, f, 'UniformOutput',false);
err=max( abs( bsxfun(@minus,validity{1},vertcat(validity{:})) ),[],2 );
%Display in order of speed:
idxTime=1:size(timings,1);
[sortedTime,sortMap]=sort(timings);
fprintf('Test Results:\n\n');
fprintf('func%d\r runTime: %g,\r error:%g\n',[idxTime(sortMap);sortedTime.';err(sortMap).'])
end
%% Functions below
function xnorm = func1(X)
nanix = isnan(X);
nx = sum(~nanix); % count the number of non-NaN values in each column for calculating mean
X(nanix) = 0;
bsxminus = sum(X)./nx;
X = bsxfun(@minus, X, bsxminus);
X(nanix) = 0;
xnorm = sqrt(sum(X.^2));
end
function xnorm2 = func2(X)
nanix = isnan(X);
nx = sum(~nanix);
X(nanix) = 0;
Xsum = sum(X);
Xmean = Xsum./nx;
X2 = sum(X.^2);
Xmean2 = Xmean.^2.*nx;
XXmean = Xsum.*Xmean;
xnorm2 = sqrt( X2+Xmean2-XXmean-XXmean ); % avoid bsx
end
function xnorm3 = func3(X)
nanix = isnan(X);
nx = sum(~nanix);
X(nanix) = 0;
Xsum = sum(X);
X2sum = sum(X.^2); % 50% time consumed here
xnorm3 = sqrt( X2sum - Xsum.^2./nx ); % avoid bsx
% relative error
end
%Try column version of setting to zero, Replace element-wise powers with multiplication, add realsqrt
function xnorm4 = func4(X)
[r,c] = size(X);
nanix = isnan(X);
nx = sum(~nanix);%prob can't speedup without Mex
Xtemp=X(:);
Xtemp(nanix(:))=0;
X=reshape(Xtemp,r,c); %Sometimes faster if linear?
Xsum = sum(X);
X2sum = sum(X.*X);
xnorm4 = realsqrt( X2sum - (Xsum.*Xsum)./nx ); % avoid bsx
end
%Now using indexing instead of setting to zero
function xnorm5 = func5(X)
nnanix = ~isnan(X);
nx = sum(nnanix);
% Xtemp=X(nnanix); %40% of time lost
Xtemp=X(:);
Xtemp=Xtemp(nnanix(:)); %drops by 20%, still slower than =0
%Now have two vectors, one of values, and one of # of values per
%original column... avoid cells!
ind(1,numel(Xtemp)) = 0;%hack pre-allocation
ind(cumsum(nx(end:-1:1))) = 1;
ind = cumsum(ind(end:-1:1)); %grouping index for accumarray
Xsum = accumarray(ind(:), Xtemp(:)); %sum each "original" columns
X2sum = accumarray(ind(:), Xtemp(:).*Xtemp(:));%./nx(:) for mean directly
xnorm5 = realsqrt( X2sum - (Xsum.*Xsum)./nx(:) ).'; % avoid bsx
end
%Now using indexing instead of setting to zero
function xnorm6 = func6(X)
nnanix = ~isnan(X);
nx = sum(nnanix); %only Mex could speed up
% Xtemp=X(nnanix); %40% of time lost
Xtemp=X(:);
Xtemp=Xtemp(nnanix(:)); %drops by 20%, still slower than =0
numNoNans=numel(Xtemp);
%Now have two vectors, one of values, and one of # of values per
%original column... avoid cells! It's like a "run length encoding"
%Almost like "run-length encoding" FEX: RunLength, for MEX
ind(1,numNoNans) = 0;%hack pre-allocation
ind(cumsum(nx(end:-1:1))) = 1;
ind = cumsum(ind(end:-1:1)); %// generate grouping index for accumarray
Xmean = (accumarray(ind(:), Xtemp(:))./nx(:)).'; %means of each col
%"Run-length decoding" %Repelem is fastest if 2015b>
idx(1,numNoNans)=0;
idx([1 cumsum(nx(1:end-1))+1]) = diff([0 Xmean]);
XmeanRepeated = cumsum(idx);
XminusMean = Xtemp(:).'-XmeanRepeated; %Each col subtracted by col mean
XSumOfSquares = accumarray(ind(:), XminusMean.*XminusMean); %sum each "original" columns
xnorm6 = realsqrt(XSumOfSquares).';
end