# Parallelized Algorithms for computing large Similarity Matricies

Fast-similarity-matrix repository

In this post I’m going to demonstrate three algorithms for computing large symmetric similarity matrices (in this case, 20000 x 20000). The first is a baseline approach that we will use as a benchmark to measure the performance of the latter two. Algorithm 1 and 2 compute the matrix by distributing the work across multiple parallel workers. Algorithm 1 achieves a **6x speed up**. Algorithm 2 achieves a **14.3x speed up** but is more memory-intensive than algorithm 1. I use python’s `concurrent.futures`

, `concurrent.futures.ProcessPoolExecutor()`

,
and the computation was performed on a Ryzen 9, 12 core, 24 thread processor, 3.8GHz base.

Algorithm | Description | Time elapsed (s) | relative speed | speed up |
---|---|---|---|---|

0 | loop through each pair | 1376 | 1 | 0 |

1 | send inner loop to workers | 229 | 0.1664 | 6.00 |

2 | sum component-wise distance matricies | 96 | 0.0698 | 14.33 |

# Parameters used for all algorithms

I will be computing the similarity matrix using a synthetic dataset of 20000 uniformily distributed points in 4-dimensional space.

```
N = 20000 #number of data points
D = 4 #dimensionality of data
tau = 0.05 #optional similarity cutoff
dat=np.random.rand(N,D)
```

The `tau`

will be used for creating sparsity in the matrix so that it can be stored in sparse matrix format. This is optional but will reduce the memory used to store the matrix.

# Baseline algorithm: cycle through all pairs of data points

See `Fast_similarity_matrix.ipynb`

In this naive approach, we simply loop through all elements in the lower triangular part of the matrix and record the distance separating each pair of points. This assumes the matrix should be symmetric, which is the case for most types of similiarity matrices. It also assumes the elements on the diagonal are not needed.

```
W = np.zeros((N,N), dtype=float) #similarity matrix to be filled in
#fill lower triangular part only
for i in np.arange(1,N):
for j in np.arange(i):
W[i,j] = F.dist(dat[i,:],dat[j,:])
```

We then apply the gaussian similarity function and threshold the resulting similarity values by `tau`

, effectively setting the most distant points to have zero similarity. This obviously discards some information on the global structure of the data. The effect of doing this is the subject of another post! The benefit of creating sparsity is that we can reduce the amount of memory used to store the matrix. Here I use the scipy COOrdiate sparse matrix format.

```
sigm = lambda x : np.exp(-x**2/(2*0.3**2)) #gaussian similarity function or kernel
Q = sigm(W)
Q[Q==1]=0 #set upper triangle to 0 (it should hold 1's currently)
#threshold to increase sparsity, sacrificing some information
Qthresh = Q.copy()
Qthresh[Qthresh<=tau] = 0
#convert to COOrdinate format sparse matrix
Qthresh_coo = coo_matrix(Qthresh)
```

### Memory use

Matrix | description | size (.nbytes) (GB) | size (sys.getsizeof) (bytes) |
---|---|---|---|

Q | dense similarity matrix | 3.2 | 3.2E+09 |

Q_coo | sparse version of Q, coo format | 3.19984 | 48 |

Qthresh | thresholded dense matrix | 3.2 | 3.2E+09 |

Qthresh_coo | sparse version of Qthresh | 1.37335 | 48 |

# Algorithm 1: Each processor computes the distances for a single data point

In this method we allow each CPU worker to compute the distances for a different data point. That is, the inner loop of the baseline algorithm is delegated to one worker so that the rows of the matrix are computed in parallel. This algorithm calls the method `row_similarity`

in `Fast_similarity_matrix_func.py`

. The function and the data needed by each CPU worker is distributed using `concurrent.futures.ProcessPoolExecutor()`

.
To do this we must create lists of each parameter and a corresponding list of copies of the data set: `[dat for _ in range(N)]`

.

```
ROW = list(np.arange(1,N))
TAU = [tau for _ in range(N)]
```

Each worker will receive one element of these lists at the same index. There may be a way to avoid creating a list of copies of the dataset, but I haven’t gone that far with this method. `ROW`

