Performance degradation

[gyi]

I've constructed a small script that demonstrates performance degradation in successive calls to solvesos. The degradation is more pronounced on the larger problems I'm working on, but is also present in this test case. Calling clear all
resets execution time, which is impracticle withing a script. Any idea if this a matlab artifact (e.g. delayed garbage collection)?

Thanks for your help.

Matlab version: 8.3.0.532 (R2014a)

Yalmip('version'): 20140619

% variables
x = sdpvar(2,1);
t = sdpvar(1);
P = sdpvar(2);
P(1,1) = polynomial(t, 2);
P(2,1) = polynomial(t, 2);
P(1,2) = polynomial(t, 2);
P(2,2) = polynomial(t, 2);
[s1,s1c] = polynomial([x;t], 2);
[s2,s2c] = polynomial([x;t], 2);
params = [coefficients(P,t); s1c; s2c];
E =[+1.000466767369124-0.002408620734945*t^1-0.184725034833517*t^2+0.005128298848572*t^3, -0.000046676736913-0.374759137926505*t^1+0.018472503483352*t^2-0.000512829884857*t^3
   +0.001555891230420+0.991971264216851*t^1-0.615750116111723*t^2+0.017094329495239*t^3, +0.999844410876958-1.249197126421685*t^1+0.061575011611172*t^2-0.001709432949524*t^3 ];


% solver options
acc = 1e-9;
opts = sdpsettings('verbose',0,'solver','mosek');
opts.mosek.MSK_DPAR_INTPNT_CO_TOL_PFEAS = acc;
opts.mosek.MSK_DPAR_INTPNT_CO_TOL_DFEAS = acc;
opts.mosek.MSK_DPAR_INTPNT_CO_TOL_REL_GAP = acc;


%constraints
T = [1e-3 1];
bounds = (t-T(1))*(T(2)-t);
lyap = x'*(P-eye(2))*x - s1*bounds;
stab = -x'*(E'*P*E - P)*x - s2*bounds;
F = [sos(stab), sos(lyap), sos(s1), sos(s2)];


for k=1:200
    s = tic;
    solvesos(F,[],opts,params);
    p = toc(s); fprintf('k=%d, elapsed=%3.3fs\n',k,p);
end

[Johan Löfberg]

I cannot reproduce this. For how long are you running?

k=1, elapsed=1.282s
k=2, elapsed=0.436s
k=3, elapsed=0.433s
k=4, elapsed=0.434s
k=5, elapsed=0.426s
k=6, elapsed=0.416s
k=7, elapsed=0.422s
k=8, elapsed=0.450s
k=9, elapsed=0.447s
k=10, elapsed=0.432s
k=11, elapsed=0.437s
k=12, elapsed=0.435s
k=13, elapsed=0.437s
k=14, elapsed=0.418s
k=15, elapsed=0.429s
k=16, elapsed=0.425s
k=17, elapsed=0.423s
k=18, elapsed=0.434s
k=19, elapsed=0.434s
k=20, elapsed=0.425s
k=21, elapsed=0.427s
k=22, elapsed=0.444s
k=23, elapsed=0.484s
k=24, elapsed=0.502s
k=25, elapsed=0.424s
k=26, elapsed=0.435s
k=27, elapsed=0.440s
k=28, elapsed=0.430s
k=29, elapsed=0.446s
k=30, elapsed=0.441s
k=31, elapsed=0.430s
k=32, elapsed=0.436s
k=33, elapsed=0.415s
k=34, elapsed=0.433s
k=35, elapsed=0.436s
k=36, elapsed=0.426s
k=37, elapsed=0.421s
k=38, elapsed=0.433s
k=39, elapsed=0.443s
k=40, elapsed=0.422s
k=41, elapsed=0.437s
k=42, elapsed=0.426s
k=43, elapsed=0.434s
k=44, elapsed=0.441s
k=45, elapsed=0.461s
k=46, elapsed=0.512s
k=47, elapsed=0.496s
k=48, elapsed=0.482s
k=49, elapsed=0.434s
k=50, elapsed=0.513s
k=51, elapsed=0.464s
k=52, elapsed=0.431s
k=53, elapsed=0.417s
k=54, elapsed=0.441s

