如何在 Matlab 中加快对分位数的调用？

Posted 2023-03-05

【中文标题】如何在 Matlab 中加快对分位数的调用？

【英文标题】：How can I speed up this call to quantile in Matlab?

【发布时间】：2012-01-25 23:14:16

【问题描述】：

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

函数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

有什么办法可以加快速度吗？代码可以向量化吗？

【问题讨论】：

您使用什么分析器来获取该配置文件？ 10x

>> f=@(x)exp(x); profile on; x=f(1); profile viewer

【参考方案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 解决方案稍慢。

【讨论】：

谢谢，这让我比安德烈的建议有了一点加速——现在看的是 14 秒而不是 18 秒。当我有时间时，我会将配置文件结果添加到我的答案中。

【参考方案2】：

我认为你应该使用histc

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

您可以在 matlab 的文档中看到：

说明

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

【讨论】：

我一定是不清楚。我想知道输入中的每个条目都属于哪个存储桶。例如，如果我给出输入 x=[3 1 4 6 7 2] 和 nLevels=2 我希望看到输出 [1 1 2 2 2 1] 因为第一个、第二个和第六个条目落在第一个桶中（即最小的一半，因为我们只有 2 个桶) 和其他人落入第二个桶。

【参考方案3】：

我做了一些改进（包括在另一个答案中受 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 秒，这代表我的代码部分比 MathWorks 快了 150%。

最后，根据 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

以及分析结果：

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

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

【参考方案4】：

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

n = hist(Y,nbins)

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

编辑：在评论之后：

可以使用histc

的第二个输出参数

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

【讨论】：

不，我想知道物品落入了哪个桶。例如，对于输入向量 (3, 1, 4, 6, 7, 2)，如果我使用 nLevels=2 在其上运行 levels 函数，我希望看到输出 (1, 1, 2, 2 , 2, 1)。

太好了，谢谢！我不得不在列上编写一个循环，以便单独处理每一列，但这似乎表现得更好。我会用详细信息更新我的答案。已投票。

【参考方案5】：

这个怎么样

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 的开始。不过我觉得没必要。

【讨论】：

这会将整个矩阵的条目放入存储桶中，但不会单独处理列。但是，你给了我一个加快我的代码速度的想法——如果它有效，我会告诉你的！

@ChrisTaylor 我错过了你的问题中应该对每一列进行操作的部分。修改上面的内容以这种方式运行应该不难。如果您需要更多帮助，请告诉我。

谢谢，我在我自己的答案中使用了您在循环的每个阶段添加逻辑数组的想法（如下）。它给了我一个不错的加速，所以我非常感激！已投票。

不客气。我将建议在循环内使用矩阵 q 的相应第 i 行的 repmat 来进行 >=q 测试。 X>=repmat(q(i,:),size(X,1),1);但是，我对 bsxfun() 不熟悉，它似乎可以自己处理所有这些。 ——