[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
     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.



More information about the SciPy-User mailing list