[Johan Löfberg]

OK, saw now you ran till 200. I saw an increase from 0.5s to 0.6s. Is that what you are talking about.

Generally, the most common way to avoid any internal build-up of variables etc is to issue a yalmip('clear') when setting up new models. In your case though, it appears as if you want to run a loop where you set up a complex polynomial model, so then it might be inefficient to model from scratch every iteration. I think I need to see the complete code to be able to give any hints

[gyi]

Yes, that's the increase. It's minor for this test case (which was the smallest example I came up with). The main code (more bloat) exhibits the same behavior but with more rapid buildup. For example,

Cursor: 0.57229, psatz[0.00100, 0.57229] [solv=0.295s] Bad[2]

Cursor: 0.65390, psatz[0.00100, 0.65390] [solv=0.296s] Bad[2]

Cursor: 0.73551, psatz[0.00100, 0.73551] [solv=0.296s] Ok

Cursor: 0.81712, psatz[0.00100, 0.81712] [solv=0.294s] Ok

Cursor: 0.89873, psatz[0.00100, 0.89873] [solv=0.295s] Ok

Cursor: 0.98035, psatz[0.00100, 0.98035] [solv=0.294s] Ok

Cursor: 1.06196, psatz[0.00100, 1.06196] [solv=0.300s] Bad[2]

Cursor: 1.14357, psatz[0.00100, 1.14357] [solv=0.298s] Ok

Cursor: 1.22518, psatz[0.00100, 1.22518] [solv=0.304s] Ok

For this simple problem the results should all be "Ok", but I'm getting sol.problem=2 results. When I graph the solution results (sos residuals and min(eig(sos-grammians)) vs. the p-satz upper bounds), I'm seeing sensitivity of the solver accuracy (between 1e-10 and 1e-8 ... very narrow). To get a better picture, I expanded the problem space, requiring more iterations, leading me to the issue at hand where a 1000 iterations should complete in 5 minutes, but is taking hours. I tried sedumi also, and it has a similar degradation problem (except sedumi is not able to find any solutions where mosek does). I can send you the code if it will help.

It sounds like there are several immediate remedies:

