Compare Speeds of coneprog Algorithms

This example shows the solution times for coneprog with various problem sizes and all algorithms of the LinearSolver option. The problem is to find the distance of a point to an ellipsoid where the point is in n dimensions and the ellipsoid is represented by a cone constraint with m rows for the constraint matrix. Choose (n,m) = i*(100,20) for i from 1 to 10. The define_problem helper function at the end of this example creates the problem for specified values of m, n, and the seed for the random number generator. The function creates pseudorandom cones with 10 entries of 1 in each matrix row and at least two entries of 1 in each column, and ensures that the first matrix column is a (dense) column of 1s.

Prepare Problem Data

Set the parameters for the problem generation function.

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

Set the experiment to run for ten problem sizes.

numExper = 10;

Create the complete list of LinearSolver option values.

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

For this data, the 'auto' setting causes coneprog to use the 'prodchol' linear solver, so you obtain essentially the same results for those two values.

Create structures to hold the resulting timing data and the number of iterations each run takes. The 'normal-dense' algorithm requires special handling because a structure name cannot contain a hyphen.

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

Warm Up Solver

To obtain meaningful timing comparisons, run solve (which calls coneprog) a few times without timing the results. This "warm-up" prepares the solver to use data efficiently, and prepopulates the internal just-in-time compiler.

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

Run Solver

Run the solver on all the problems while recording the solution times and the number of iterations the solver takes.

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

Display Results

Display the timing results. The probsize column indicates the problem size as "m x n", where m is the number of cone constraints and n is the number of variables.

timetable = struct2table(time)

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

    "20x100"       0.25884    0.083297     0.043986    0.032592    0.036523     0.077243  
    "40x200"      0.032718     0.33583      0.15218    0.031943    0.032378      0.08593  
    "60x300"      0.067792     0.47739      0.30618     0.03186    0.034864       0.1286  
    "80x400"      0.042137     0.92409      0.73938    0.049023    0.048405      0.14189  
    "100x500"     0.051248      1.6447       1.2491    0.051082    0.064689       0.2011  
    "120x600"      0.11088       3.338       2.7731    0.053814     0.11316      0.36706  
    "140x700"      0.14294        9.26       7.3395     0.12978     0.15002       0.6798  
    "160x800"      0.20539       10.39       11.146    0.093721     0.21549      0.73684  
    "180x900"      0.19003      12.692       13.101     0.17081     0.25221       1.2588  
    "200x1000"     0.14401      9.7961       11.964     0.10683     0.15439      0.76893  

The shortest times appear in the auto, schur, and prodchol columns. The auto and prodchol algorithms are identical for the problems, so any timing differences are due to random effects. The longest times appear in the augmented column, while the normal column times are intermediate. The normaldense times are smaller than the normal and augmented times for all but the smallest problems.

Are the differences in the timing results due to differences in speed for each iteration or due to the number of iterations for each solver? Display the corresponding table of iteration counts.

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     

For this problem, the number of iterations is not clearly correlated to the problem size, which is a typical characteristic of the interior-point algorithm used by coneprog. The number of iterations is nearly the same in each row for all algorithms. The schur and prodchol iterations are faster for this problem than those of the other algorithms.

Helper Function

The following code creates the define_problem helper function.

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