matrix storage, many rows or many columns?

tl;dr Should I store a matrix of Nx3 elements row-wise or column-wise, if I'm going to access 3 elements at a time? Does it matter?

I'm setting up a numerical simulation with matrices of x,y,z positions in 3 dimensions, and I wonder whether it's better to store them row-wise or column-wise in R, Matlab, and Armadillo (C++). Maybe it doesn't matter at all.

The simulation considers N particles, say N=5000, and each is described by 3 coordinates, 3 angles, and 3 sizes. So I've 3 matrices of Nx3 elements. In the code (setting up a linear system of 3N equations), I will need to access positions, sizes, and angles one particle at a time. Therefore, I'm tempted to store those as columns, i.e. those matrices would have 3 rows, and N columns. This is because R, Matlab and Armadillo all have a column-major ordering, so my thinking is that accessing a column will be (ever so slightly) faster. On the other hand, maybe there is a penalty of having many columns (memory addresses)? I do not know.

By-the-way, I'm mentioning three programming languages because I'll be writing the code in all three (it's not much effort, since the syntax is pretty much the same).

• This question is sufficiently broad that it's not really specific to any language or matrix library - it's mostly about the memory hierarchy and how memory access patterns affect performance. – Kirill Nov 14 '15 at 1:39

3 Answers

Which storage scheme is optimal depends very strongly on the access pattern. You already said that you're accessing them one particle at a time, but are you going to access the particles sequentially in order, or randomly? Another conventional name for this issue is struct-of-arrays vs array-of-structs.

The key thing here is to know how the memory hierarchy works on a CPU (see, e.g., https://www.scss.tcd.ie/Jeremy.Jones/CS3021/5%20caches.pdf, as well as the many questions on this topic on Stack Overflow). Reads from main memory have a very large latency, which will lead to long stalls in the execution, unless this is mitigated by prefetching. Each read from main memory will read one cache line worth of data, placing it in cache, and the reads from cache are comparatively free.

So for example if you have a matrix $(a_{ij})_{1\leq i,j\leq 3}$, and it is stored in column-major order in memory as $$a_{11}, a_{21}, a_{31}, a_{12}, a_{22}, a_{32}, a_{13}, a_{23}, a_{33},$$ then reading $a_{11}$ (or any of $a_{11},\ldots,a_{23}$) will place all the eight elements from $a_{11}$ to $a_{23}$ in cache (assuming correct alignment and a typical 64-byte cache line and an element size of 8 bytes), making them free to access.

If you have several different arrays, a struct-of-arrays, like $$x_{1}, \ldots, \qquad y_{1}, \ldots, \qquad z_{1}, \ldots$$ and you read $x_1, y_1, z_1$, you will effectively read the first eight elements of each array. If you end up using the next seven elements immediately, that's fine. If not, then that wasted 7/8 of your memory bandwidth by reading 24 elements instead of 3, which is inefficient.

If instead you store them as array-of-structs, $$x_1, y_1, z_1, x_2, y_2, z_2, \ldots$$ then reading $x_1,y_1,z_1$ will read the first eight elements. If you use those next elements, there is no great difference, otherwise you waste 5/8 (reading 8 instead of 3) of your memory bandwidth - note this is still better than reading 24 elements as in the previous example.

With nine 8-byte elements, if you know you'll access them consecutively there should be no big difference either way, as above. If you will access them in random order, then reading one set of nine elements will lead to either $9\times 8 = 72$ elements being read from memory (if stored as struct-of-arrays) or $\lceil \frac{9}{8}\rceil\times 8 = 16$ elements (if stored as array-of-structs) - this could be a very large performance gain.

The difference between column-major and row-major, I feel, is a bit of red herring here - the key is knowing the access pattern, and making sure the data layout matches the access pattern. The next most important thing would be prefetching, which is a technique for hiding the relatively large memory access latency.

A separate related issue to consider is how you access data that you already read previously, and whether it is likely to still remain in the cache or whether it was purged to make way for other data, causing cache misses - an array of $5000\times 9\times 8$ bytes is $360\,\mathrm{KB}$ and will likely fit entirely in either the L2 or L3 cache.

It gets even more complicated if your program involves multiple threads running on different cores, with potentially separate lower-level caches per core.

Note, incidentally, that for non-sequential access patterns, having nine elements per particle is somewhat inefficient in itself if a cache line is 8 elements long as reading nine elements would be essentially equivalent to reading 16 elements - one technique is to split your data into often-used versus rarely-used parts so that the often-used parts fit better.

There is no penalty for having either many rows or columns. Not sure what you mean by a penalty for "many memory addresses", but there is no such penalty.

Some of the things I said depend on exact properties of the architecture your program will run on. Using an interpreted language like R or Matlab can also affect this severely - for C++ it should hold. In general, it is best to use the array-of-structs pattern, unless you have some specific reasons not to.

• thanks a lot for this thorough discussion. The particles are accessed sequentially, so I guess my current strategy of 3 rows, N columns is better than the alternative. I also have a set of N matrices of dimension 3x3, and I was wondering whether storing them in a cube (3D array) would make sense, but it doesn't sound like it because each slice of 9 elements wouldn't fit in the cache (and they're complex values, to make things worse). – baptiste Nov 14 '15 at 3:40
• @baptiste Re your last sentence: cache line size (typically 64 bytes) is different from cache size - cache will contain many cache lines (e.g. L1 cache might be 64KB = 1024 cache lines). Cache lines are like units into which the cache is divided. Units of data can easily occupy multiple cache lines - my remark only meant that for random access patterns if there is an "overhang" then some memory bandwidth might be wasted. – Kirill Nov 14 '15 at 4:31

R, Matlab, and Armadillo (Mat class) matrices use column-major order so you should create 3xN matrices so that each set of 3 values is stored consecutively in memory.

A bit more info:
You should store your matrices such that elements which are accessed together are stored sequentially in memory. The memory order of matrix storage is language, data structure, and possibly even implementation dependent. The two storage orders used for matrices are referred to as "row-major" or "column-major" order. In row-major order elements in the same row are stored consecutively in memory and vice versa in column-major.

The advantage of choosing the correct arrangement is decreased cache misses and decreased memory latency. Efficient cache usage can sometimes lead to large speedups.

A final note: C and C++ 2D arrays are row-major, but Armadillo's Mat class is column-major since the storage layout was designed to be compatible with BLAS and LAPACK which use Fortran's column-major storage order. Matlab derives much of it underpinnings and syntax from Fortran so it naturally uses column-major storage order. I'm unsure of the history behind R's storage order, but perhaps someone will comment with that information.

• thanks, that confirms my preferred choice of columns-first – baptiste Nov 14 '15 at 3:50

Column-major order seems to be more natural. For example suppose if you want to save movie to file picture by picture then you are using column order, and that is very intuitive and nobody would save it in row-major order.

If you are programmer in C/C++ you should use some higher level libraries for matrices (Eigen, Armadillo,...) with default column-major order. Only maniac would use raw C pointers with row-major order, although C/C++ offers something that reminds matrix indexing.

For simplicity everything with row-major order should be considered as at least strange formed. Slice by slice is simply natural order and it means column-major order (like Fortran). Our fathers/mothers had a very good reasons why they chose it.

Unfortunately before it became clear several interesting libraries were created in row-major order, probably due to lack of experience.

To clarify recall the definition of row-major order where right index vary faster in one step through memory eg A(x,y,z) it is z-index, it means that in memory pixels from different slices are adjacent, what we wouldn't want. For movie A(x,y,t) the last index is time t. It is not hard to imagine that it is simply impossible to save movie in row-major mode.