`coneprog` アルゴリズムの速度の比較

ライブスクリプトを開く

この例では、さまざまな問題のサイズ、および LinearSolver オプションのすべてのアルゴリズムについて、coneprog の求解時間を示します。この問題は、点から楕円までの距離を明らかにするというものです (ただし点は n 次元で表現され、楕円は制約行列が m 行の錐制約で表現される)。i が 1 ～ 10 の場合、(n,m) = i*(100,20) を選択します。この例の最後に記載されている補助関数 define_problem は、指定された m と n、そして乱数発生器のシードで問題を生成します。この関数は、行列の各行に 1 が 10 個、各列に 1 が 2 個以上ある疑似ランダム錐を生成し、行列の開始列がいくつもの 1 を含む (1 が密集した) 列になるようにします。

問題のデータの準備

問題を生成する関数のパラメーターを設定します。

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

10 通りのサイズの問題に対して実行する実験を設定します。

numExper = 10;

LinearSolver オプションの値をすべて列挙します。

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

このデータの場合、'auto' 設定により、coneprog は 'prodchol' 線形ソルバーを使用することになるため、これら 2 つの値からは実質的に同じ結果が取得されます。

結果として取得できたタイミングデータと各実験の反復回数を保持するための構造体を作成します。構造体の名前にハイフンを含めることはできないため、'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 が呼び出される) を数回実行します。この "ウォームアップ" により、ソルバーは効率的にデータを使用できるようになり、内部の Just-In-Time コンパイラへのデータ入力が自動で実行されます。

[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.21963    0.074885     0.031169    0.033047    0.035503     0.087617  
    "40x200"      0.055387     0.25089     0.090956    0.035106    0.027389      0.06859  
    "60x300"      0.037744     0.37969      0.11412    0.040638    0.057119      0.11951  
    "80x400"      0.035908     0.79509      0.22753     0.02981    0.035143      0.13754  
    "100x500"     0.039585      1.4861      0.48035    0.034093    0.039057      0.20784  
    "120x600"      0.16931      3.0589      0.88939    0.038203    0.078317      0.47463  
    "140x700"      0.11269      9.0631       2.0746     0.17058     0.14926       0.9722  
    "160x800"      0.18219      8.7183       2.5639    0.055328     0.16489      0.67185  
    "180x900"      0.10618      9.9944       3.2133      0.1645     0.11613       0.9407  
    "200x1000"     0.18115      7.8044       2.9427    0.055434     0.10911      0.92372

最短時間は 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 アルゴリズムの速度の比較