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.24218 0.088809 0.057279 0.047298 0.056655 0.10572
"40x200" 0.031536 0.28948 0.15877 0.031328 0.033341 0.089172
"60x300" 0.053939 0.41269 0.31192 0.045295 0.051191 0.13935
"80x400" 0.044627 0.84497 0.64379 0.034844 0.039875 0.13325
"100x500" 0.046401 1.3594 1.1704 0.053358 0.10585 0.23048
"120x600" 0.10021 3.0271 2.6949 0.054769 0.11981 0.35816
"140x700" 0.14428 8.9091 7.0499 0.15114 0.16224 0.59638
"160x800" 0.14779 9.6596 8.8541 0.092719 0.22785 0.75742
"180x900" 0.2272 11.017 10.823 0.15018 0.20452 0.77623
"200x1000" 0.13669 8.1434 10.498 0.097073 0.1499 0.67545
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