Why does bsxfun produce different result than brute force mean?

I would like to apply a 2D mask to a 3D matrix, then average the non-masked values. In the example below, the mean is calculated in two ways, (a) with bsxfun and (b) in a loop which sums values and divides by the count. The results aren't the same, though. Can you spot my error, which I suspect is in approach (b)?

nRow = 200; nCol = 100; nChannel = 10;
bigMatrix = rand(nRow, nCol, nChannel); % 3D matrix of values that we want to mask
mask = ones(nRow, nCol); % Make a 2D mask
mask(1:25, 26:50) = NaN; % Mask set to NaN where values should be excluded

%%(a) Use bsxfun to substitute NaNs for masked values in bigMatrix
% Since mask values are either NaN or 1, can mask by multiplying 
maskedByMult = bsxfun(@times, bigMatrix, cast(mask, 'like', bigMatrix));
% Average the non-nan values
meanByBsxfun = squeeze(mean(mean(maskedByMult, 'omitnan'), 'omitnan')); 

%%(b) Mask by loop
nGoodValues = 0; % Count of unmasked values used in later calculation of their mean
sumValues = zeros(nChannel, 1); % Sums will be added to this
for iRow = 1:length(nRow)
    for iCol = 1:length(nCol)
        if isfinite(mask(iRow, iCol)) % Check if this x,y location is masked out
            nGoodValues = nGoodValues + 1;
            sumValues = sumValues + squeeze(bigMatrix(iRow, iCol, :));
        end
    end
end
meanByLoop = sumValues/nGoodValues; % Get the mean

%%Show they differ
figure('Name', 'Compare two ways of getting the mean');
plot(1:10, meanByBsxfun); hold on;
plot(1:10, meanByLoop, '--');
legend('Using bsxfun', 'Using a loop');

Answer 1

[David Goodmanson]

Hello KAE,

I believe that each method has a small problem. In the second method, since length(nRrow) = length(nCol) = 1, the for loops only execute one time. So you need

for iRow = 1:nRow          % etc.

In the first method, denote maskedbyMult by A, and let N = the number of elements that are not nan, which in this case is 200*100 - 625. Then

mean(mean(A),'omitnan')   is not equal to    sum(sum(A),'omitnan')/N    which is what you want.

for example,

a = magic(3);  
a(1,2) = nan

a =
     8   NaN     6
     3     5     7
     4     9     2

mean(mean(a,'omitnan'))
ans = 5.6667

sum(sum(a,'omitnan'))/8
ans = 5.5000

After those changes both methods agree, and both agree with

bigMatMask = repmat(mask,1,1,10).*bigMatrix;
mean3 = squeeze(sum(sum(bigMatMask,2,'omitnan'),1)/(200*100-625))

% if you have a newer version of Matlab with explicit expansion,
% the first line can be    bigMatMask = mask.*bigMatrix;

[KAE]

Thanks so much, since I have used mean(mean(x)) erroneously in lots of other code. By the way here are the time differences on my machine,

• bxsfun (a) 0.003004 s
• loop (b) 0.108191 s
• repmat without implicit expansion since I am using R2015b (c) 0.003103 s

And here is the corrected full code if that helps anyone.

nRow = 200; nCol = 100; nChannel = 10;
bigMatrix = rand(nRow, nCol, nChannel); % 3D matrix of values that we want to mask
mask = ones(nRow, nCol); % Make a 2D mask
mask(1:25, 26:50) = NaN; % Mask set to NaN where values should be excluded
nIn = nRow*nCol - sum(isnan(mask(:))); % How many non-masked values are in matrix

%%(a) Use bsxfun to substitute NaNs for masked values in bigMatrix
tic;
% Since mask values are either NaN or 1, can mask by multiplying
maskedByMult = bsxfun(@times, bigMatrix, cast(mask, 'like', bigMatrix));
% Average the non-nan values
meanByBsxfun = squeeze(sum(sum(maskedByMult, 'omitnan'))/nIn);
toc;

%%(b) Mask by loop
tic;
nGoodValues = 0; % Count of unmasked values used in later calculation of their mean
sumValues = zeros(nChannel, 1); % Sums will be added to this
for iRow = 1:nRow
    for iCol = 1:nCol
        if isfinite(mask(iRow, iCol)) % Check if this x,y location is masked out
            nGoodValues = nGoodValues + 1;
            sumValues = sumValues + squeeze(bigMatrix(iRow, iCol, :));
        end
    end
end
meanByLoop = sumValues/nGoodValues; % Get the mean
toc;

%%(c) Repmat
tic;
bigMask = repmat(mask,1,1,nChannel).*bigMatrix;
meanByRepmat = squeeze(sum(sum(bigMask,2,'omitnan'),1)/nIn);
toc;

%%Show they are the same
figure('Name', 'Compare two ways of getting the mean');
plot(1:10, meanByBsxfun); hold on;
plot(1:10, meanByLoop, '--');
plot(1:10, meanByRepmat, ':');
legend('Using bsxfun', 'Using a loop', 'Using repmat');

[David Goodmanson]

nice to see the code comparison.

mean(mean( )) of course does work out in lots of cases such as rectangular matrices without any nans.

By the way I meant to say implicit expansion, not explicit expansion.

[KAE]

If someone wants to do another operation besides mean on a masked matrix, here is an example using median,

        % Reshape matrix into [spatial dimension x channel]
        temp = reshape(maskedByMult, [nRow*nCol nChannel]); 
        medianByBsxfun = median(temp, 'omitnan');