Why Do Row-Major/Column-Major and Iteration vs. Indexing Tests Yield Unexpected Performance in Python?

The reason why people usually say that is because of cache.

When a computer access to memory at address 1000000, it also access memore at address 1000001, 1000002, 1000003, ..., 1000255 (typical cache line size), and stores it in cache, which is way faster than usual RAM.

If, later, it needs what is at address 1000100, that access will be faster, since it will find it in the previous cache line.
Unless, in the meantime, that cache line was removed from the cache. Since each time it writes a cache line in the cache, it is removing a previous one ; cache size is limited (typically a few MB) ; probably the older one.

So, the more often you use memory that is already in cache, the faster.

And that occurs if you access your data in the same order they are stored in memory.

We can try it with an example

PYTHONCopied
import numpy as np
from numba import njit
import timeit

@njit 
def somRow(M):
    s=0
    h,w=M.shape
    for i in range(h):
        for j in range(w):
            s+=M[i,j]
    return s

@njit 
def somCol(M):
    s=0
    h,w=M.shape
    for j in range(w):
        for i in range(h):
            s+=M[i,j]
    return s

M=np.random.uniform(size=(5000,5000))
print("M")
print("som row =", somRow(M))
print("som col =", somCol(M))
print("som row timing", timeit.timeit(lambda: somRow(M), number=100))
print("som col timing", timeit.timeit(lambda: somCol(M), number=100))

P=M.T
print("P")
print("som row =", somRow(P))
print("som col =", somCol(P))
print("som row timing", timeit.timeit(lambda: somRow(P), number=100))
print("som col timing", timeit.timeit(lambda: somCol(P), number=100))

On my computer, it takes 2.5 seconds for somRow and 8 seconds for somCol (for 100 iterations, not just one), on matrix M. While on matrix P, it takes 8 seconds for somRow and 2.5 for somCol (matrix P is stored in column major order, since P[i+1,j] is what immediately follows P[i,j] in memory. Because P, aka M.T use just the same memory as M, with different strides)

So cache effect is not negligible. Its a factor more than 3 on performance.

And that answers to "why do people say that we need to process in row major order data that are stored in row major order"

Note that this is just an application of the more general rule "it is more efficient to process data in the order they are in memory".

Another example that does not involved matrices and "row/columns" order

PYTHONCopied
import numpy as np
from numba import njit

T=np.random.uniform(size=(25000000,))

@njit
def som1(T):
    s=0
    for i in range(25000000):
        j=(i*25000001)%25000000
        s+=T[j]
    return s

@njit
def som2(T):
    s=0
    for i in range(25000000):
        j=(i*25001547)%25000000
        s+=T[j]
    return s

print("som =", som1(T), som2(T))
# The two sums are identical, because 
print("time1", timeit.timeit(lambda: som1(T), number=100))
print("time2", timeit.timeit(lambda: som2(T), number=100))

I use a bizarre way to compute index j. In som1, j is just i, but written in a way that forces (I hope) the compiler to do the actual computation, so that timing for som2 is not unjustly biased by the j computation.
In som2, j is not i, but since 1547 and 25000000 are prime together, the sum remains the same: all indices are iterated, just, not in the canonical order (which is also the order in which they are in memory).

Now some remarks about this code that are relevant to your question

1. I need to compute the sum a first time before I start timing. For two reasons. The first one is because I use numba, and first computation in numba includes compilation of the code. I don't want to include compilation time in my timings. But secondly, and more importantly to you, because if I didn't, my result would be biased in favor of the second computation, whatever it is. Because for the second computation, memory has already been read for the first computation, and is therefore more likely to be already in cache. Especially if the array size is not that big. See my previous answer to another question = reason for the timing difference experienced by asker was simply that the second computation was faster ; just swapping the order of the experiments changed completely the conclusions of the timings !

Your experiment fails to address that bias. You perform two experiments on the same data (index vs direct iteration). That is a bias against the first one, that have to put the data in cache (note that I doubt that this is the reason of the result. But, well, when you perform an experiment that is based upon the very fact that cache effect are not negligible, since that is the postulate of the book you read, you should make sure that the same cache effect does not alter the experiment result!).

2. Since I mentioned numba, you may wonder why I use it. There are two reasons. The first one is because I need to be sure of the order in which data are accessed. With numpy function you cannot be sure. Most of some are more complex that just a sum. And more importantly, numpy is often smarter than I am; it is aware of those cache effect, and often use the best order.
3. The other reason is because if I didn't use numba, what I would have timed would have been the pure python overhead of the computation. The sum that is taking 2.5 seconds vs 8 seconds (for 100 computations) on my computer, takes 450 seconds for both without numba (so 4.5 seconds per iteration). Cache effect becomes negligible before the elephant in the room: python is a VERY slow language. The art of making efficient code in python is about not coding in python, but having most of the important code run in the internal code of some fast libraries (not coded in python). Like numpy for example. If you iterate your numpy array with pure python for loop (with or without index) you've already lost the performance challenge. The trick with numpy is to always find way to avoid for loops and rely on numpy functions.

Said otherwise: it is completely irrelevant to wonder about row order and column order when you code in pure python. That is optimizing milliseconds in code that takes dozens of seconds.

And that is also the reason why the "without index" code is faster. What you are timing are pure python overhead. Your code is only made of pure python overhead. And the cost of maintaining a pure python variable j is way more important than any cache consideration that your book is referring to.

You are using Python.

Most of the execution time of your program is spent by the Python interpreter.

Memory access costs (cache misses etc) are negligible compared to that.

In summary, the measurement did not and cannot measure what you wanted to measure.

---

With compiled languages, or using numba in Python, you have a chance of seeing how one access pattern causes more cache misses, which would make it slower than the other pattern.