Implementing a Symmetric Matrix in Python

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.

Symmetric Matrix

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.

Storing Symmetric Matrix in One-Dimensional 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.

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\)

sopticek