In computer science, symmetric matrices can be utilized to store distances between objects or represent as adjacency matrices for undirected graphs. The main advantage of using a symmetric matrix in comparison with a classic matrix lies in smaller memory requirements. In this post, a Python implementation of such a matrix is described.

### Introduction

A matrix is called symmetric if \(a_{ij}\) is equal to \(a_{ji}\). Thanks to this rule, an \(N \times N\) symmetric matrix needs to store only \((N + 1) \cdot \frac{N}{2}\) elements instead of \(N^2\) elements needed to be stored in case of a classic matrix. An example of such a matrix is shown below.

The matrix diagonal can be seen as a mirror. Every element above this mirror is reflected to an element under this mirror. Therefore, the elements above the diagonal do not have to be stored. As mentioned previously, symmetric matrices can be used to represent distance or adjacency matrices. In the following part of this post, a Python implementation of a symmetric matrix is explained step by step along with its usage.

### Creating a Matrix

In this and subsequent sections, I show a particular usage first and then I show the implementation. The following source code shows how to create a \(4 \times 4\) symmetric matrix:

>>> from symmetric_matrix import SymmetricMatrix as SM >>> sm = SM(4)

To make this code runnable, the `SymmetricMatrix`

class has to be implemented. For now, only one special method has to be written, particularly the `__init__()`

method, which takes a single parameter called `size`

. This parameter specifies the number of rows. There is no need to pass the number of columns since symmetric matrices are square. `__init__()`

firstly checks if the provided `size`

is valid. If it is not, the `ValueError`

exception is raised. Otherwise, `size`

of the matrix is stored and the data storage for the matrix, a list in this case, is initialized. The following code shows the implementation:

