Efficient pairwise logical operation on all combinations of vectors?

3 次查看（过去 30 天）
Chris
Chris 2018-5-28
评论： Chris 2018-5-28
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.
  2 个评论
John D'Errico
John D'Errico 2018-5-28
Note that this is not valid MATLAB syntax:
x = bsxfun(@gt, X(i,:), X), 2); % 19.353 s (100 calls)
You have a spare set of closing parens there, and an extra ,2, which suggests this is a copy paste error.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2018-5-28
Note that a big part of the problem may be that you are just doing a LOT of tests, and grabbing pretty large chunks of RAM to do them.
That bsxfun call you do will create a matrix of size 37000x10000 for each i. Since it is logical, that is not too bad, but still each element is stored as one byte, thus roughly .37 GB to allocate, etc.
You do not say what version of MATLAb you are using. Current releases of MATLAB allow expansion without BSXFUN, so I had to check the time on that instead.
X = rand(37000,10000);
timeit(@() bsxfun(@gt, X(:,1), X))
ans =
0.51082
timeit(@() X(:,1) > X)
ans =
0.59996
That in R2017b. (Rats. I thought I bad gotten around to doing the upgrade. Or, maybe I have, and I am just running the wrong release. Whatever.)
Anyway, big computations can take big time. We get spoiled by fast computers (and yours seems about 2.5 times faster than mine), expecting everything to happen immediately.
You ask about whether this is doable using linear algebra. That depends on your definitions of what is big. So, yes, you MAY be able to do this more quickly, if you are willing to be a little loose on how you do the computations.
[Xdp,Xtags] = sort(ones(1,37000)*X,'descend');
See that this will be IMMENSELY faster than what you are doing, but it uses a different way of looking at which columns are big. So the 10 largest columns, by one measure are:
Xtags(1:10)
ans =
9493 4864 9289 54 8955 1644 3724 4277 8419 608
The measure I chose was simply which columns have the largest dot product with a vector of ones, which I could have written using sum, but writing it as a dot product makes the idea more explicit.
Xdp(1:10)
ans =
18713 18694 18693 18689 18687 18683 18681 18676 18676 18675
You do say that size of the elements is important. (Hey, sometimes size does matter.) Are the largest elements important? If so, then this might be a better measure:
[Xdp,Xtags] = sort(sum(X.^2),'descend');
Really, it depends on how you define size in your context. Which of these two vectors are larger?
x1 = [1 1 10]
x2 = [2 2 2]
So in effect, you MIGHT be able to choose a column p-norm to do your work. If you really do want to count the number of elements that are larger, then using a small value of p will give you a decent measure. Thus:
p = 0.1;
[Xdp,Xtags] = sort(sum(X.^p),'descend');
Xtags(1:10)
ans =
4864 6561 8031 1563 2228 7500 769 9493 2595 9689
So by this measure, column 4864 is large compared to the other columns.
I'm not sure what your real need is here, so it is difficult to know if your problem can be restructured for more speed. Or if it is really important to do exactly what you are doing now, as you are doing it.
  1 个评论
Chris
Chris 2018-5-28
Thank you so much for this thoughtful reply. I seeing your range of interesting suggestions, I should have been more clear about my goal. The objective is to determine whether specific elements in Xi are larger than typically seen in the sample (all p observations of that element, across that row in X). I'm not interested in the size of the vector for this analysis. But I'll keep these tips in mind, because you've provided some solutions/idioms I would not have come up with myself.
So yeah, my problem is large, and so it takes some time. :) That makes sense. I just implemented a very efficient cross-product analysis last week where Matlab solves (Y' * X) / n, where X is something like 1000 x 37000 and Y was 1000 x 100, in a fraction of a second. I was optimistic that this problem could also be formulated in such a way that Matlab could crush it like that.
Thanks again.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by