1. Call yalmip('clear') after every N iterations (shouldn't be too much of a burden)

2. Use this test case in place of the main code (for now)

3. Try optimizer() for T(2)? ...though I don't see it being an affine parameterization

Thanks

[Johan Löfberg]

Feels like you are talking about different things here (numerics and setup overhead)

I hope you aren't setting up the model from scratch inside an iteration. Sounds like you only make minor changes, so only that should be inside the loop, if possible. The drawback from using yalmip('clear') in the loop is that it would mean you thus redefine the model from scratch, which typically is inefficient

You are welcome to send the code to see if it can be done differently to avoid a blow-up of internal variables inside YALMIP

[gyi]

Closing out this thread: yalmip('clear') did the trick. I didn't need to clear every iteration, but around every 10-25 iterations.

For my particular problem, I did not have to build the model from scratch. But I would note that if I had a constraints as function of sdpvar, then a clear will wipes them out. In order to save my model between clears, I had to use coefficients() to separate the coefficients of the model from the basis (then only the basis gets wiped out on a clear).

[Johan Löfberg]

What does your code look like? In the most recent versions of YALMIP, I've extended the optimizer framework to support SOS problems with varying data

[gyi]

Here's a re-implementation of the above problem with a) yalmip('clear') and b) matrix of coefficients (mode_coefficients=1). The loop can obviously be made more efficient, but is kept this way to emphasize the effect of yalmip('clear') (or lack thereof)

function [] = SolvesosPerformanceIssue1()
    clc
    yalmip('clear');
    ss = tic;


    % Coefficient-ized polynomial matrix
    % The cost of calculating this matrix is a significant duration of the
    % problem
    E_ = [ 0.005128298849  -0.1847250348   -0.002408620735  1.000466767          
           0.0170943295    -0.6157501161    0.9919712642    0.00155589123          
          -0.0005128298849  0.01847250348  -0.3747591379   -4.667673691e-05      
          -0.00170943295    0.06157501161  -1.249197126     0.9998444109];


    mode_coefficients = 1;


    for k=0:99
        s = tic;


        % Reset
        if ~mod(k,1)
            yalmip('clear');
            p = toc(s); fprintf('[clear] duration=%3.3fs\n',p); s=tic;        
        end


        % Variables

        x = sdpvar(2,1);
        t = sdpvar(1);
        P = sdpvar(2);
        P(1,1) = polynomial(t, 2);
        P(2,1) = polynomial(t, 2);

        P(1,2) = P(2,1);

        P(2,2) = polynomial(t, 2);
        [s1,s1c] = polynomial([x;t], 2);
        [s2,s2c] = polynomial([x;t], 2);
        params = [coefficients(P,t); s1c; s2c];

        if mode_coefficients
            E = polymat(t,E_);
        else

            E =[+1.000466767369124-0.002408620734945*t^1-0.184725034833517*t^2+0.005128298848572*t^3, -0.000046676736913-0.374759137926505*t^1+0.018472503483352*t^2-0.000512829884857*t^3
               +0.001555891230420+0.991971264216851*t^1-0.615750116111723*t^2+0.017094329495239*t^3, +0.999844410876958-1.249197126421685*t^1+0.061575011611172*t^2-0.001709432949524*t^3 ];

        end


        % Constraints
        T = [.1 1];

        bounds = (t-T(1))*(T(2)-t);
        lyap = x'*(P-eye(2))*x - s1*bounds;
        stab = -x'*(E

'*P*E - P)*x - s2*bounds;

        F = [sos(stab), sos(lyap), sos(s1), sos(s2)];

        % Solve
        opts = sdpsettings('verbose',0,'solver','mosek');
        sol = solvesos(F,[],opts,params);
        p = toc(s); fprintf('k=%d, elapsed=%3.3fs\n',k,p);
    end
    fprintf('total elapsed=%3.3fs\n',toc(ss));
   
    function [P] = polymat(var,coefficients);
        d = size(coefficients,2)-1;
        n = sqrt( size(coefficients,1) );
        basis = var.^(d:-1:0)';
        polynomials = coefficients*basis;
        P = reshape(polynomials,n,n);
    end


end

The timing results are below. Clearing for every iteration provides the best performance. And determining the polynomial matrix E from coefficients (mode_coefficients==1) and a univariate basis is more efficient than computing it every time. The benefit of the matrix E_ (of type double) is that it persists outside of yalmip('clear'), while sdpvar E does not.

		mode_coefficients	0	1
			elapsed (s)	elapsed (s)
clear (mod 100)			64.20	62.15
clear (mod 10)			43.71	42.80
clear (mod 1)			41.76	40.70
optimizer			3.13	3.13

A similar problem is constructed using optmizer(); the elapsed time to generate results with this method (code snippet below) is 3.13 s, an order of magnitude improvement over the other methods.

S = optimizer(F,[],opts,[d],params);
p = toc(ss); fprintf('[setup] elapsed=%3.3fs\n',p); s = tic;


% Solve repeatedly
s = tic;
for k=0:99
    sol = S{0};
    p = toc(s); fprintf('k=%d, elapsed=%3.3fs\n',k,p); s = tic;
end
fprintf('total elapsed=%3.3fs\n',toc(ss));

However, I get NaN for sol, which I am currently tracking down.

[Johan Löfberg]

I cannot run the optimizer version (d?) so I cannot help you to debug this

[gyi]

I upgraded to the most recent version of yalmip (20141030), and the NaN issue is resolved.