About memory saving issue on large scale LCQP problems

[yuxiao]

I am trying to solve a linearly constrained quadratic programming problem:

%---------------------------------------------------------------------------------------------------------------------------------

% For simplification, the known variables are omitted

% variables

e_p = sdpvar(N,n);

e_q = sdpvar(N,n);

e_v = sdpvar(N,n);

e_va = sdpvar(N,n);

C = sdpvar(2*n,1);

e_Y = [e_p e_q];

e_X = [e_va e_v];

e_B_T = sdpvar(n,n);

e_G_T = sdpvar(n,n);

e_Q_d = sdpvar(n,1);

e_P_d = sdpvar(n,1);

e_A = sdpvar(2*n,2*n,'full');

%constriants

Constraints = [];

Constraints = [Constraints,(Y+e_Y)' == A*X' + e_A*X' + A*e_X' + C*ones(1,N)];

Constraints = [Constraints, e_va(:,ref) == zeros(N,1)];

Constraints = [Constraints, -max(abs(data.P_noised))'<=e_P_d<=max(abs(data.P_noised))'];

Constraints = [Constraints, -max(abs(data.Q_noised))'<=e_Q_d<=max(abs(data.Q_noised))'];

Constraints = [Constraints, A_sum(zero.para==1) == 0];

%objective

if (is_sym == 1)

    Objective = enlarge * sum(e_p.*e_p*(w_p_2.^-1) + e_q.*e_q*(w_q_2.^-1) + e_v.*e_v*(w_v_2.^-1) + e_va.*e_va * (w_va_2.^-1));

else

    Objective = enlarge * sum(e_p.*e_p*(w_p_2.^-1) + e_q.*e_q*(w_q_2.^-1) + e_v.*e_v*(w_v_2.^-1) + e_va.*e_va * (w_va_2.^-1))...

    + lamda*sum(sum(abs(w_A.*e_A)));

end

%----------------------------------------------------------------------------------------------------------------------------------------------------------------------------

The scale of the LCQP problem is N=1000 and n=118, but the necessary memory is approximately 26GB.

Is there any methods  save the memory?

Thanks a lot!

[Johan Löfberg]

supply a reproducible example, i.e something that runs and shows the same behaviour

[yuxiao]

Hi Johan！

I have written a demo to show this:

%% This is a demo for LCQP problem

%% set the parameters

N = 1000;

n = 2*118;

noise = 0.01;

lambda = 10;

%% set the real value

rng(1);

X_hat = (rand(N,n)-0.5)*2;

rng(2);

A_hat = (rand(n,n)-0.5)*5;

rng(3);

C_hat = rand(n,1)-0.5;

