8

我有一个带有一个相当明显瓶颈的 MATLAB 例程。我已经对函数进行了概要分析，结果在函数中使用了 2/3 的计算时间levels：

在此处输入图像描述

该函数levels接受一个浮点矩阵并将每一列拆分为nLevels桶，返回一个与输入大小相同的矩阵，每个条目都替换为其所属桶的编号。

为此，我使用quantile函数来获取存储桶限制，并使用循环将条目分配给存储桶。这是我的实现：

function [Y q] = levels(X,nLevels)
% "Assign each of the elements of X to an integer-valued level"

p = linspace(0, 1.0, nLevels+1);

q = quantile(X,p);
if isvector(q)
    q=transpose(q);
end

Y = zeros(size(X));

for i = 1:nLevels
    % "The variables g and l indicate the entries that are respectively greater than
    % or less than the relevant bucket limits. The line Y(g & l) = i is assigning the
    % value i to any element that falls in this bucket."
    if i ~= nLevels % "The default; doesnt include upper bound"
        g = bsxfun(@ge,X,q(i,:));
        l = bsxfun(@lt,X,q(i+1,:));
    else            % "For the final level we include the upper bound"
        g = bsxfun(@ge,X,q(i,:));
        l = bsxfun(@le,X,q(i+1,:));
    end
    Y(g & l) = i;
end

我能做些什么来加快速度吗？代码可以向量化吗？

4

5 回答 5

4

如果我理解正确，您想知道每个桶中有多少物品。采用：

n = hist(Y,nbins)

虽然我不确定它是否有助于加速。这样更干净。

编辑：在评论之后：

您可以使用histc的第二个输出参数

[n,bin] = histc(...) 还返回一个索引矩阵 bin。如果 x 是向量，则 n(k) = >sum(bin==k)。对于超出范围的值，bin 为零。如果 x 是 M×N 矩阵，则

于 2011-12-22T09:50:15.510 回答
2

我认为你应该使用histc

[~,Y] = histc(X,q)

正如您在 matlab 的文档中看到的：

描述

n = histc(x,edges) 计算向量 x 中落在边向量中元素之间的值的数量（必须包含单调非递减值）。n 是包含这些计数的长度（边）向量。x 的任何元素都不能是复数。

于 2011-12-22T09:51:58.230 回答
2

这个怎么样

function [Y q] = levels(X,nLevels)

p = linspace(0, 1.0, nLevels+1);
q = quantile(X,p); 
Y = zeros(size(X));
for i = 1:numel(q)-1    
    Y = Y+ X>=q(i);
end

这导致以下结果：

>>X = [3 1 4 6 7 2];
>>[Y, q] = levels(X,2)

Y =

     1  1  2  2  2  1

q =

     1  3.5  7

您还可以修改逻辑线以确保值小于下一个 bin 的开始。但是，我认为没有必要。

于 2011-12-22T13:33:06.363 回答
1

我进行了一些改进（包括在另一个答案中受 Aero Engy 启发的改进），这些改进导致了一些改进。为了测试它们，我创建了一个包含一百万行和 100 列的随机矩阵来运行改进的函数：

>> x = randn(1000000,100);

首先，我运行未修改的代码，结果如下：

在此处输入图像描述

请注意，在 40 秒中，其中大约 14 秒用于计算分位数 - 我不能指望改进这部分例程（我假设 Mathworks 已经对其进行了优化，尽管我想假设这样做会...... )

接下来，我将例程修改为以下，应该会更快，并且行数更少！

function [Y q] = levels(X,nLevels)

p = linspace(0, 1.0, nLevels+1);
q = quantile(X,p);
if isvector(q), q = transpose(q); end

Y = ones(size(X));

for i = 2:nLevels
    Y = Y + bsxfun(@ge,X,q(i,:));
end

使用此代码的分析结果是：

在此处输入图像描述

所以它快了 15 秒，这代表了我的代码部分的 150% 加速，而不是 MathWorks。

最后，根据 Andrey 的建议（再次在另一个答案中），我修改了代码以使用histc函数的第二个输出，它将条目分配给 bin。它不会独立处理列，所以我不得不手动遍历列，但它似乎表现得非常好。这是代码：

function [Y q] = levels(X,nLevels)

p = linspace(0,1,nLevels+1);

q = quantile(X,p);
if isvector(q), q = transpose(q); end
q(end,:) = 2 * q(end,:);

Y = zeros(size(X));

for k = 1:size(X,2)
    [junk Y(:,k)] = histc(X(:,k),q(:,k));
end

分析结果：

在此处输入图像描述

我们现在在函数之外的代码中只花费了 4.3 秒quantile，这比我最初编写的代码快了大约 500%。我花了一些时间写这个答案，因为我认为它变成了一个很好的例子，说明如何结合使用 MATLAB 分析器和 StackExchange 从代码中获得更好的性能。

我对这个结果很满意，当然我会继续很高兴听到其他答案。在这个阶段，主要的性能提升将来自于提高当前调用的部分代码的性能quantile。我无法立即看到如何执行此操作，但也许这里的其他人可以。再次感谢！

于 2011-12-22T14:08:25.263 回答
1

您可以sort列并除以+舍入反向索引：

function Y = levels(X,nLevels)
% "Assign each of the elements of X to an integer-valued level"
[S,IX]=sort(X);
[grid1,grid2]=ndgrid(1:size(IX,1),1:size(IX,2));
invIX=zeros(size(X));
invIX(sub2ind(size(X),IX(:),grid2(:)))=grid1;
Y=ceil(invIX/size(X,1)*nLevels);

或者您可以使用tiedrank：

function Y = levels(X,nLevels)
% "Assign each of the elements of X to an integer-valued level"
R=tiedrank(X);
Y=ceil(R/size(X,1)*nLevels);

quantile令人惊讶的是，这两种解决方案都比+解决方案稍慢histc。

于 2011-12-22T14:18:42.147 回答