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[ 
               {p, 1, P}];  

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:

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

 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.

  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – Kuba
    Jun 23, 2017 at 5:06

1 Answer 1


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]}];
    Parallelization -> True, RuntimeAttributes -> {Listable}

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

pp=30; nn=3000;

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}

  • $\begingroup$ 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. $\endgroup$
    – Fred
    Jun 23, 2017 at 21:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.