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