Fasten the speed of assigning values to a large tensor

Hi! Now I have a for-loop assigning values to some entries of a large tensor at every iteration. But it is too slow since the tensor is large. Now I am thinking if there is any way to fasten this speed. One simple example is shown below. Here I used 'ndSparse' function to create a four dimensional tensor instead of 'zeros' function. The link of 'ndSparse' is https://www.mathworks.com/matlabcentral/fileexchange/29832-n-dimensional-sparse-arrays. ndSparse generalizes 'sparse' function from sparse matrix to high dimensional sparse tensor.

% Simple example
M = 500;
A = ndSparse.build([M,M,M,M]);
n = 10;

for i=1:n
    disp(i)
    idx1 = M/n*(i-1)+1:M/n*i;
    idx2 = M/n*(i-1)+1:M/n*i;
    idx3 = M/n*(i-1)+1:M/n*i;
    idx4 = M/n*(i-1)+1:M/n*i;
    tmp = randn([M/n,M/n,M/n,M/n]);
    A(idx1,idx2,idx3,idx4) = tmp;
end

Is there any way to make it faster?

Thank you very much!

[Matt J]

What are the typical magnitudes of n, m, and M? Also, when you say you are "creating a sparse all zero 6-D tensor", don't you just mean you are doing this:

A=zeros(M,M,M,M,M,M)

If there is a reason you have not posted it this way, you should tell us.

[Matt J]

Also, are the idx1,idx2,idx3 always assigning into the first 3 dimensions, or is that n-dependent?

[Jason Yang]

@Matt JYes, always. In every iteration, I get different idx1,idx2,idx3.

[Jason Yang]

@Matt J M=7e+2, n=1e+2. m=1e1. The magnitude is not quite accurate since I omit some details.

Yes. I did something just like

A = zeros(M,M,M,M,M,M)

But I create A in a sparse way using 'ndSparse' function provided in the community, which generalizes 'sparse' function from sparse matrix to sparse tensor.

A = ndSparse.build([M,M,M,M,M,M])

Sorry I didn't make it clear!

[Jason Yang]

Hi @Matt J@Jan, I updated my question with a working example. Hope I can make my question clear. Could you give me some help?

Answer 1

[Jan]

Avoid creating index vectors explicitly:

    idx1        = a1:b1;
    A(idx1, :)  = tmp;  % Slow
    A(a1:b1, :) = tmp;  % Faster    

In the first case, each element of idx1 is tested, if it is a valid index:

% For each index something like this is performed:
if idx(1) == flor(idx(1)) && idx(1) > 0 && idx(1) < size(A, 1)

With a1:b1 Matlab checks the validity of the first and last element only.

Please post a minimal working example. Before the readers can test, if their suggestion is faster, they currently have to implement your code based on guesses, what

A = create a sparse all zero 6-D tensor of size (M,M,M,M,M,M);

or

idx1 = a1:b1; (b1-a1=m)

explicitly means.

[Jason Yang]

Thanks for your answer, I tried but it seemed that

A(idx1, :)  = tmp;

is faster. Below is one example. Here I used 'ndSparse' function to create a four dimensional tensor instead of 'zeros' function. The link of 'ndSparse' is https://www.mathworks.com/matlabcentral/fileexchange/29832-n-dimensional-sparse-arrays. ndSparse generalizes 'sparse' function from sparse matrix to high dimensional sparse tensor.

% Simple example
M = 500;
A = ndSparse.build([M,M,M,M]);
n = 10;
tic
for i=1:n
    disp(i)
    idx1 = M/n*(i-1)+1:M/n*i;
    idx2 = M/n*(i-1)+1:M/n*i;
    idx3 = M/n*(i-1)+1:M/n*i;
    idx4 = M/n*(i-1)+1:M/n*i;
    tmp = randn([M/n,M/n,M/n,M/n]);
    A(idx1,idx2,idx3,idx4) = tmp;
end
toc

tic
for i=1:n
    disp(i)
    tmp = randn([M/n,M/n,M/n,M/n]);
    A(M/n*(i-1)+1:M/n*i,M/n*(i-1)+1:M/n*i,M/n*(i-1)+1:M/n*i,M/n*(i-1)+1:M/n*i) = tmp;
end
toc

[Jan]

@Jason Yang: You are right, ndSparse is infact an exception, because this user defined class is not supoorted by Matlab's JIT acceleration.

Answer 2

[Matt J]

You can save some time if you pre-allocate the array:

% Simple example
M = 400;
n = 10;

tic
    A = ndSparse.spalloc([M,M,M,M],(M/n)^4*n); %Pre-allocate
    %A=ndSparse.build([M,M,M,M]);
    for i=1:n
        idx1 = M/n*(i-1)+1:M/n*i;
        idx2 = M/n*(i-1)+1:M/n*i;
        idx3 = M/n*(i-1)+1:M/n*i;
        idx4 = M/n*(i-1)+1:M/n*i;
        tmp = randn([M/n,M/n,M/n,M/n]);
        A(idx1,idx2,idx3,idx4) = tmp;
    end
toc%Elapsed time is 1.282638 seconds.

I don't think there are great gains to be made, however. Most of the computation is being spent just generating the entries for the array:

clear B

tic; 
    [B{1:4}]=ndgrid(1:M/n);
    B = reshape( cat(5,B{:}) ,[],4); 
    
    [C,vals]=deal(cell(n,1));
    
    for i=1:n
        C{i}=B+M/n*(i-1);
        vals{i} = randn(prod([M/n,M/n,M/n,M/n]),1);
    end
    
     C=cell2mat(C); vals=cell2mat(vals);
    
    I=sub2ind([M,M,M,M],C(:,1),C(:,2),C(:,3));
    J=C(:,4);
toc; %Elapsed time is 1.499094 seconds.

tic
 S=sparse(I,J,vals, M^3,M, (M/n)^4*n);
toc; %Elapsed time is 0.344042 seconds.

tic;
 A=ndSparse(S,[M,M,M,M]);
toc;%Elapsed time is 0.016531 seconds.