比较 `coneprog` 算法的速度

打开实时脚本

此示例显示了具有不同问题规模的 coneprog 和 LinearSolver 选项的所有算法的解时间。该问题是找到一个点到一个椭圆体的距离，其中该点位于 n 维度，并且椭圆体由具有 m 行的圆锥约束表示，约束矩阵为。从 1 到 10 中选择 (n,m) = i*(100,20) 作为 i。此示例末尾的 define_problem 辅助函数为 m、n 的指定值以及随机数生成器的种子创建了问题。该函数创建伪随机锥，其中矩阵每行有 10 个 1 的条目，每列有至少两个 1 的条目，并确保矩阵第一列是 1 的（密集）列。

准备问题数据

设置问题生成函数的参数。

n = 100;
m = 20;
seed = 0;

将实验设置为针对十种问题规模进行运行。

numExper = 10;

创建 LinearSolver 选项值的完整列表。

linearSolvers = {'auto','augmented','normal','schur','prodchol','normal-dense'};

对于该数据，'auto' 设置导致 coneprog 使用 'prodchol' 线性求解器，因此对于这两个值您可以获得基本相同的结果。

创建结构来保存结果计时数据和每次运行所需的迭代次数。'normal-dense' 算法需要特殊处理，因为结构名称不能包含连字符。

time = struct();
s = " ";
time.probsize = repmat(s,numExper,1);
% Initialize time struct to zeros.
for solver_i = linearSolvers
    if strcmp(solver_i,'normal-dense')
        time.normaldense = zeros(numExper, 1);
    else
        time.(solver_i{1}) = zeros(numExper, 1);
    end
end

iter = struct();
iter.probsize = repmat(s,numExper,1);
for solver_i = linearSolvers
    if strcmp(solver_i,'normal-dense')
        iter.normaldense = zeros(numExper, 1);
    else
        iter.(solver_i{1}) = zeros(numExper, 1);
    end
end

预热求解器

为了获得有意义的时间比较，请运行 solve（调用 coneprog）几次，但不对结果进行计时。这种“预热”可以让求解器高效地使用数据，并预先填充内部的即时编译器。

[prob, x0] = define_problem(m, n, seed);
options = optimoptions('coneprog','Display','off');
for i = 1 : 4
    sol = solve(prob,x0,'options',options);
end

运行求解器

对所有问题运行求解器，同时记录解时间和求解器所进行的迭代次数。

for i = 1:numExper
    % Generate problems of increasing size.
    [prob, x0] = define_problem(m*i, n*i, seed);
    time.probsize(i) = num2str(m*i)+"x"+num2str(n*i);
    iter.probsize(i) = num2str(m*i)+"x"+num2str(n*i);
    % Solve the generated problem for each algorithm and measure time.
    for solver_i = linearSolvers
        if strcmp(solver_i,'normal-dense')
            options.LinearSolver = 'normal-dense';
            [~,~,~,output] = solve(prob,x0,'options',options);
            time.normaldense(i) = toc;
            iter.normaldense(i) = output.iterations;
        else
            options.LinearSolver = solver_i{1};
            tic
            [~,~,~,output] = solve(prob,x0,'options',options);
            time.(solver_i{1})(i) = toc;
            iter.(solver_i{1})(i) = output.iterations;
        end
    end
end

显示结果

显示计时结果。probsize 列表示问题大小为 "m x n"，其中 m 是锥约束的数量，n 是变量的数量。

timetable = struct2table(time)

timetable=10×7 table
     probsize       auto      augmented     normal      schur      prodchol    normaldense
    __________    ________    _________    ________    ________    ________    ___________

    "20x100"       0.23528     0.10614     0.047865    0.038705    0.038369     0.079893  
    "40x200"      0.037537     0.30514      0.16188    0.033747    0.033241     0.078169  
    "60x300"      0.037295     0.46518      0.31949    0.050156    0.058039      0.11347  
    "80x400"      0.043227     0.89629      0.68154    0.035713    0.041157      0.13298  
    "100x500"     0.047905      1.2003       1.2287    0.041476    0.053212      0.18065  
    "120x600"      0.12594      2.6806       2.7509    0.053601     0.10175       0.3098  
    "140x700"      0.12933      9.2429       7.4331     0.17999     0.17228      0.59548  
    "160x800"      0.13766      9.9432       9.4517    0.099113     0.19963      0.77809  
    "180x900"      0.16249      11.273       11.292     0.22298     0.15892      0.77537  
    "200x1000"     0.23294      7.7597        10.97     0.10126     0.20274      0.73988

最短时间出现在 auto、schur 和 prodchol 列中。auto 和 prodchol 算法对于这些问题是相同的，因此任何计时差异都是由于随机效应造成的。最长的时间出现在 augmented 列，而 normal 列的时间介于两者之间。除了最小的问题外，所有问题的 normaldense 时间都小于 normal 和 augmented 时间。

计时结果的差异是由于每次迭代的速度差异还是由于每个求解器的迭代次数差异？显示相应的迭代计数表。

itertable = struct2table(iter)

itertable=10×7 table
     probsize     auto    augmented    normal    schur    prodchol    normaldense
    __________    ____    _________    ______    _____    ________    ___________

    "20x100"        7         7           7        7          7            7     
    "40x200"       10        10          10       10         10           10     
    "60x300"        7         7           7        7          7            7     
    "80x400"        7         7           7        7          7            7     
    "100x500"       7         7           7        7          7            7     
    "120x600"      18        10          10       10         18           10     
    "140x700"      17        18          18       19         17           19     
    "160x800"      15        15          15       15         15           15     
    "180x900"      12        13          13       12         12           13     
    "200x1000"      9         9           9        9          9            9

对于这个问题，迭代次数与问题大小没有明显的相关性，这是 coneprog 使用的内点算法的典型特征。对于所有算法来说，每一行的迭代次数几乎相同。对于该问题来说，schur 和 prodchol 迭代比其他算法的迭代速度更快。

辅助函数

以下代码创建 define_problem 辅助函数。

function [prob, x0] = define_problem(m,n,seed)
%%% Generate the following optimization problem
%%%
%%% min t
%%% s.t.
%%%     || Ax - b ||   <= gamma
%%%     || x - xbar || <= t
%%%
%%% which finds the closest point of a given ellipsoid (||Ax-b||<= gamma)
%%% to a given input point xbar.
%%%

rng(seed);

% The targeted total number of nonzeros in matrix A is
% 10 nonzeros in each row
% plus 2 nonzeros in each column
% plus a dense first column.
numNonZeros = 10*m + 2*n + m;
A = spalloc(m,n,numNonZeros);

% For each row generate 10 nonzeros.
for i = 1:m
    p = randperm(n,10);
    A(i,p) = 1;
end

% For each column generate 2 nonzeros.
for j = 2:n
    p = randperm(m,2);
    A(p,j) = 1;
end

% The first column is dense.
A(:,1) = 1;

b = ones(m, 1);
gamma = 10;

% Find a point outside of the ellipsoid.
xbar = randi([-10,10],n,1);
while norm(A*xbar - b) <= gamma
    xbar = xbar + randi([-10,10],n,1);
end

% Define the cone problem.
prob = optimproblem('Description', 'Minimize Distance to Ellipsoid');
x = optimvar('x',n);
t = optimvar('t');

prob.Objective = t;
prob.Constraints.soc1 = norm(x - xbar) <= t;
prob.Constraints.soc2 = norm(A*x - b) <= gamma;

x0.x = sparse(n,1);
x0.t = 0;

end

另请参阅

coneprog

比较 coneprog 算法的速度