Sparse gpuArray indexing limitations

When trying to extract a sub-matrix of a sparse gpuArray, I get,

>> S=gpuArray(sprand(100,100,0.1));
>> S(1:5,1:5)

Sparse gpuArray matrices only support referencing whole rows or whole columns.

OK, but here, I am indexing a block of complete rows and get the same error,

>> S(1:5,:)

Sparse gpuArray matrices only support referencing whole rows or whole columns.

Aside from this, what is the internal hurdle that makes subsrefing a sparse matrix so hard on the GPU? This has been a limitation for quite a while now.

[Walter Roberson]

I do not have easy access to a gpu system at this time.

I wonder if

S(1,:)
S(:,1)

would work? That is, perhaps the limit is on indexing more than one complete row or column ?

[Matt J]

Yes, single row/column extraction works,

>> S=gpuArray.speye(3);
>> S(1,:)

ans =

1×3 sparse gpuArray double row vector (1 nonzero)

(1,1) 1

>> S(:,1)

ans =

3×1 sparse gpuArray double column vector (1 nonzero)

(1,1) 1

But I wonder why that would be the limit. If you can do one, surely you should be able to do more?

Answer 1

[dpb]

Technically it appears to be a bug based on the doc at gpuArray doesn't have any explicit restrictions on addressing by either full rows or columns, agreed. However, I suspect this has uncovered limitations on actual useage owing to the unique nature of how the GPU implements sparse storage. Details are lacking in the MATLAB doc's, but I venture it's owing to how the GPU stores the sparse matrix; it uses Compressed Sparse Column (CSC) or Compressed Sparse Row (CSR) layouts and since zeros are not stored, storage is represented by three arrays:

1. Values: The actual non-zero numbers.
2. Inner Indices: The exact row indices (for CSC) or column indices (for CSR).
3. Outer Pointers: A map specifying the physical offset index where each new column or row starts in the arrays above.

Because different rows/columns have a varying number of non-zero elements, the distance between elements is completely dynamic. there is no fixed "stride" to allow for direct addressing without the lookup.

Consequently, trying to pull a sub-matrix (e.g., A(2:5, 2:5)), requires thousands of parallel threads to fetch that data simultaneously. To get to the initial element A(2, 2) a GPU thread must:

1. Look up the Outer Pointer for column 2.
2. Read the start pointer for column 3.
3. Linear-search through every non-zero index inside that block just to see if row 2 exists.

Because the data are tightly packed, a thread handling row 3 cannot know where its data begins until the threads handling rows 0, 1, and 2 have resolved how many elements they contain. This creates a strictly sequential data dependency and forcing a massively parallel architecture like a GPU to execute a sequential pointer-chasing search destroys performance, making arbitrary sub-array slicing structurally infeasible.

The reason MATLAB allows full rows or full columns (depending on if it is utilizing CSR or CSC under the hood) is because the starting and ending bounds of a full slice are explicitly defined by the Outer Pointers array. The GPU can read the single pointer for the requested column and instantly copy that entire contiguous chunk of memory into a new vector.

I don't know if the MATLAB implementation documents whether CSR or CSC is used or if it may depend upon the structure of the sparsity, perhaps. I could see that when creating a particular sparse GPU array, it must choose one or the other and it then could run into the indexing problem if one asks for the opposite slice direction.

I wonder if the example might be an illustration of that; if you try the other direction

S(:,1:5)

on the same array it would succeed?

In short, GPUs are good if you can load 'em up and then they can operate on the data in parallel; anything else is likely to be nothing but a bottleneck that totally defeats the purpose.

[Matt J]

It makes me wonder then, how sparse matrix multiplication A*B is done if submatrices are so difficult to extract. Seemingly, each thread multiplies a whole row of A by a whole column of B?

[dpb]

It uses the NVIDIA libraries which have specialized code designed specifically to use the capabilities built into the NVIDIA hardware and algorithms that calculate partial results using outer products and merge them in shared memory. The architecture also uses internal structured formats which condense storage and that align smaller sparse sections with the GPU's dense array storage tensor core hardware. So, as I understand it, the sparse matrices eventually get mapped indirectly from sparse to dense(+) in a zillion little blocks. Then parallelism gives the performance.

<A paper that may provide more insight> than the summaries I've pulled from and edited down significantly from NVIDIA developer forum and just asking google for background info.

It's a real feat of engineering, indeed, inside there.

In going back to the original question and trying to confirm if the previous surmise was at all on track, does the full row index form work on the given array for you? If that did, I think it would confirm the particular internal sparse storage used; I've yet to find anything that indicates it it always uses one or the other inside, though.

If it worked that way and not the other, it would also be reason that the documentation should be updated to reflect that slices only work in one direction based on the internal storage orientation and perhaps restricted to only one or the other universally if there is always a preferred storage orientation. The error handling might also need some revisting.

(+) I presume but did not find it stated explicitly without more extensive digging that only subarrays with a nonzero element are mapped, not reconstructing the entire array since the others wouldn't add anything to the outer product being calculated.