Sparse Array Compression and Decompression upgrade to v0.7+

General Usage
, ,
1

Summary: I can “compress” and “decompress” a sparse array using deprecated functions. I want to upgrade the code to v0.7+. (I want to use this compression/decompression to numerically evaluate the Jacobian of a function of a vector of size ~106.)

  1. I have a function f(x), and a boolean sparse matrix AM that represents the Jacobian’s sparsity pattern. AM is a sort of an adjacency matrix (thus the name AM) of dependencies between different indices of x. AM has non-zero sub indices i and j, and non-zero linear indices ij:

    i, j, _ = findnz(AM)
ij = find(AM)

  2. I also have a home-made vector of colors col_idx which groups columns of the matrix AM by color. (The vector col_idx has integer values.) Two columns (say, j1 and j2) that share the same color (i.e., col_idx[j1] == col_idx[j2]) do not have a common non-zero. (I think the exact term is structural orthogonality.) This allows me to construct the corresponding indices of my compressed objects:

    icol = i
jcol = col_idx[j]
ijcol = sub2ind(icol, jcol)

  3. Once these 6 linear and sub indices (i, j, ij, icol, jcol, and ijcol) are created, I can easily

    • compress using

      compress(Jac, icol, jcol, ij) = sparse(icol, jcol, Jac[ij])

    • and decompress using

      decompress(smallJac, i, j, ijcol) = sparse(i, j, smallJac[ijcol])

Now, with Julia v0.7+, the functions findnz, find, and sub2ind are deprecated, so I am trying to upgrade this code to v0.7+. Deprecation of those function suggest using the CartesianIndex and LinearIndex types, but those are not straightforward to use with sparse matrices. I have been scratching my head for almost a week on related topics now, so any help is much appreciated!

2

The findnz function is not deprecated but part of the SparseArrays stdlib, so I think this becomes:

AM = sparse([1,2,3],[1,2,3],ones(3))
i, j, _ = findnz(AM)
LI = LinearIndices(AM)
ij = LI[i,j]
icol = i
jcol = col_idx[j]
ijcol = LI[icol,jcol]
1 Like
3

Thanks for your reply! However, I think

is not the correct answer because LI[i,j] is a 3x3 array while I want a 3x1 array like ind2sub(i,j).

Instead, I think

ij = LI[CartesianIndex.(i, j)]

works. This ij can then also be used for getting at the values of a sparse matrix S via S[ij].

Thanks for your reply though as it helped be find the solution! I’ll report back if this works all the way through :slight_smile:

More details: I think the issue stemmed from the suggestion in ind2sub(i,j)'s deprecation warning. ind2sub(i,j) may be wrongly interpreted because the suggestion:

julia> ij = sub2ind(i,j)
┌ Warning: `sub2ind(A::AbstractArray, I...)` is deprecated, use `(LinearIndices(A))[I...]` instead.
│   caller = top-level scope at none:0
└ @ Core none:0

is what you suggested (i.e., ij = LI[i,j]) and does not correspond.
Instead, it should have been the one that comes when using the size of the matrix in the arguments:

julia> ij = sub2ind(size(AM), i, j)
┌ Warning: `sub2ind(inds::Union{DimsInteger, Indices}, I1::AbstractVector{T}, I::AbstractVector{T}...) where T <: Integer` is deprecated, use `(LinearIndices(inds))[CartesianIndex.(I1, I...)]` instead.
│   caller = top-level scope at none:0
└ @ Core none:0
2 Likes
4

So after a bit of work, I finally found a decent solution, posting if anyone is interested… Although the solution above with ij = LI[CartesianIndex.(i, j)] worked, I thought I should note that once the compressed sparse matrix smallJac is built, accessing it via smallJac[ijcol] can be quite slow for large sizes. This is because CSC sparse matrices are stored by columns, so that the order of indices are perturbed by compression and decompression. It is thus faster (better?) to figure out the permutation of indices that occurs during compression and use it, via

_, _, perm = sparse(icol, jcol, 1:length(icol))
iperm, jperm = i[perm], j[perm]
decompress(smallJac, iperm, jperm) = sparse(iperm, jperm, smallJac.nzval)

This is faster because accessing the non-zero values of smallJac via its .nzval field is fastest, right? (Jac.nzval can also be used for the compression part, which might be faster and avoids linear indices altogether.)

5

I don’t know if it’s applicable for your use case, but maybe you could also write an AbstractArray implementation that wraps your compressed matrix, eliminating the need for compression/decompression altogether and doing the index translation on the fly?

6

That would be nice but I would need some help! Could you give me some hints on how to do that?

7

Basically you just need to define a new type that holds your compressed array, make it inherit from AbstractArray and then implement the 4 functions described here:

https://docs.julialang.org/en/stable/manual/interfaces/#man-interface-array-1

1 Like