How to store an symmetrical matrix?

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!