class SymmetricMatrix: def __init__(self, size): if size <= 0: raise ValueError('size has to be positive') self._size = size self._data = [0 for i in range((size + 1) * size // 2)]

It is worth noting the size of the `_data`

storage used to store the matrix. It is smaller than \(size^2\). To explain the computation of the number of elements, suppose that we have a \(N \times N\) symmetric matrix. To save space, only elements under and on the diagonal need to be saved. Therefore, for the first row only one element has to be stored, for the second row two elements are saved and so on. Accordingly, for the \(N\)-th row, \(N\) elements need to be saved. If we sum all elements that need to be saved from all rows, we get the following result:

$$1 + 2 + \cdots + N = (1 + N) \cdot \frac{N}{2}$$

### Getting Matrix Size

It would be nice to have a possibility to use a standard Python way for gaining the matrix size, which is the `len()`

function. Therefore, to obtain the matrix size, we wish that the following code could be used:

>>> print(len(sm)) 4

To actuate the previous code, another magic method has to be implemented. This method is `__len__()`

and its only responsibility is to return the `_size`

attribute:

def __len__(self): return self._size

### Accessing the Matrix

Until now, we were able to create a symmetric matrix with all elements initialized to zero and get its size. However, this is not very useful in real life. We also need to write to and read from the matrix. Since we want the usage of the matrix be as much comfortable and natural as possible, the subscript operator `[]`

will be used when accessing the matrix:

>>> sm[1, 1] = 1 >>> print('sm[1, 1] is', sm[1, 1]) sm[1, 1] is 1

#### Writing

Firstly, let us focus on writing to the matrix. In Python, when an assignment to `sm[1, 1]`

is executed, the interpreter calls the `__setitem__()`

magic method. To achieve the expected behaviour, this method has to be implemented in `SymmetricMatrix`

. Since only elements under and on the diagonal are stored and the whole matrix is saved in a one-dimensional data storage, a correct index to this storage needs to be calculated. For now, assume that the `_get_index()`

method returns this index. Later on, the implementation of this method will be shown. Now, when we have the index, we can use the `__setitem__()`

method provided by the underlying storage that can be called simply as `self._data[index] = value`

:

def __setitem__(self, position, value): index = self._get_index(position) self._data[index] = value

#### Reading

For obtaining an element from the matrix, we will proceed in a similar way. Therefore, another magic method, particularly the `__getitem__()`

method, has to be implemented. Similarly as in the previous case, to get the desired element from the matrix, the position has to be converted to a proper index to the underlying storage. This service is done by the `_get_index()`

method for which the last part of this section is devoted. When we have the correct index, the element on this position in the underlying storage is returned:

def __getitem__(self, position): index = self._get_index(position) return self._data[index]

#### Calculating Indexes

Now, it is time to show how `_get_index()`

is implemented. The passed position is a pair of the form `(row, column)`

. The source code of this method can be broken down into two steps that have to be executed in the provided order:

- convert a position above the diagonal into a proper position below the diagonal and
- calculate the correct index into the underlying storage.

If the given position, `(row, column)`

, is above the diagonal, then `row`

is swapped with `column`

, since every element above the diagonal has its counterpart exactly at the `(column, row)`

position. To clarify the second part, particularly the calculation of the index into the used storage, the above picture and the following table will be used:

Row | Size of All Previous Rows | Matrix Position | Calculated Index |
---|---|---|---|

\(1\) | \(0\) | \((0, column)\) | \(0 + column\) |

\(2\) | \(0 + 1\) | \((1, column)\) | \(0 + 1 + column\) |

\(3\) | \(0 + 1 + 2\) | \((2, column)\) | \(0 + 1 + 2 + column\) |

… | … | … | … |

\(row + 1\) | \(0 + 1 + 2 + 3 + \cdots + row\) | \((row, column)\) | \(0 + 1 + 2 + 3 + \cdots + row + column\) |

Note that for the first row, the column part of the `(row, column)`

pair is sufficient to use as index to the underlying storage. For the second row, the number of elements in the previous row and column part of the `(row, column)`

pair is enough. In the case of the second row, the calculated index is \(1 + column\), since the previous row contains only one element. For the third row, the situation is a little bit complicated because the elements from all the previous rows have to be summed. So, the index for the `(2, column)`

position is \(1 + 2 + column\). Therefore, for the `(row, column)`

position the correct index is \(1 + 2 + 3 + \cdots + row + column\). Thanks to the finite arithmetic progression, this expression can be simplified as follows:

$$0 + 1 + 2 + \cdots + row + column = (0 + row) \cdot \frac{row + 1}{2} + column$$

Finally, the implementation of calculating the index into the underlying storage is shown in the following source code:

def _get_index(self, position): row, column = position if column > row: row, column = column, row index = (0 + row) * (row + 1) // 2 + column return index

### Using an Own Type of a Storage

Now, we have a working implementation of a symmetric matrix. Since the main motivation for using this type of matrix is memory efficiency, the question that may emerged is if a more memory efficient implementation can be made. This leads us to think if the used `list`

is the best data structure for the storage. If you are familiar with the Python implementation of `list`

, you may know that `list`

does not contain elements that you insert into it. Indeed, it contains pointers to these elements. Hence, the memory requirements are higher for `list`

than, for example, for `array.array`

that stores the elements directly. Of course, there are other data structures that are more memory efficient than `list`

. So, the question is which one should be used. Suppose that we chose `array.array`

instead of `list`

during the symmetric matrix implementation. Later, this matrix needs to be shared between several processes. Certainly, it will not work since `array.array`

is not supposed to be shared by different processes.

Therefore, a better solution when choosing the underlying data structure is leaving space for users to choose the type of the storage according to their requirements. If no special demands are present then `list`

can be used as the default storage type. Otherwise, the user passes his storage type during the matrix creation like in the following example:

>>> def create_storage(size): ... return multiprocessing.Array(ctypes.c_int64, size) ... >>> sm = SM(3, create_storage) >>> sm[1, 2] = 5

The above `create_storage()`

returns an array holding 64b integers that can be shared by different processes.

To implement this improvement only small changes are necessary in the `__init__()`

method. Firstly, one parameter, namely `create_storage`

, is added with default value set to `None`

. If an argument for this parameter is not passed, then `list`

will be used as the storage type. Otherwise, a function that takes one parameter, particularly the size of the storage, and returns the created storage is expected:

def __init__(self, size, create_storage=None): if size <= 0: raise ValueError('size has to be positive') if create_storage is None: create_storage = lambda n: [0 for i in range(n)] self._size = size self._data = create_storage((size + 1) * size // 2)

### Benchmark

To provide a comparison between the introduced symmetric matrix and a matrix created via the `numpy`

module, I have written a benchmark script that uses a \(4000 \times 4000\) matrix to show memory requirements and average access times for the implemented symmetric matrix and the `numpy`

matrix. When creating a symmetric matrix, `array.array()`

is used as the underlying storage. To create the `numpy`

matrix, `numpy.zeros()`

is called. The elements in both matrices are 64b integers.

#### Memory Usage

Firstly, memory usage is compared. The `asizeof.asizeof()`

function from the `pympler`

module computes the sizes of the created matrices. The result of this experiment can be seen in the table below.

Matrix Type | Memory Usage |
---|---|

Symmetric Matrix (via `array` ) |
61.05 MB |

`numpy` Matrix |
122.07 MB |

We can see that the symmetric matrix can save approximately 50% of memory space.

#### Access Time

Next, access times for writing to the entire matrix are computed for both matrix types. This computation is performed five times and then the average result is calculated. The experiments ran on an Intel Quad-Core i7-4700HQ (6M Cache, 2.40 GHz) processor. From the following table, we can see that the average access time for the implemented symmetric matrix is much worse than the average access time for the `numpy`

matrix:

Matrix Type | Access Time |
---|---|

Symmetric Matrix (via `array` ) |
11.26 sec |

`numpy` Matrix |
2.00 sec |

The reasons behind the slow access time for the symmetric matrix can be revealed by the `cProfile`

module. Before running the script with the `cProfile`

module, only the relevant parts were present. Therefore, the first part comparing memory requirements and all parts using the `numpy`

code are not included in the profiling.

$ python -m cProfile -s cumtime benchmark.py ... Ordered by: cumulative time ncalls tottime percall cumtime percall filename:lineno(function) 130/1 0.005 0.000 74.098 74.098 {built-in method exec} 1 0.000 0.000 74.098 74.098 benchmark.py:8(<module>) 1 0.002 0.002 74.028 74.028 benchmark.py:91(main) 1 25.421 25.421 74.026 74.026 benchmark.py:56(benchmark_access_time) 80000000 24.473 0.000 48.290 0.000 symmetric_matrix.py:22(__setitem__) 80000000 23.817 0.000 23.817 0.000 symmetric_matrix.py:30(_get_index) 1 0.000 0.000 0.315 0.315 symmetric_matrix.py:9(__init__) 1 0.315 0.315 0.315 0.315 benchmark.py:21(create_array) ...

For understanding the above output, only three columns are important for us, namely `ncalls`

, `cumtime`

and `filename:lineno(function)`

. The first one, named `ncalls`

, represents how many times the function from `filename:lineno(function)`

was called. The `cumtime`

column informs us about the cumulative time spent in this function and all sub-functions during all calls. As can be seen from the output, the time is spent mostly in `__setitem__()`

and `_get_index()`

. The overhead is due to internal workings of Python and computing indexes to the underlying storage. Thus, this symmetric matrix implementation is suitable in circumstances where memory usage is a bigger problem than processor power.

### Source Code

The complete source code of the implemented `SymmetricMatrix`

class, alongside with unit tests and the benchmark script, is available on GitHub. All code was written, tested and profiled in Python 3.4.