sparse and gpuArray/sparse errors very odd! (with MWE)

Hi everyone!

When I run the following code, in sections, (MATLAB 2019b), I get the errors indicated in the comments, and it is extremely frustrating, because it is total non-sense!

Sadly, I think these errors are due to flaws in MATLAB.

If anyone experienced these before, and know some workarounds (still involving sparse however), then do please share.

I hope that staff can shed some light on this if noone else can.

%%
sparse(412761000,1,gpuArray(1),6728535000,1) % As 6728535000 > 412761000 this should work
% Error using gpuArray/sparse
% Value out of range.

%% However works on CPU
[i,j,v] = find(sparse(412761000,1,1,6728535000,1)) %<--- on CPU
sparse(i,j,gpuArray(v),6728535000,1) % breaks on GPU again

% Error using gpuArray/sparse
% Value out of range.
% 
% Error in Untitled2 (line 7)
% sparse(i,j,gpuArray(v),6728535000,1)

%% Made my computer go completely nuts and ran for 2 minutes (but succeeded in the end)! (on Quadro P1000 GPU)!:
[i,j,v] = find(sparse(412761000,1,(1),6728535000,1))
sparse(i,j,gpuArray(v))

% It gave the following error (but not if it is run several times in a row):

% Operation terminated by user during
% parallel.internal.shared.buildDisplayHelper>iFirstNNumericDisplayHelper (line
% 60)
% 
% 
% In parallel.internal.shared.buildDisplayHelper>iBuildDisplayHelper (line 33)
%     dh = iFirstNNumericDisplayHelper( ...
% 
% In parallel.internal.shared.buildDisplayHelper (line 24)
%     dh = iBuildDisplayHelper( x, transferDenseFcn, transferSparseFcn,
%     xClassName, xName, N );
% 
% In dispInternal (line 10)
%     dh = parallel.internal.shared.buildDisplayHelper( ...
% 
% In gpuArray/display (line 21)
% dh = dispInternal( obj, thisClassName, objName );

% else the error looks like this:
% Out of memory on device. To view more detail about available memory on the
% GPU, use 'gpuDevice()'. If the problem persists, reset the GPU by calling
% 'gpuDevice(1)'.

More strange stuff:

%% However the following does NOT work on CPU:
sparse(1,412761000,1,1,6728535000)

% Error using sparse
% Requested 1x6728535000 (50.1GB) array exceeds maximum array size preference.
% Creation of arrays greater than this limit may take a long time and cause MATLAB
% to become unresponsive. See array size limit or preference panel for more
% information.
% 
% Error in Untitled2 (line 41)
% sparse(1,412761000,1,1,6728535000)

%% Also does not work with spalloc
S = spalloc(1,6728535000,0) 


% Error using sparse
% Requested 1x6728535000 (50.1GB) array exceeds maximum array size preference.
% Creation of arrays greater than this limit may take a long time and cause MATLAB
% to become unresponsive. See array size limit or preference panel for more
% information.
% 
% Error in spalloc (line 16)
% s = sparse([],[],[],m,n,nzmax);
% 
% Error in Untitled2 (line 53)
% S = spalloc(1,6728535000,0) 


%% Just for fun:
%% The following does work but the "all zero sparse" takes up 3.3 GB memory!!!!:
S = spalloc(1,412761000,0) 
whos S

% S =
% 
%    All zero sparse: 1×412761000
% 
%   Name      Size                        Bytes  Class     Attributes
% 
%   S         1x412761000            3302088024  double    sparse  

Let's make it even more strange!:

%% Both these work: (the subscripts match the indices used before!)
sparse(275174,1500,1,4485690,1500)
sparse(275174,1500,gpuArray(1),4485690,1500)

% ans =
% 
%              (275174,1500)                    1
% 
% 
% ans =
% 
%              (275174,1500)                    1

%% The following is also weird:
%% This does NOT work:
sparse([275174 275174],1500,gpuArray(1),4485690,1500) % <-- On GPU

% Error using gpuArray/sparse
% Vectors must be the same lengths.
% 
% Error in WeirdSparseErrors (line 99)
% sparse([275174 275174],1500,gpuArray(1),4485690,1500)

%% But this does:
sparse([275174 275174],1500,1,4485690,1500) % <--- on CPU

%% However this also does (but is not always a possibility):
sparse([275174 275174],[1500 1500],gpuArray(1),4485690,1500)

Thanks in advance!

//AJ

Answer 1

[Edric Ellis]

@Matt J has already provided some great information here, just adding a few more thoughts on your specific cases.

Case 1:

sparse(412761000,1,gpuArray(1),6728535000,1)

This fails because the required dimensions of the sparse gpuArray exceed intmax('int32'). Unfortunately, the error message you received didn't make this clear.

Case 2:

[i,j,v] = find(sparse(412761000,1,(1),6728535000,1))
sparse(i,j,gpuArray(v))

While the gpuArray gets constructed correctly, displaying it takes a long time. (The display involves a call to gather, and that can be slow because internally sparse gpuArrays are stored in "Compressed Sparse Row" format, whereas on the CPU, MATLAB uses "Compressed Sparse Column").

Case 3:

sparse([275174 275174],1500,gpuArray(1),4485690,1500)

This should work, and the fact that it does not is a bug in the gpuArray implementation of sparse. (The bug is that the scalar expansion of the j parameter is not being supported).

Thanks for reporting these!

[Alexander Jensen]

Thanks a lot Edric, for, once again, great answers!

Matt J did give good information, yes, I was just lacking the last bits of the puzzle, which you filled in ;)

Just one last question. Since gpuArrays use CSR, and CPU uses CSC, would it then make sense to have the larger dimension as the first dimension with sparse gpuArrays, and as the second dimension with CPU?

EDIT:

in regards to:

This should work, and the fact that it does not is a bug in the gpuArray implementation of sparse. (The bug is that the scalar expansion of the j parameter is not being supported).

I hope the "fix", however, does not slow down the creation of sparse matrices on the GPU :)