Y_hat = (A_hat*X_hat' + C_hat)';

%% set the noised data

rng(4,'v5uniform');

X = X_hat + normrnd(0,noise,size(X_hat));

rng(5,'v5uniform');

Y = Y_hat + normrnd(0,noise,size(Y_hat));

rng(6,'v5uniform');

A = A_hat + normrnd(0,noise,size(A_hat));

%% set the variables

e_Y = sdpvar(N,n);

e_X = sdpvar(N,n);

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

C = sdpvar(n,1);

%% set the constriants

Constraints = [];

Constraints = [Constraints,(Y+e_Y)' == A*X' + e_A*X' + A*e_X' + C*ones(1,N)];

%% set the objective function

Objective = sum(sum(e_X.*e_X + e_Y.*e_Y))...

    + lambda*sum(sum(abs(e_A)));

options = sdpsettings('solver','gurobi');

%% Solve the problem

sol = optimize(Constraints,Objective,options);

data.e_X = value(e_X);

data.e_Y = value(e_Y);

[Johan Löfberg]

This is going to be extremely slow with all the high-level modelling and symbolic manipulation required here

To begin with, sum(sum(e.*e)) etc is incredibly inefficient as you will force YALMIP to work with close to half a million symbolic monomial expressions. Same thing with sum(sum(abs()) as you force YALMIP to introduce 236^2 models of abs operators.

You should manually represent the quadratic term using a norm version, or maybe even a low-level cone operator. The abs stuff is more efficiently described as a 1-norm

Still takes time though, takes 1 minute for YALMIP to massage the model before sending to gurobi

[Johan Löfberg]

sdpvar t1 t2

Constraints = [Constraints,norm(e_X(:)) <= t1,norm(e_Y(:)) <= t2];

Objective = t1^2+t2^2 + lambda*norm(e_A,1);

[Johan Löfberg]

and of course, this problem is much likely much much more efficiently solved with a dedicated solver for this kind of problems

[yuxiao]

Thanks for the response!

I have added your code into the initial one, and find that the results get even worse:

Set the parameters N=100, n=10, and lambda=1.

The initial one solves the problem at:

Factor NZ: 1.446e+5

Solving time: 0.30s

But the modified one solves the problem at:

Factor NZ: 1.218e+5

Solving time: 1.51s

Thus the modification makes the problem slower without saving any memories. I wonder if there is any wrong with this LCQP problem?

By the way, is there any dedicated solver recommended?

Thanks a lot!

[Johan Löfberg]

There is absolutely nothing wrong with your model. It is just a fact that you cannot use YALMIP to setup quadratic models with half a million symbolic terms as the overhead grows too big.

Hence, your options are

1. Accept that it takes longer to solve when in a form where massive symbolic quadratics are avoided

2. Test another solver, such as mosek which is very good on SOCP models

3. Interface the solver manually to avoid the symbolic overhead on quadratics

4.  Don't solve this using a standard solver. It appears to be the archetypal case of something you would develop a specialized routine for

[yuxiao]

Another question:

Is the formulation of norm(e_X(:)) equivalent to sum(sum(e_X.*e_X )) in yalmip?

Because e_X is a matrix, the matrix norm is the maximum value of each column norm e_X(:,i), not the sum of it.

So should I use the Frobenius norm instead?

The same question with 1-norm representation 'lambda*norm(e_A,1)'.

Thanks a lot!

[Johan Löfberg]

No, norm(e_X(:))^2 is equal to sum(sum(e_X.^e_X)), that why t^2 is is minimized

The L1-penalty should be norm(e_A(:),1) to be consistent with your original code

and here's a hack which exploits some highly optimized code inside YALMIP to deal with very simple quadratic functions

e = sdpvar(N*n*2,1);

e_Y = reshape(e(1:N*n),N,n);

e_X = reshape(e(N*n+1:end),N,n);

...

Objective = sum(e.*e) + lambda*norm(e_A(:),1);
...

with that, gurobi will be the bottle-neck

[yuxiao]

Thanks for the response, that really helps!

And I will try some low-level cone representations or try other solvers for large-scale improvements.

[Johan Löfberg]

 Cone operator won't help you. You already come to the conclusion that the SOCP representation was slow when using gurobi. THe YALMIP over-head in norm vs cone for your large problems is negligable compared to the time the solver time

You should probably solve this using a specialiserad solver with some least-squares+thresholding/proximal/admm stuff or what ever it is called

[yuxiao]

I ‘ve tried mosek and it really takes less time! And cone representation really makes it faster in large scale problems!

And I still have a question: what is the difference between norm(A) and norm(A(:))?

Recall that you suggest me to add the code:

sdpvar t1 t2

Constraints = [Constraints,norm(e_X(:)) <= t1,norm(e_Y(:)) <= t2];

Objective = t1^2+t2^2 + lambda*norm(e_A,1);

And both norm(A) and norm(A(:)) are used.

[Johan Löfberg]

norm(A) is the frobenius norm.  norm(A(:)) is norm(A,2)

[Johan Löfberg]

meant the other way around

>> norm(A)

ans =

    2.2857

>> norm(A(:))

ans =

    2.3624

>> norm(A,2)

ans =

    2.2857

>> norm(A,'fro')

ans =

    2.3624