Two ways of specifying quadratic function, very different performance

[Hendrata Dharmawan]

Hello,

I attached a script below. My problem is a simple quadratic programming, but in the real world scenario I might add some linear constraints and linear objective terms. In the script I was merely isolating the quadratic component of it. 

There are two ways of specifying quadratic function: x' Q x = (Hx)' where Q = H'H. The thing is, depending on how I specify it, the performance can vary greatly, even though mathematically they are supposed to be the same. The script below illustrates the difference. The "second way" is much faster than the "first way." The difference can be observed with n as low as 200, even though for my problem I usually need n around 2000 ish. 

This is not a big deal as I've discovered the difference, and I can structure my code around it. But I was just wondering what the cause of the discrepancy may be, so that in the future when I use YALMIP, I can be more aware of these potential performance hot spots and not make the same mistake again.

Thank you!

Hendrata

-----

clear;

profile clear;

yalmip('clear');

n = 500;

T = 10;

sigma = zeros(T, n, n);

for i = 1:T

    H = 0.1 + 0.03*randn(n, n);

    thisSigma = H * H';

    sigma(i,:,:) = thisSigma;

end

%% First Way

x = sdpvar(n,1);

C = sdpvar(n,n,'full');

% objective: quadratic function

aux = C * x;

objective = aux' * aux;

% constraint

constraints = [];

fprintf('First way: starting to compile problem\n');

tic

options = sdpsettings('solver','cplex','verbose',0,'savesolveroutput',1);

P = optimizer(constraints,objective,options, C,x);

fprintf('First way: time taken to compile problem is : %s\n', hms(toc));

fprintf('First way: starting to solve problem\n');

for i = 1:T

    thisC = reshape(sigma(i,:,:), [n,n]);

    H = cholcov(thisC);

    solution = P(H);

end

fprintf('First way: time taken to solve problem is : %s\n', hms(toc));

%% Second Way

yalmip('clear');

x = sdpvar(n,1);

C = sdpvar(n,n,'full');

% objective: quadratic function

aux = sdpvar(n,1);

objective = aux' * aux;

% constraint

constraints = [aux == C*x];

fprintf('Second way: starting to compile problem\n');

tic

options = sdpsettings('solver','cplex','verbose',0,'savesolveroutput',1);

P = optimizer(constraints,objective,options, C,x);

fprintf('Second way: time taken to compile problem is : %s\n', hms(toc));

fprintf('Second way: starting to solve problem\n');

for i = 1:T

    thisC = reshape(sigma(i,:,:), [n,n]);

    H = cholcov(thisC);

    solution = P(H);

end

fprintf('Second way: time taken to solve problem is : %s\n', hms(toc));

[Johan Löfberg]

Some solvers solve x'*H'*H*x fast, and some prefer to have a sparse objective and have y'*y with constraint y == H*x. THat's just a fact.

However, in YALMIP, x'*H'*H*x will be absolutely horrible when H'*H is dense, for lrge dimensions. That' means YALMIP will have to reason and manipulate with O(n^2) monomial terms. By simpy writing y'*y and a linear constraint y == H*x, the symbolics will be trivial

Setting up an optimizer object with x'*C'*C*x with C a parameter will be beyond horrible for anything but very small problems. This is the object YALMIP has to manipulate for the case n = 3. Now imagine n = 10, or 100...

>> C = sdpvar(3,3,'full');x = sdpvar(3,1);

>> sdisplay(x'*C'*C*x)

C(1,1)^2*x(1)^2+C(2,1)^2*x(1)^2+C(3,1)^2*x(1)^2+2*C(1,1)*C(1,2)*x(1)*x(2)+2*C(2,1)*C(2,2)*x(1)*x(2)+2*C(3,1)*C(3,2)*x(1)*x(2)+2*C(1,1)*C(1,3)*x(1)*x(3)+2*C(2,1)*C(2,3)*x(1)*x(3)+2*C(3,1)*C(3,3)*x(1)*x(3)+C(1,2)^2*x(2)^2+C(2,2)^2*x(2)^2+C(3,2)^2*x(2)^2+2*C(1,2)*C(1,3)*x(2)*x(3)+2*C(2,2)*C(2,3)*x(2)*x(3)+2*C(3,2)*C(3,3)*x(2)*x(3)+C(1,3)^2*x(3)^2+C(2,3)^2*x(3)^2+C(3,3)^2*x(3)^2

[Hendrata Dharmawan]

Hi Johan,

Thank you! In reality, because H is a cholesky decomposition of a real positive definite matrix, it's always triangular. Is there a way to exploit this fact? For example when I said C = sdpvar(n,n), is there a way to tell YALMIP that half of these will be zero?

Your explanation about symbolic manipulation makes perfect sense I just thought that when I define z = 2x + 5y for example, you would create a new sdpvar object and "tack on" that linear constraint as part of "z", but that appears not to be the case. Then the discrepancy of performance makes total sense.

Thank you,

Hendrata

[Johan Löfberg]

Structure/sparsity can be exploited, but in your question it is essentially redundant, as your problem is the quadratic growth w.r.t n, so n^2 vs (n/2)^2 is just noise

https://yalmip.github.io/sparseoptimizer

Moreover, if you have C and quadratic objective x*C'*C*x, you would not have C as the parameter, but simply Q in x'*Q*x, and then you call the optimizer object with the data Cdata'*Cdata. That way, the object is at least linear in the parameter in contrast to quadratic, making it less challenging symbolically (still huge though)

But once again, with these large problems and C / C'*C dense, you will not get it fast