I have an algorithm I'm trying to implement more efficiently. These are the basic steps:
Let X be an n-by-p matrix. Let vector Xi be the i-th column of X, and vector Xj be the j-th column of X. For every combination of Xi and Xj, I want to test whether each element of Xi is larger than the corresponding element in Xj. Because Xi > Xj is different from Xj > Xi, I am interested in all permutations.
I am working in Matlab, so my current approach is to loop over the p columns of X. On each loop, I select the vector Xi, repeat it p times to form an n-by-p matrix (lets call it XI), and then do XI > X. The most efficient way I can think to write this is:
bsxfun(@gt, X(:,i), X);
I have profiled the code, and this operation is far and away the slowest part.
I know that if I could re-frame this operation as matrix arithmetic, it would be much faster. But if that is possible, I do not know how.
Is it possible do this operation faster? If not this operation specifically, maybe a bit more context. What I really want to do is identify whether elements in Xi are exceptionally large relative to the rest of the sample. The next step is to count the number of true values in the logical matrix resulting from (XI > X). The following are the lines of code of chief interest, written on 3 lines to make profiling easier. X in this case is 37,000 x 10,000
x = bsxfun(@gt, X(:,i), X); % 19.748 s (100 calls)
y = sum(x ,2); % 9.819 s (100 calls)
b = y >= rth; % 0.007 s (100 calls)
I thought it might make a difference to transpose X at the outset, since summing over columns rather than rows might be faster. But expanding over rows instead of columns might be slower… here are the numbers with transposed X. No real difference:
x = bsxfun(@gt, X(i,:), X); % 19.353 s (100 calls)
y = sum(x ,1); % 11.549 s (100 calls)
b = y >= rth; % 0.007 s (100 calls)
It may be that this operation just takes time … but I have the sense that it should be possible to do this faster. I am not interested in writing my own C-extension, but if there is a way to do with better within the current Matlab toolkit I would love to hear how.
Best Answer