Optimize filling of sparse matrix: "This sparse indexing expression is likely to be slow."

Hi,

I have a sparse matrix for which I'm replacing certain values in a foor-loop structure. It all goes on the lines of

A = sparse(n,n);
for i = 1:size(Im,1);
  for j = 1:size(Im,2);
      for k = 1:size(Im,3);
        row = ind(i,j,k);
        if(Labels(i,j,k) == -1)
           continue 
        elseif(Labels(i,j,k) == 1) || (Labels(i,j,k) == 2)
            A(row,row) = 1;
            continue
        end

          [sx, sx_value, sy, sy_value, sz, sz_value] = myFunction(Labels,i,j,k);

          A(row,ind(i+sx,j,k)) = A(row,ind(i+sx,j,k)) + sx_value;
          A(row,ind(i,j+sy,k)) = A(row,ind(i,j+sy,k)) + sy_value;
          A(row,ind(i,j,k+sz)) = A(row,ind(i,j,k+sz)) + sz_value;
          end
        end
    end
  end

Basically what I have is a 2D matrix Im and a (large) number n over which I create a sparse matrix A. I then loop over all dimensions (perhaps not optimal, I know) and for certain situations (that is certain rows of n) I call a function, and set values of A with the output of that function, on the lines of

A(i,j) = myValue;

For large values of n (e.g. 500, 1000, etc.) the code is very slow. I ran some profiling and the specific lines where I fill values in A indicates as

This sparse indexing expression is likely to be slow.

On the help I saw some suggestions of not doing this to sparse matrices but rather initiating e.g. A = zeros(...) and 'sparsifying' it in the end. The problem when n is large though, is that I cannot allocate enough memory for A = zeros(...).

Anyone got any suggestions of how to improve the code? Any help would be really great!

Answer 1

[Steven Lord]

In this case I would construct the i, j, and v vectors containing row and column indices and values respectively inside the loop. Call sparse once at the end to turn those vectors into the sparse matrix. See the "Accumulate Values into Sparse Matrix" example on the documentation page for the sparse function for a demonstration of this technique.

[David]

Thanks for the quick response. Just asking to make sure I understand the concept. So let's say I get the indices (i and j) and know what value each combination of indices should have ((i,j) --> v) I can simply assign A as A = sparse(i,j,v,n,n)? I'll try to test it asap and let you know

[Walter Roberson]

Yes, this is a good approach provided you can afford to construct the index vectors.

Another approach that would help instead would be to use spalloc to allocate approximate final number of used entries (better to overestimate a bit). Sparse array can have available slots and filling one of those is done in-place. When there are no available slots then the entire array has to be copied.

[David]

Just finished some tests. Both suggestions does improve the code, but it seems like Steven's initial thought gave the fastest speed-up (for the matrix sizes I am working with). For completeness a code example is below:

I = ones(6*n, 1); J = ones(6*n, 1); V = zeros(6*n, 1); 

iter = 1;
for i = 1:size(Im,1);
  for j = 1:size(Im,2);
      for k = 1:size(Im,3);
        row = ind(i,j,k);
        if(Labels(i,j,k) == exterior)
           continue 
        elseif(Labels(i,j,k) == inlet) || (Labels(i,j,k) == walls)
            I(iter) = row; J(iter) = row; V(iter) = 1; % Set pos. and value
            dirichlet(row) = 1; 
            iter = iter+1;
            continue
        end
        [sx, sx_value, sy, sy_value, sz, sz_value] = get_Ablock_stencil(Labels,i,j,k);
        % Set x-direction pos (I,J) and values (V)
        I(iter:iter+length(sx)-1) = row*ones(length(sx),1); 
        J(iter:iter+length(sx)-1) = ind(i+sx, j, k);
        V(iter:iter+length(sx)-1) = (sx_value./(PixDim(1)^2))';
        iter = iter+length(sx);
        % Set y-direction pos (I,J) and values (V)
        I(iter:iter+length(sy)-1) = row*ones(length(sy),1);
        J(iter:iter+length(sy)-1) = ind(i, j+sy, k)';
        V(iter:iter+length(sy)-1) = (sy_value./(PixDim(2)^2))';
        iter = iter+length(sy);
        if(image_is_3D)
            % Set z-direction pos (I,J) and values (V)
            I(iter:iter+length(sz)-1) = row*ones(length(sz),1);
            J(iter:iter+length(sz)-1) = ind(i, j, k+sz)';
            V(iter:iter+length(sz)-1) = (sz_value./(PixDim(3)^2))';
            iter = iter+length(sz);
        end
      end
  end
end
% Generate sparse matrix at the end
A = sparse(I,J,V);

For n ~ 500-1000 I had a speed-up of almost 80%! Another nice feature (I realized) is that if I,J is ever repeated in the input, sparse simply adds up V (which is exactly what I wanted). For completeness, doing spalloc(n,n,3*n) (or similar) gave ~20% speed-up in this case.

Note (for whoever reads this): initializing with 6*n is an arbitrary choice, based on a-priori info about A.

Again: thanks for all input.

[Monisha S R]

@David Do you use spalloc before using sparse to construct the matrix?