How to improve this code?

Hi,

I need to fill a very sparse matrix Pjj of size n-by-n where n is the total number of pixels in an image. The function receives as inputs, a weights matrix W (n-by-n sparse matrix) which has multiple non-consecutive diagonals (~9 or more, k), a features vector F (~25 or more, sift keypoints linear indices), the value n for the total number of pixels, the value k for the number of diagonals, and the constant beta which is used in the loop.

How would you improve this code?

    function Pjj = matrix_inner(W, F, n, k, beta)
        P = F(:,3);     % keep all the linear indices of keypoints
        D = sum(W, 2);  % pre-calculate the sum for every rows (n)
        T = zeros(k,n); % pre-allocate avalue vector
        L = T;          % pre-allocate an indices vector

        % for every rows of Pjj
        % find the indices for the values on the diagonals on row i
        % create an indices vector to fit the values in L
        % keep the diagonals indices on row i
        % if any feature indice equals to row i, put Z value, or 0
        % keep the Z values for all the diagonals of row i 
        for i = 1:n
            N = find(W(i,:));
            M = 1:length(N);
            L(M,i) = N;
            Z = ((1-(beta*any(P==i)))*W(i,N))/D(i);
            T(M,i) = Z;
        end

        % find the Y (row) and X (col) indices from L to fit T
        Y = repmat((1:n), [k 1]); 
        Y = Y(L>0);
        X = L(L>0);

        % create the n-by-n Pjj sparse matrix with Y, X, T
        Pjj = sparse(Y,X,T(L>0),n,n);
    end

I timed parts of the proceess. It takes about ~0.02-0.04 to get N, and ~0.0005 for the Z. It would take more than one hour to process a 320x480 image... The bottleneck is clearly the n-loop, but I cannot think of another way to do this.

I would appreciate any comments on my code, thank you!

[Jean-Philippe]

Would it be possible to reverse the problem and loop through the features instead?

[Jean-Philippe]

Would it be a possibility to use spdiags to get all the data from W and then try to find the values for the rows in W corresponding to the indices in F?