However, if the scalar expansions happens on the MEX/C side (if it is actually MEX/C - it's only a guess), I would imagine that this would actually be faster, than explicit expansion on the MATLAB side, after which a larger array is passed to the sparse-generator, or am I wrong?

[Edric Ellis]

Yes - on the GPU, sparse rows take up less storage than sparse columns; the other way around on the CPU.

Answer 2

[Matt J]

Ultimately, the reason is that you are creating matrices that violate various internal array size restrictions. There are differences between the size restrictions on the CPU versus the GPU, because they use different libraries, one BLAS the other CUBLAS. See also,

https://www.mathworks.com/matlabcentral/answers/375728-maximum-variable-size-on-gpu-exceeded-on-one-machine-and-not-others

%% The following does work but the "all zero sparse" takes up 3.3 GB memory!!!!:

This happens because a sparse matrix stores additional data besides just the table of matrix entries. Roughly speaking, this additional memory is proportional to the number of matrix columns, so if you have a very long sparse vector, it can be more economical to store it as a column vector than as a row. Example:

>> S=spalloc(1,412761000,0);  St=S.';

>> whos S St

  Name   Size            Kilobytes     Class    Attributes
                                                          
  S      1x412761000       3224696     double   sparse    
  St     412761000x1         0.031     double   sparse    

[Alexander Jensen]

Although the restrictions are different on CPU/GPU:

sparse(275174,1500,gpuArray(1),4485690,1500)

still works, which contains just as many elements as the vectorized one that was not possible to create.

While the restrictions are different due to the different libraries used, the CPU one can atleast be set, while the gpuArray one does not seem to. gpuArray works in mysterious ways!

[Matt J]

Although the restrictions are different on CPU/GPU ... still works, which contains just as many elements as the vectorized one that was not possible to create.

Some of the size restrictions are per dimension, so by reshaping the matrix to have smaller individual dimensions, you fall within the limits. Also, there are memory limits as well. You have far fewer columns in this second example and as I mentioned, columns take up more memory than rows.

[Alexander Jensen]

However, it is kinda silly that it's out of our hands, and that it is per dimension, and not just memory restricted. And besides, I'm not even remotely close to the (GPU) memory limit.

[Matt J]

I would not say that it is “silly”. There has to be an upper limit to matrix dimensions because there is a limit to the range of integers that a computer can represent. I would agree though that it is strange for CUBLAS to use a signed integer limit intmax(‘int32’) for an inherently positive thing like a matrix dimension rather than unsigned integer limit.

[Alexander Jensen]

You might be right, that some limit in terms of the range of representable integers must be in place, but yeah, signed int is kinda strange :)

Thanks a lot for the answers though! It gave some nice insight into things.

[Joss Knight]

cuBLAS just does what the BLAS does, and my understanding is that the BLAS uses signed integers so that it can do integer arithmetic with the sizes. For instance, translating the problem to 8-bit, if you were allowed to specify size up to 255, then you wouldn't be able to store the offset between, say, the size in some dimension for two input arrays in a generic way if the difference were more than 128 - not without checking which is the larger every time. The correct solution is to move to a 64-bit signed integer. cublas's support for that is not complete mainly because GPU memory sizes didn't warrant it - but that of course has changed now.