Search code examples
matlabfor-loopsparse-matrix

Accelerating the index in huge sparse matrix during loop in MATLAB


I need to construct a huge sparse matrix in iterations. The code is as follow:

function Huge_Matrix = Create_Huge_Matrix(len, Weight, Index)

    k = size(Weight,1);

    Huge_Matrix = spalloc(len, len,floor(len*k));


    parfor i = 1:len
        temp = sparse(1,len);
        ind = Index(:,i);
        temp(ind) = Weight(:,i);
        Huge_Matrix(i,:) = temp;
    end


    Huge_Matrix = Huge_Matrix + spdiags(-k*ones(len,1),0,len,len);


end

As is shown, len is size of the height * weight of the input image, for 200*200 image, the len is 40000! And I am assigning the Weight into this huge matrix according the position stored in Index. Even though I use parfor to accerlate the loop, the speed is very slow.

I also try to create full matrix at first, it seems that the code can becomes faster, but memory is limited. Is there any other way to speed up the code? Thanks in advance!


Solution

  • As @CrisLuengo says in the comments, there is probably a better way to do what you're trying to do than to create a 40kx40k matrix, but if you have to create a large sparse matrix, it's better to let MATLAB do it for you.

    The sparse function has a signature that takes lists of rows, columns and the corresponding values for the nonzero elements of the matrix:

    S = sparse(i,j,v) generates a sparse matrix S from the triplets i, j, and v such that S(i(k),j(k)) = v(k). The max(i)-by-max(j) output matrix has space allotted for length(v) nonzero elements. sparse adds together elements in v that have duplicate subscripts in i and j.

    If the inputs i, j, and v are vectors or matrices, they must have the same number of elements. Alternatively, the argument v and/or one of the arguments i or j can be scalars.

    So, we can simply pass Index as the row indices and Weight as the values, so all we need is an array of column indices the same size as Index:

    col_idx = repmat(1:len, k, 1);
    Huge_Matrix = sparse(Index, col_idx, Weight, len, len);
    

    (The last two parameters specify the size of the sparse matrix.)

    The next step is to create another large sparse matrix and add it to the first. That seems kind of wasteful, so why not just add those entries to the existing arrays before creating the matrix?

    Here's the final function:

    function Huge_Matrix = Create_Huge_Matrix(len, Weight, Index)
    
       k = size(Weight,1);
        
       % add diagonal indices/weights to arrays
       % this avoids creating second huge sparse array 
       Index(end+1, :) = [1:len];
       Weight(end+1, :) = -k*ones(1,len);
       
       % create array of column numbers corresponding to each Index
       % make k+1 rows because we've added the diagonal
       col_idx = repmat(1:len, k+1, 1);
    
       % let sparse do the work
       Huge_Matrix = sparse(Index, col_idx, Weight, len, len);
    
    end