Hi,

for my C++ course at may university, we were asked to write a small jet analysis for some LHC data. While writing mine, I was confronted with this: How do you store a symmetrical $N \times N$ matrix efficiently in C++?

If you just store it as a normal matrix, the needed space will go with $N^2$. Also, every update has to be done twice. Not very nice (and that does not even take any problem with multidimensional arrays in C++ in to account)

$$\begin{pmatrix} a_{00} & a_{01} & a_{02} & a_{03} \\ a_{10} & a_{11} & a_{12} & a_{13} \\ a_{20} & a_{21} & a_{22} & a_{23} \\ a_{30} & a_{31} & a_{32} & a_{33} \end{pmatrix}$$

Now, what do we really need for this matrix?

$$\begin{pmatrix} a_{00} & \color{gray}{a_{01}} & \color{gray}{a_{02}} & \color{gray}{a_{03}} \ a_{10} & a_{11} & \color{gray}{a_{12}} & \color{gray}{a_{13}} \ a_{20} & a_{21} & a_{22} & \color{gray}{a_{23}} \ a_{30} & a_{31} & a_{32} & a_{33} \end{pmatrix}$$

Which does not look like its much but for simplicity sake I am only using an $N=4$-matrix. For my analysis, $N$ was around $10000$. So how much elements do we really need?

The first row has $1$ element, the second $2$ and so on up to the last row which has $N$ elements. Therefore we have a total of $\sum_{i=1}^N i = \frac{N^2 + N}{2}$ elements. Great! We basically need only half the elements and therefore only half the RAM!

But how do we store it? We need a way to transform this multidimensional array into a single dimensional array and then use a simple function which transforms the 2D coordinates i,j to an index p of our single dimensional array. Lets take a look at the matrix again:

$$\begin{pmatrix} \color{red}{a_{00}} & \color{gray}{a_{01}} & \color{gray}{a_{02}} & \color{gray}{a_{03}} \ \color{green}{a_{10}} & \color{green}{a_{11}} & \color{gray}{a_{12}} & \color{gray}{a_{13}} \ \color{blue}{a_{20}} & \color{blue}{a_{21}} & \color{blue}{a_{22}} & \color{gray}{a_{23}} \ \color{orange}{a_{30}} & \color{orange}{a_{31}} & \color{orange}{a_{32}} & \color{orange}{a_{33}} \end{pmatrix}$$

Now, then we write this rows just after each other like this:

$$\begin{pmatrix} \color{red}{a_{00}} & \color{green}{a_{10}} & \color{green}{a_{11}} & \color{blue}{a_{20}} & \color{blue}{a_{21}} & \color{blue}{a_{22}} & \color{orange}{a_{30}} & \color{orange}{a_{31}} & \color{orange}{a_{32}} & \color{orange}{a_{33}} \end{pmatrix}$$

and use a reference function like this: (we will require $i \geq j$ for this but if $j > i$, you can just switch both as the matrix is symmetrical)

$$p = j + \sum_{k=0}^i k = j + \frac{i^2 + i}{2}$$

Yeah! We now can allocate an single dimension array and use it as an symmetrical matrix without worrying about multiple elements!