contains the list of row indices; each worker will receive one row to process and will use the `row_similarity`

method shown later.

```
data=[]
rows=[]
cols=[]
with concurrent.futures.ProcessPoolExecutor() as executor:
results_generator = executor.map(F.row_similarity, [dat for _ in range(N)], ROW, TAU) #
#append each peice of data sent back from the workers to a list of arrays
for i in results_generator:
if i[0]: #if nonempty
data.append(np.array(i[0]))
rows.append(np.array(i[1]))
cols.append(np.array(i[2]))
# print(data)
data = np.concatenate(data, axis = 0)
rows = np.concatenate(rows, axis = 0)
cols = np.concatenate(cols, axis = 0)
Q_PPEmap = coo_matrix((data, (rows, cols)), shape=(N,N))
```

The `data`

,`rows`

, and `cols`

(used by the scipy coo format) seem to be returned in the order they are completed. We can simply append all incoming data returned by `row_similarity`

together and later convert it to coo matrix format. Below are the functions given to each worker.

```
def dist(x, y):
return np.sqrt(np.sum((x-y)**2))
def sigm(d):
return np.exp(-d**2/(2*0.3**2))
def row_similarity(dat, row, tau):
#dat: data matrix
#row: what row index this worker is in charge of looping across
#tau: similarity threshold
row_arr = []
col_arr = []
data = []
for j in np.arange(row):
if sigm(dist(dat[row, :], dat[j, :])) >= tau:
data.append(sigm(dist(dat[row, :], dat[j, :])))
row_arr.append(row)
col_arr.append(j)
return data, row_arr, col_arr
```

Note that the value of tau may affect the speed of this since we do not append similarity values below tau. A totally fair comparison with the baseline algorithm would use `tau = 0`

, but the majority of the speed up should be due to the parallelization and not to skipping some pairs of points.

# Algorithm 2: Component-wise computation

In this method we split the computation into \(D\) components, where \(D\) is the dimensionality of the data. Each worker receives one dimemsion and computes a distance matrix on that dimension only. These component-wise distance matrices are ultimately added together to obtain the full distance matrix. The method is much more memory intensive since each CPU processor must do matrix operations on large matricies. The speed up in this case comes from using vectorized operations.

Let \(O\) be the data matrix where rows are \(N\) observations in \(D\) dimensions. We will tile the \(j^{th}\) column of \(O\) horizonally to form matricies \(C\) and \(C^T\)

\(j^{th}\) column of \(O\):

\[x_j = [O_{1,j}, O_{2,j}, ...,O_{N,j}]^T\]\(C_j\) (a tiling of \(N\) repetitions):

\[C_j = [x_j, x_j, ..., x_j]\]The distance matrix can then be computed as

\[W = \sqrt{\sum_{j=1}^D (C_j - C_j^T)^2}\]We again use the `concurrent.futures.ProcessPoolExecutor()`

. However, this method will not utilize all processors unless the dimensionality of the data exceeds the processor count. Run `multiprocessing.cpu_count()`

to check. This time we call the `componentwise_similarity`

method.

```
Nlist = [N for _ in range(D)]
Dlist = [d for d in range(D)]
M = np.zeros((N,N), dtype=float)
with concurrent.futures.ProcessPoolExecutor() as executor:
result_generator = executor.map(F.componentwise_similarity, [dat for _ in range(D)], Dlist, Nlist)
for d in result_generator:
M += d
sigm = lambda x : np.exp(-x**2/(2*0.3**2))
M = sigm(np.sqrt(M))
M[M<tau]=0
M[M==1]=0
QM = coo_matrix(np.tril(M))
```

From the module `Fast_similarity_matrix_func.py`

handed to each worker:

```
def componentwise_similarity(dat, dim, N):
#dim: dimension of input vectors to take the difference in
#dat: data matrix
col_tile = np.tile(np.expand_dims(dat[:, dim], axis=1), (1, N))
row_tile = np.transpose(col_tile)
return (col_tile - row_tile)**2
```