What is the most efficient way to write 'for' loops in Matlab?

Question

I have read that if, for example, I have a double for loop that runs over the indexes of a matrix, then putting the column running index in the outer loop is more efficient. For example:

a=zeros(1000);
for j=1:1000
 for i=1:1000
  a(i,j)=1;
 end
end

What is the most efficient way to code it if I have three or more for loops?

For example:

a=zeros(100,100,100);
for j=1:100
 for i=1:100
  for k=1:100
   a(i,j,k)=1;
  end
 end
end

Answer 1

Short answer, you want to have the leftmost index on the innermost loop. In your example, the loop indices would go k, j, i and the array indices would be i, j, k. This has to do with how MATLAB stores the different dimensions in memory. For more, see #13 of this reddit post.

Answer 2

A somewhat longer answer that explains why it's more efficient to have the left most index varying most rapidly. There are two key things that you need to understand.

First, MATLAB (and Fortran, but not C and most other programming languages) stores arrays in memory in "column major order." e.g. if A is a 2 by 3 by 10 matrix, then the entries will be stored in memory in the order

A(1,1,1)

A(2,1,1)

A(1,2,1)

A(2,2,1)

A(1,3,1)

A(2,3,1)

A(1,1,2)

A(2,1,2)

...

A(2,3,10)

This choice of column major order is arbitrary- we could just easily adopt a "row major order" convention, and in fact that's what is done in C and some other programming languages.

The second important thing that you need to understand is that modern processors don't access memory one location at a time, but rather load and store "cache lines" of 64 or even 128 contiguous bytes (8 or 16 double precision floating point numbers) at a time from memory. These chunks of data are temporarily stored in a fast memory cache and written back out as needed. (In practice the cache architecture is now quite complicated with as many as 3 or 4 levels of cache memory, but the basic idea can be explained with a one-level cache of the sort that computers had in my younger days.)

Now, suppose that $A$ is an array with 10,000 rows and columns, and I'm looping over all of the entries.

If the loops are nested so that the innermost loop updates the row subscript, then the array entries will be accessed in the order A(1,1), A(2,1), A(3,1), ... When the first entry A(1,1) is accessed, the system will bring a cache line containing A(1,1), A(2,1), ..., A(8,1) into the cache from main memory. The next 8 iterations of the innermost loop work on this data without any additional main memory transfers.

If in the alternative, we structure the loops so that the column index varies in the innermost loop, then the entries of A would be accessed in the order A(1,1), A(1,2), A(1,3), ... In this case, the first access would bring A(1,1), A(2,1), ..., A(8,1) into the cache from main memory, but 7/8 of these entries wouldn't be used. The access to A(1,2) in the second iteration would then bring another 8 entries in from main memory, and so on. By the time the code got around to working on row 2 of the matrix, the A(2,1) entry might well be flushed out of the cache to make way for other needed data. As a result, the code is generating 8 times as much traffic as necessary.

Some optimizing compilers are capable of automatically restructuring loops to avoid this problem.

Many numerical linear algebra algorithms for matrix multiplication and factorization can be optimized to work efficiently with the row-major or column-major ordering scheme depending on the programming language. Doing this the wrong way can have a significant negative impact on performance.