[SciPy-User] efficient way to store and use a 4D redundant matrix

Emanuele Olivetti emanuele@relativita....
Wed Mar 16 10:41:07 CDT 2011

```Hi Everybody,

I have a 4D matrix "A" where entries differing only by a permutation of
indices are equal, i.e. A[1,2,3,4] = A[3,4,2,1]. Since each set of
of 4 indices I have then 24 permutations, my goal would be to avoid
storing 24 times the necessary space. I am looking at scipy.sparse but
I have no clear understanding whether it could handle my case.

Any idea?

For sake of clarity here is a simple code computing a matrix A like
the one above. I usually have A of shape (100,100,100,100) which requires
~800Mb when the dtype is double. Of course the non-redundant part is
just ~4Mb so you may understand my interest in this issue.

----
import numpy as np
from itertools import permutations
try:
form itertools import combinations_with_replacement
except ImportError: # if Python < v2.7 that function is not available...
form itertools import product
def combinations_with_replacement(iterable, r):
pool = tuple(iterable)
n = len(pool)
for indices in product(range(n), repeat=r):
if sorted(indices) == list(indices):
yield tuple(pool[i] for i in indices)

rows = 10
columns = 20
x = np.random.rand(rows,columns)+1.0
A = np.zeros((rows, rows, rows, rows))

indices_rows = range(rows)
for i,j,k,l in combinations_with_replacement(indices_rows, 4):
tmp = (x[i,:]*x[j,:]*x[k,:]*x[l,:]).sum()
for a,b,c,d in permutations([i,j,k,l]):
A[a,b,c,d] = tmp
----

In case you wonder which kind of operations do I need to do on A,
they are usual manypulations (slicing, reshaping, etc.), and common
linear algebra.

Best,

Emanuele

```