# Fast manipulation of large dimensional arrays

I am dealing with a multiplication of two- and three- dimensional arrays. Since the arrays are large, I am compiling a function to deal with it. The code of the function, with a working example, is as follows:

(*array dimensions*)
pp = 30; nn = 3000;

(*sample data*)
a = Table[1., {p, 1, pp}, {n1, 1, nn}, {n2, 1, nn}];
b = Table[1., {p, 1, pp}, {n1, 1, nn}];
c = Table[1., {p, 1, pp}, {n1, 1, nn}];
d = Table[1., {p, 1, pp}, {n1, 1, nn}];

(*defining compiled function*)
With[{P = pp}, (*localizing P to avoid MainEvaluate*)
y=Compile[{{a, _Real, 3}, {b, _Real, 2}, {c, _Real, 2}, {d, _Real, 2}},
Module[{num, den, i1, out},
i1 = b c d ;
num = a i1;
den = Table[Total[num[[p]]], {p, 1, P}];
out = Table[
Transpose[Transpose[num[[p]]]/den[[p]]],
{p, 1, P}];
out
]
]
];


All inputs are real and positive numbers. The dimensions of the input variables are as follows:

• a: pp x nn x nn
• b,c,d: pp x nn

The variable P is also compiled within the module, and there are no MainEvaluate in the code.

To be clear, the purpose of the last instruction in the compiled function is to:

• take "page p" of num (I am thinking of a 3-dimensional array as a book, with pages, and tables NxN in each page);
• divide each column by the correspondent row of den, to obtain again an NxN matrix (this is why I am using a double Transpose instruction).
• create pp pages of such matrices.

In my application, pp~30 and nn~3000. It takes about 9-10 seconds to run this function on my (fairly fast) computer:

 DateString[]
test = y[a, b, c, d];
DateString[]

Out[363]= "Thu 22 Jun 2017 09:10:43"

Out[365]= "Thu 22 Jun 2017 09:10:54"


Since I will need to run it thousands of times, this is a pretty inefficient code. I am looking for a way to substantially speed up the computation, possibly avoiding the Transpose functions.

I have found that the most expensive part it the last instruction, that seems to take around 4 seconds by itself. A couple of seconds go before the first instruction and after the last computation. Three-four seconds are taken by the first three instructions.

Any suggestions are greatly appreciated.

• Comments are not for extended discussion; this conversation has been moved to chat. – Kuba Jun 23 '17 at 5:06

Using parallelization as suggested by @aardvark2012 and avoiding Transpose provides a nice speed up:

y2 = Compile[{{a, _Real, 2}, {bcd, _Real, 1}},
Module[{num = a (bcd), den},
den = Total[num];

Do[num[[i]] /= den, {i, Length[a]}];
num
],
Parallelization -> True, RuntimeAttributes -> {Listable}
];


Instead of comparing performance on unpacked arrays of of 1s, I will use packed random arrays:

pp=30; nn=3000;
SeedRandom[1];

a = RandomReal[1, {pp, nn, nn}];
b = RandomReal[1, {pp, nn}];
c = RandomReal[1, {pp, nn}];
d = RandomReal[1, {pp, nn}];


Now, for a timing/memory comparison:

r1 = y[a, b, c, d]; //MaxMemoryUsed //AbsoluteTiming
r2 = y2[a, b c d]; //MaxMemoryUsed //AbsoluteTiming
MinMax[r1 - r2]


{17.2071, 6625489808}

{1.95376, 288822456}

{-3.46945*10^-18, 3.46945*10^-18}

• Thank you. I am finding that some of the gains in time emerge because of the use of packed arrays. Using unpacked arrays moves my computations from 10.2 to 7 seconds; declaring the variables as packed arrays moves them from 8.3 to 2.04. Overall, my speed increases with a ratio of 5 (as opposed to around 8 in your case, if I read correctly). That may be due to different processors, I imagine. – Fred Jun 23 '17 at 21:37