12
$\begingroup$

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
$\endgroup$
  • 4
    $\begingroup$ For loops are very slow in MATLAB. You should avoid explicit loops in MATLAB whenever possible. Instead, usually a problem can expressed in terms of matrix/vector operations. That is the MATLABic way. There are also a lot of built-in functions to initialise matrices, etc. For instance, there is a function, ones(), that will set all elements of a matrix to 1 (by extension, to any value by multiplication (a scalar multiplied by the all-ones matrix)). It also works on 3-D arrays (which I think covers the example here). $\endgroup$ – Peter Mortensen Jul 28 at 16:33
  • 3
    $\begingroup$ @PeterMortensen By what factor (roughly) the efficiency of loops in Matlab is smaller compared to C and Python? And why is that? Also, havn't the efficiency of loops in Matlab gotten better in the last several years? $\endgroup$ – TensoR Jul 28 at 17:56
  • 3
    $\begingroup$ @PeterMortensen “usually a problem can be expressed in terms of matrix/vector operations” – for certain values of “usually”, yes. IMO it's more accurate to say that people working in Matlab and the like have a decades-long culture of ignoring all the things that can't be done with matrix/vector operations, so much that everything looks like a nail to them for that hammer. And we shouldn't just say “for loops in Matlab are slow” but “Matlab is slow” (it only happens to be linked to a fast library of LA primitives written in C and Fortran). $\endgroup$ – leftaroundabout Jul 29 at 9:33
  • 5
    $\begingroup$ The performance of for loops is controversial: matlabtips.com/matlab-is-no-longer-slow-at-for-loops $\endgroup$ – ohreally Jul 29 at 10:58
  • $\begingroup$ @leftaroundabout True. Being concerned about speed in an interpreted (or semi-interpreted) language is a pretty clear indication you have an X-Y problem where the actual solution is "don't use this language". The exception of course is if you're using code generation in Simulink, but then the question is what C the code generator builds and how efficient that is. $\endgroup$ – Graham Jul 30 at 9:04
18
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ Or use the built-in function ones(). $\endgroup$ – Peter Mortensen Jul 28 at 16:28
  • 4
    $\begingroup$ @Peter OP's example is almost certainly just a toy example of a for loop that does something and not the actual use case. $\endgroup$ – Matt Jul 28 at 17:09
  • $\begingroup$ @Matt You are correct. $\endgroup$ – TensoR Jul 28 at 17:57
11
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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