Efficient calculations with symmetric matrix with effective memory

Question

I need to do a lot of calculations like this:

X.dot(Y).dot(Xt)

X = 1 xn matrix

Y = symmetric nxn matrix, each element has one of 5 values ​​(for example, 0, 0.25, 0.5, 0.75, 1)

Xt = nx 1, transpose of X, i.e. X[np.newaxis].T

X and Y are dense. I have a problem for large n, I can not store and use the matrix Y due to memory problems. I am limited to using a single machine, so distributed computing is not an option.

It occurred to me that Y has 2 functions that could theoretically reduce the amount of memory needed to store Y:

• Elements Y are covered by a small list of values.
• Y is symmetrical.

How can I put this into practice? I was looking for a repository of symmetric matrices, but as far as I know, all the number multiplications of a matrix require "unpacking" the symmetry to create a regular nxn matrix.

I understand that numpy is for in-memory computing, so I am open to alternative python-based solutions that are not dependent on numpy. I am also willing to sacrifice speed for memory efficiency.

Answer 1

UPDATE: I found using a format that fills 3 matrix elements in one byte, actually quite quickly. In the example below, the speed penalty is less 2xcompared to direct multiplication using @, while the space savings are greater than 20x.

>>> Y = np.random.randint(0, 5, (3000, 3000), dtype = np.int8)
>>> i, j = np.triu_indices(3000, 1)
>>> Y[i, j] = Y[j, i]
>>> values = np.array([0.3, 0.5, 0.6, 0.9, 2.0])
>>> Ycmp = (np.reshape(Y, (1000, 3, 3000)) * np.array([25, 5, 1], dtype=np.int8)[None, :, None]).sum(axis=1, dtype=np.int8)
>>> 
>>> full = values[Y]
>>> x @ full @ x
1972379.8153972814
>>> 
>>> vtable = values[np.transpose(np.unravel_index(np.arange(125), (5,5,5)))]
>>> np.dot(np.concatenate([(vtable * np.bincount(row, x, minlength=125)[:, None]).sum(axis=0) for row in Ycmp]), x)
1972379.8153972814
>>> 
>>> timeit('x @ full @ x', globals=globals(), number=100)
0.7130507210385986
>>> timeit('np.dot(np.concatenate([(vtable * np.bincount(row, x, minlength=125)[:, None]).sum(axis=0) for row in Ycmp]), x)', globals=globals(), number=100)
1.3755558349657804

. .

, np.bincount :

>>> Y = np.random.randint(0, 5, (1000, 1000), dtype = np.int8)
>>> i, j = np.triu_indices(1000, 1)
>>> Y[i, j] = Y[j, i]
>>> values = np.array([0.3, 0.5, 0.6, 0.9, 2.0])
>>> full = values[Y]
>>> # full would correspond to your original matrix,
>>> # Y is the 'compressed' version
>>>
>>> x = np.random.random((1000,))
>>>
>>> # direct method for reference 
>>> x @ full @ x
217515.13954751115
>>> 
>>> # memory saving version
>>> np.dot([(values * np.bincount(row, x)).sum() for row in Y], x)
217515.13954751118
>>>
>>> # to save another almost 50% exploit symmetry
>>> upper = Y[i, j]
>>> diag = np.diagonal(Y)
>>> 
>>> boundaries = np.r_[0, np.cumsum(np.arange(999, 0, -1))]
>>> (values*np.bincount(diag, x*x)).sum() + 2 * np.dot([(values*np.bincount(upper[boundaries[i]:boundaries[i+1]], x[i+1:],minlength=5)).sum() for i in range(999)], x[:-1])
217515.13954751115

Answer 2

Y, ​​ numpy.array int, @PaulPanzer, , : 27 64 , 64/log2 (5) = 27,56...

:

import numpy as np

row = np.random.randint(5, size=100)

# pad with zeros to length that is multiple of 27
if len(row)%27:
    row_pad = np.append(row, np.zeros(27 - len(row)%27, dtype=int))
else:
    row_pad = row

row_compr = []
y_compr = 0
for i, y in enumerate(row_pad):
    if i > 0 and i % 27 == 0:
        row_compr.append(y_compr)
        y_compr = 0
    y_compr *= 5
    y_compr += y

# append last 
row_compr.append(y_compr)
row_compr = np.array(row_compr, dtype=np.int64)

:

row_decompr = []
for y_compr in row_compr:
    y_block = np.zeros(shape=27, dtype=np.uint8)
    for i in range(27):
        y_block[26-i] = y_compr % 5
        y_compr = int(y_compr // 5)
    row_decompr.append(y_block)

row_decompr = np.array(row_decompr).flatten()[:len(row)]

Y:

assert np.allclose(row, row_decompr)