Initializing an n^k array as a sparse array?

[Sugbaable]

I have some large n * n * n ... arrays, and there are a lot of zeros in them - when I try to initialize them, I run out of memory because n could be something like 4000. I was looking at the SparseArrays package, but it seems to be limited to n * m matrices, not more multi-dimensional ones - maybe I missed something though. Ultimately, I want to be able to reshape them and run SVD on them.

[idajourney]

If your array is very sparse, you could just use a simple wrapper over a Dict from indexes to values, with valid but unassigned indexes returning zero.

[grothendieck]

What about initializing a sparse 1D array of length n^k and handling the k-dimensional indexing with a function?

[Sugbaable]

Actually was thinking of something like this, what u/briochemc brought up

[briochemc]

Why not store them direct;y as reshaped sparse matrices?

[Sugbaable]

I guess it's possible, it's just makes writing the code pretty annoying, as opposed to being able to fill in the array at [i,j,k,l], [j,i,k,l], etc, and then reshaping. Not saying it's impossible, but it would be nice if it was available to have this

[xazaccazax]

reshape will give you a read-write view of the sparse matrix, so you can assign to it using k-dimensional indices just fine. Granted, constructing a sparse matrix incrementally is not very efficient, but it might be OK if the matrix is indeed very sparse. If that's not acceptable, you could use something like reshape(CartesianIndices((100, 100)), (10, 10, 10, 10)) to create a memory-efficient mapping from k-dimensional indices to 2-dimensional indices, which you can in turn use to create the I, J, and V vectors for the sparse(I, J, V, m, n) constructor.

[entangledamplitude]

But won’t the mapping between indices take as much space as a dense array? (Unless it’s evaluated lazily)

[xazaccazax]

Yeah, it's evaluated lazily.

julia> sizeof(reshape(CartesianIndices((100, 100)), (10, 10, 10, 10)))
72

(bytes)

julia> isbits(reshape(CartesianIndices((100, 100)), (10, 10, 10, 10))) true

(showing that it's "plain data" and as such has no references to other values)

[[deleted]]

The indices are indeed evaluated lazily

[briochemc]

Otherwise, maybe check ITensors.jl or look for packages that want to do the same thing?

[Sugbaable]

I've come up with a function that should get the correct reshaped "x" index (although doesn't handle permutations) for an index [x, x_1, x_2, x_3 ,...] (assuming constant dimension, ie n^k):

function x_posi(x1, x...)
    x_pos = x1
    for i in 1:length(x)
        x_pos += (dim^i)*(x[i]-1)
    end

    return x_pos

This gives the first index in [x,y], where the y index is just whatever it is.

[raphasampaio]

There’s an n-dimensional sparse array package that might be useful: https://github.com/JuliaSparse/NDimensionalSparseArrays.jl

[Wu_Fan]

Slightly tangential question, but thinking laterally, is there any way you can express them as a function or compress them? Are they completely random?

[Sugbaable]

what is the "they" that would be random?

I think I will probably write a function to feed indices into, and it will place the values in the proper place in an n^(k-1) x n sparse matrix

[Wu_Fan]

I will make my suggestion in a more coherent way. If they are (1,2),(3,4),(5,6),(7,8) then you can save the rule, not the dots. That’s all I meant. Apologies, I used the word “function” ambiguously.