# Speed up computation of sum from large matrix

I have a 3x3 matrix (252^3) of data (densities) and I want to compute a correlation xi from nearest neighbours as well as a sum involving a check whether the density is in a certain bin. The computation works in principle (so for subsamples, when step size is about 10) but is painfully slow. I am not an expert on handling data with Mathematica, so I am pretty sure there is some solution that works way better than my first try. For the correlation function I think that basically need an analogue of AbsoluteCorrelationFunction for my matrix.

xi=
1/3*Sum[data[[i,j,k]] (data[[i + 1,j,k]] +
data[[i,j + 1,k]] +
data[[i,j,k + 1]]), {i, 1, 251, 1}, {j, 1,
251, 1}, {k, 1, 251, 1}]/26^3 - 1;


For the following quantity I am not sure if there is another faster Mathematica routine. I need to check whether the density is in a bin of width deltarho around rho and if yes add the values of the neighbouring densities

brho[rho_]:=(Sum[If[Abs[data[[i,j,k]] - rho] <
deltarho/2, (data[[i + 1,j,k]] +
data[[i,j + 1,k]] +
data[[i,j,k + 1]]), 0], {i, 1, 251, 1}, {j, 1,
251, 1}, {k, 1, 251, 1}]/
Sum[Boole[
Abs[data[[i,j,k]] - rho] < deltarho/2], {i, 1,
251, 1}, {j, 1, 251, 1}, {k, 1, 251, 1}] - 1)


I tried replacing Sum with ParallelSum but at least for the subsample it makes things slower.

Are there any ideas how to speed this up?

edit: improved readability & corrected typo, thanks to comment of Jason

• Just first impressions - to improve readability in the code for yourself and others, you can replace things like data[[i]][[j]][[k]]  with data[[i,j,k]] and , {i, 1, 251, 1}, {j, 1, 251, 1}, {k, 1, 251, 1}] with ,{i, 251}, {j, 251}, {k, 251} but this is just a style point and doesn't address your main question. – Jason B. Apr 15 '16 at 9:35
• I tried running your definition of xi on a random matrix dims = {251, 251, 251}; data = RandomReal[1, dims]; and it takes a long time. You could get the same thing virtually instantly by a judicious use of Total and RotateLeft and Drop – Jason B. Apr 15 '16 at 9:52
• Also, Can I assume that {i, 1, 251, 10} is a typo in the brho code? It's the only 10 in there so I assume it should be a 1 – Jason B. Apr 15 '16 at 10:06

Well, basically your code is designed with a less efficient algorithm.

Bear in mind that Mathematica generally treats a matrix the "same" as a number, so don't spend time on computation conducted at the number level.

Then we can convert your code to be running at the matrix level, where the speed is much increased.

len = 20;
data = RandomReal[{1, 10}, {len, len, len}];
xi=1/3*Sum[data[[i,j,k]] (data[[i+1,j,k]]+data[[i,j+1,k]]+data[[i,j,k+1]]),{i,len-1},{j,len-1},{k,len-1}]/26^3-1;//AbsoluteTiming
(*{0.030215, Null}*)
xxi=Total@Flatten@(data[[1;;len-1,1;;len-1,1;;len-1]](data[[2;;len,1;;len-1,1;;len-1]]+data[[1;;len-1,2;;len,1;;len-1]]+data[[1;;len-1,1;;len-1,2;;len]]))/3/26^3-1;//AbsoluteTiming
(*{0.00039, Null}*)
xi == xxi
(*True*)


Now we increase len to 252 as required by your scenario:

len = 252;
data = RandomReal[{1, 10}, {len, len, len}];
xxi=Total@Flatten@(data[[1;;len-1,1;;len-1,1;;len-1]](data[[2;;len,1;;len-1,1;;len-1]]+data[[1;;len-1,2;;len,1;;len-1]]+data[[1;;len-1,1;;len-1,2;;len]]))/3/26^3-1;//AbsoluteTiming
(*{0.653753, Null}*)


See, the performance is not bad.

I think now you can improve the second part of your problem with the idea above.

Edit: As for the second part of your problem, there is no need to perform testing over each element: a holistic testing works as well:

Clear[brho]
brho[data_List,rho_,deltarho_]:=Module[{c,sum},
c=UnitStep[deltarho/2-Abs[data[[1;;-2,1;;-2,1;;-2]]-rho]];
sum=data[[2;;-1,1;;-2,1;;-2]]+data[[1;;-2,2;;-1,1;;-2]]+data[[1;;-2,1;;-2,2;;-1]];
Total@Flatten[c*sum]/Total@Flatten[c]
]


Give it a try:

brho[data,10,3]