Custom array index, multiple index one value

Hello there,
my question is about the optimal way to realize the following data structure in Julia:

I have a couple of very large rank 4 arrays from a numerical simulation. Typical, these have around 500^4 complex valued elements. However, only about 10% of these element are actually independent of each other, due to symmetries. I have found a way to compute all symmetry classes and compute these objects for one representative of each class.
However, there is no straight forward way to only use these indices in subsequent computations. This is why I though about overloading the getindex functions, which would look up the representative of any given index, using a Dict.
Due to the very large array sizes, this will require a lot of memory operations (for the lookup table) and I was wondering if there is a better option, maybe even a package which is already implemented and tested. I don’t think the structure of these objects is a good fit for sparse arrays.

Best regards

Why can’t you make a custom type that only stores these 10% of values and a getindex function which translates an index through symmetries to an index into these 10%? You can use a sparse matrix to store the subset of values you actually need.

3 Likes

Thank you for your reply. As discussed above, the lookup requires me to store a lookup table. In most cases it would be 5-10GB in size. Since no obviously better approach has been suggested I will try and work around memory limitations using mmap.

Well, you said you computed all symmetry classes. Is there any chance that you can compute the representative from a given index on the fly? Or some kind of invariant over symmetry classes that you could leverage?

Edit: you can also run a hybrid approach if you have a mix of easily computable and complex symmetry relations, so pack the easier ones into a custom getindex and solve the rest with a very reduced lookup table)

I would also consider, what you want to do with this array.

If you would like to use it as an operator then single index access might not be required. Instead you could try using the symmetries and implement its action on a vector. Similar to what you would do when fast Fourier transforming, where the actual transformation matrix is never explicitly created. You could look into

https://github.com/JuliaSmoothOptimizers/LinearOperators.jl

for inspiration.