algorithm - 向量化Pareto Front算法

首先，这是我的设置：
x是包含第一代价函数值的n x 1向量。
y是另一个包含第二个成本函数值的n x 1向量。
a是一个包含要检查的m x 1和x索引的y向量，用于从算法中有选择地排除值。除非需要，否则可以用1:n替换。
可以安全地假设所有组合(x,y)都是唯一的。
任务是找到pareto最优值对集(x,y)，即所有不受支配的值对。如果存在另一对“cc>”，这样的一对称为“支配”，这样(u,v)和其中一个比较是严格的：u <= x && v <= y。换言之，如果另一对在一个值上更好而不是在另一个值上更差，则一对被支配。
到目前为止，我的研究已经产生了三种工作算法，不幸的是它们都依赖于循环。下面是它们的工作原理，以及我用u < x || v < y、x和y长度运行它们的时间：
按升序排序。将第一对添加到pareto集。
循环通过a。将每一对添加到pareto集合，其中1e8低于前一对pareto的x。
运行时间为80.204052秒。
_
找到x。把这对加到pareto集合中。
选择y低于先前添加的对的y的所有对。
转到步骤1，除非步骤2导致空集合。
运行时间为2.993350秒。
_
遍历所有对min(x)。
用y移除所有对y。
运行时间为105.924814秒。
现在我要做的是创建一个矢量化算法。它不必基于上面的任何一个，但我找不到任何其他的工作算法。我能做的就是：

ap = a(y < min(y(x == min(x))) | x < min(x(y == min(y))));

它通常会找到所有pareto最优对，但包括在(x,y)或(u,v)处不受该对支配的所有对，即使其中一对支配另一对。我之所以这么说，通常是因为如果只有一个全局最优对支配着另一对，那么它就完全失败了。用x >= u && y >= v替换min(x)可以解决第二个问题，但可以找到更多的支配对（只有一个更差的值）。我还使用了与上面相同的计时器：
运行时间为0.800385秒。
下面是一个测试脚本，我用它来检查算法的性能，可以随意使用

for i=1:25
    x = randi(8,10,1);
    y = randi(8,10,1);
    a = 1:10;
    ap = a(y < min(y(x == min(x))) | x < min(x(y == min(y)))); %// algorithm here
    figure(1);
    subplot(5,5,i);
    plot(a,x,'b',a,y,'r',ap,x(ap),'b.',ap,y(ap),'r.','MarkerSize',20);
    axis([0,11,0,9]);
    set(gca,'XGrid','on','YGrid','on','XTick',1:10,'YTick',0:8);
    figure(2);
    subplot(5,5,i);
    plot(x,y,'b.',x(ap),y(ap),'ro','MarkerSize',10);
    axis([0,9,0,9]);
end

最佳答案

所以，如果速度是主要特征（在正确之后），那么我发现更快循环版本的递归版本要快30%以上：

>> testPareto(1e8);
Recursive:
Elapsed time is 4.507267 seconds.
Loop:
Elapsed time is 6.136856 seconds.
Vector:
Elapsed time is 7.246806 seconds.

同样，时间取决于机器，甚至可能取决于Matlab的版本。代码如下：

function testPareto(dim)

x = rand(dim, 1);
y = rand(dim, 1);

tic;
rR = paretoRecursive(x, y);
disp('Recursive:');
toc;

tic;
rL = paretoLoop(x, y);
disp('Loop:');
toc;

tic;
rV = paretoVector(x, y);
disp('Vector:');
toc;

end

function result = paretoLoop(x, y)
    result = zeros(numel(x), 2);
    off = 1;
    loop = true;
    while loop
        xmin = min(x);
        ymin = min(y(x == xmin));
        yfilter = y < ymin;
        result(off, :) = [xmin ymin];
        off = off + 1;
        if any(yfilter)
            x = x(yfilter);
            y = y(yfilter);
        else
            loop = false;
            result(off:end, :) = [];
        end
    end
end

function result = paretoRecursive(x, y)
    xmin = min(x);
    ymin = min(y(x == xmin));
    yfilter = y < ymin;
    if any(yfilter)
        result = [xmin ymin; paretoRecursive(x(yfilter), y(yfilter))];
    else
        result = [xmin ymin];
    end
end

function result = paretoVector(x, y)
    xmin = min(x);
    xfilter = x == xmin;
    ymin = min(y(xfilter));
    yfilter = y < ymin;
    if any(yfilter)
        [x, ind] = sort(x(yfilter));
        y = y(yfilter);
        y = y(ind);
        yfilter = [true; y(2:end) < cummin(y(1:end-1))];
        result = [xmin x(yfilter)'; ymin y(yfilter)']';
    else
        result = [xmin ymin];
    end
end