MATLAB Subtract Matrices More Efficiently

Question 1

While I still hold onto the fact that you should refactor your data, you can also use the power of datenum to help you. The following will make your life a lot easier, than trying to manually enter the days:

clear all; clc;

% Init Data
dataAll = rand(1437,159,1258);

% Starting date of sampling.  Note that this assumes each day there was a sample, and only one sample
startDate = datenum('01-01-2010');
dateList = [0:1257] + startDate;
[yr, mn, ~, ~, ~, ~] = datevec(dateList);

% extract data depending on month that sample was taken
jan = dataAll(:,:,mn == 1);
feb = dataAll(:,:,mn == 2);
mar = dataAll(:,:,mn == 3);
... % and so forth

This should get the desired results you're looking for. From here you can do your resulting calculations:

dataAnom_1 = bsxfun(@minus, jan, dataClim(:,:,1));

With your updated information from the comments you can do the following to separate your data according to month and year:

jan2010 = dataAll(:,:,(mn == 1 & yr == 2010);
feb2010 = dataAll(:,:,(mn == 2 & yr == 2010);
... % and so forth

Question 2

I think you need to rethink how you approaching the way your structuring your data here. Instead of creating massive single row arrays (dataAnom1 and dataAll) I would use a better structured matrix.

At it's simplest you could used a scheme like this: 31x12xNumYears which would produce something like this:

Data = NaN(31x12xNumYears); % Blank Init
Data(1:31,1,1) = Rand(31,1); % January's populated values
Data(1:28,2,1) = Rand(28,1); % February's populated values
... % and so forth

The advantage here is that matrix operations is far easier to do, and you have a better understanding of what the data actually represents. For example, given that dataClim is the monthy average over 30 years (matrix should be 12x30) and dataAll is the daily reading (matrix should be 31x12x30) you can do the following:

subValues = NaN(31,12,30);
for yr = 1:30
    for mn = 1:12
        subValues(1:31,mn,yr) = dataAll(1:31,mn,yr) - dataClim(mn,yr);
    end
end

With the added information that you've given me, I think this is the type of structure that you might be looking for: days x months x years x 3 where three represents the lat, long, and dataValue of your data. So for example:

test = rand(31,12,30,3);
lat = test(1:end,1:end,1:end,1); 
long = test(1:end,1:end,1:end,2);
data = test(1:end,1:end,1:end,3);

Question 3

Couldn't you do

% January
dataAnom1 = bsxfun(@minus, dataAll(:,:,[1:31 363:393 728:758 1094:1124]), dataClim(:,:,1));

? I think this is the same.

And if so, then you can do

dataAnom1=zeros(size(dataAll,1),size(dataAll,2), 128*12);
for v=1:12
    dataAnom1(:,:,1+((v-1)*128:v*128)) = bsxfun(@minus, dataAll(:,:,[1:31 363:393 728:758 1094:1124]+(v-1*32)), dataClim(:,:,v));
end

(Indexing might be a little off)