Where spfun meets bsxfun

Question

Is there way to perform element-wise operations like bsxfun does, but only operate on the non-zero elements of a sparse matrix?

In particular, for each non-zero element in matrix A at position (i,j) I want to find the product of all of the non-zero elements in the i-th row, except for element (i,j).

For example if the i-th row looks like this:

0   5   3   0   0   4   0  0

The result should look like this:

0   12  20  0   0   15  0  0

The most obvious solution seem to be to take the product of the non-zero elements along each row, and then divide each element back out of the row product. So, in the above example, the row product is 5 x 3 x 4 = 60, and then I just divide out 5 3 and 4 at their respective locations.

Given a sparse matrix A, this is my best solution so far:

[M N] = size(A);
[row col vals] = find(A);
row_prod = accumarray(row,vals,[],@prod);
B = bsxfun(@ldivide, A, row_prod);
B = sparse(row,col,B(sub2ind([M N],row,col)),M,N);

The first three lines achieve what I want: a column vector representing the product of the non-zero elements of each row. However, there are quite a few problems with the last two lines.

bsxfun is going to return a non-sparse matrix the size of A

It is going to waste a LOT of cycles dividing by zero unnecessarily.

The result is a matrix comprised mostly of Inf or -Inf where I really just want zeros.

I can't mask out the Infs because Matlab defines zero times infinity to be NaN.

Do I just need to bite the bullet and write a for-loop for this? Or is there another way to approach it?

Answer 1

I think I have found a solution that resolves most of my concerns above. Sometimes when I have the bsxfun hammer in my hand, the whole world starts to look like a nail. I forget that some of the simple multiplications that I do with bsxfun could be just as easily (and arguably more readably) solved using matrix multiplication. While I don't think that my solution to this problem is any more readable, it is orders of magnitude more efficient than my last solution.

% 'This happens once, outside the loop, since the size'
% 'and sparsity structure of A dont change in the loop'
[M N] = size(A);
[row col] = find(A);

%% 'Inside iterative loop'

    % 'Get the product of non-zero row elements'
    row_prod = accumarray(row,nonzeros(A),[],@prod);

    % 'Use row products to compute 'leave-one-out' row products'
    B = spdiags(row_prod,0,M,M)*spfun(@(x) 1./x, A);

I would still be interested in hearing other suggestions if this can be improved.