numpy multidimensional indexing and diagonal symmetries

Question

I have a pretty large numpy array ...

power = ...
print power.shape
>>> (3, 10, 10, 19, 75, 10, 10)

which is symmetric w.r.t. the 10x10 parts, i.e the following 2-d matrices are symmetric

power[i, :, :, j, k, l, m]
power[i, j, k, l, m, :, :]

for all values of i, j, k, l, m

Can I utilize this factor-4 gain? E.g. when saving the matrix to file (50 mb with savez_compressed)

My attempt:

size = 10
row_idx, col_idx = np.tril_indices(size)
zip_idx = zip(row_idx, col_idx)
print len(zip_idx), zip_idx[:5]
>>> 55 [(0, 0), (1, 0), (1, 1), (2, 0), (2, 1)]
all_idx = [(r0, c0, r1, c1) for (r0, c0) in zip_idx for (r1, c1) in zip_idx]
print len(all_idx), all_idx[:5]
>>> 3025 [(0, 0, 0, 0), (0, 0, 1, 0), (0, 0, 1, 1), (0, 0, 2, 0), (0, 0, 2, 1)]
a, b, c, d = zip(*all_idx)
tril_part = np.transpose(s.power, (0, 3, 4, 1, 2, 5, 6))[:,:,:, a, b, c, d]
print tril_part.shape
>>> (3, 19, 75, 3025)

This seems ugly, but "works" ... once I can also get back to power from tril_part ...

I guess this yields two qurestions:

1. Better way of going from power to tril_part?
2. How to go from tril_part to power?

Edit: The "size" comment is clearly valid, but please ignore it :-) IMHO the indexing part of the question stands alone. I've been finding myself wanting to do similar indexing for smaller matrices.

Answer 1

You are on the right path. Using np.tril_indices you can indeed smartly index these lower triangles. What remains to be improved is the actual indexing/slicing of the data.

Please try this (copy and pasteable):

import numpy as np
shape = (3, 10, 10, 19, 75, 10, 10)
p = np.arange(np.prod(shape)).reshape(shape)  # this is not symmetric, but not important

ix, iy = np.tril_indices(10)
# In order to index properly, we need to add axes. This can be done by hand or with this
ix1, ix2 = np.ix_(ix, ix)
iy1, iy2 = np.ix_(iy, iy)

p_ltriag = p[:, ix1, iy1, :, :, ix2, iy2]
print p_ltriag.shape  # yields (55, 55, 3, 19, 75), axis order can be changed if needed

q = np.zeros_like(p)
q[:, ix1, iy1, :, :, ix2, iy2] = p_ltriag  # fills the lower triangles on both sides
q[:, ix1, iy1, :, :, iy2, ix2] = p_ltriag  # fills the lower on left, upper on right
q[:, iy1, ix1, :, :, ix2, iy2] = p_ltriag  # fills the upper on left, lower on right
q[:, iy1, ix1, :, :, iy2, ix2] = p_ltriag  # fills the upper triangles on both sides

The array q now contains a symmetrized version of p (where the upper triangles were replaced with the content of the lower triangles). Note that the last line contains iy and ix indices in inversed order, essentially creating a transpose of the lower triangular matrix.

Comparison on lower triangles In order to compare back, we set all the upper triangles to 0

ux, uy = np.triu_indices(10)
p[:, ux, uy] = 0
q[:, ux, uy] = 0
p[:, :, :, :, :, ux, uy] = 0
q[:, :, :, :, :, ux, uy] = 0

print ((p - q) ** 2).sum()  # euclidean distance is 0, so p and q are equal

print ((p ** 2).sum(), (q ** 2).sum())  # prove that not all entries are 0 ;) - This has a negative result due to an overflow