Model $\min \frac{1}{2} \parallel Ax-B \parallel_2 + \lambda_1 \parallel Cx \parallel_1 + \lambda_2 \parallel Dx \parallel_\infty $ into standard form

Question

I need to solve the following convex optimization problem:

$\min \frac{1}{2} \parallel Ax-B \parallel_2 + \lambda_1 \parallel Cx \parallel_1 + \lambda_2 \parallel Dx \parallel_\infty$ s.t $x > 0$

$x$ is a vector, varible. A, B, C and D are parameters.

I know there are several modeling tools like CVX and YALMIP that can automatically transfer this problem into standard form that can be solved by solvers like Gurobi and Mosek.

But the CVX and YALMIP are matlab based tools, and they are very slow.

So I want to use the C++ version of Gurobi or Mosek to solve the problem. But there is no C++ version of CVX or YALMIP. That means I have to model the problem by myself. But I don't even know where to begin, can anyone give me some clues?

Thank you!

Answer 1

I'm following up this slightly on Michaels comment that you perhaps should look at the numbers before you rule out anything or start spending days on learning and coding in some other language.

To begin with, on my machine, Gurobi takes roughly 0.3-0.4 seconds to solve a problem of this size. Hence, we are looking at some 30 hours of core computations. That number is probably not going to change much if you call from c++ or something like that, as the core solver is the same (unless they have some poor internal overhead)

As a first experiment, I will model this as high-level as possible in YALMIP (disclaimer, developed by me) just to see how large the worst-case overhead is.

% Random data
n = 169;
m = 30;
A0 = randn(m,n);
B0 = randn(m,1);
C0 = randn(n,n)/10;
D0 = randn(n,n)/10;

The code here that would have to be changed for every new instance of data (defining objective and making a call to solvesdp to solve the problem) takes 0.8 seconds. Hence, there is indeed a substantial overhead of 100%, a large relative number, which is to be expected as the solution time is so low (this means that if you would have implemented this when you asked the question, you would have been done since long ago)

x = sdpvar(n,1);
opts = sdpsettings('solver','gurobi','verbose',0);
Objective = norm(A0*x-B0,2) + norm(C0*x,1) + norm(D0*x,inf);
solvesdp([],Objective,opts);

Doing some modelling manually brings the total down to 0.7 seconds

t = sdpvar(n,1);
s = sdpvar(1);
u = sdpvar(1);
Objective = u + sum(t) + s;
solvesdp([cone(A0*x-B0,u), -t <= C0*x <= t, -s <= D0*x <= s],Objective,opts);

Caching solvers on the MATLAB path saves us yet some time down to 0.6 seconds

opts = sdpsettings('solver','gurobi','verbose',0,'cachesolvers',1);
solvesdp([cone(A0*x-B0,u), -t <= C0*x <= t, -s <= D0*x <= s],Objective,opts);

However, similar problems for which some data changes should be preferably be attacked using the optimizer construct in YALMIP. This is a way to define parameterized models which YALMIP will preprocess and then finalize once data is given and a call is to be made to a solver.

The size of this model is on the limit of the capabilities of the optimizer framework, so it takes a minute or so to construct it (might even fail on the official release as I've made some fixes lately)

A = sdpvar(m,n,'full');
C = sdpvar(n,n,'full');
D = sdpvar(n,n,'full');
B = sdpvar(m,1);
x = sdpvar(n,1);
t = sdpvar(n,1);
s = sdpvar(1);
u = sdpvar(1);
Objective = norm(A*x-B,2) + norm(C*x,1) + norm(D*x,inf);
P = optimizer([],Objective,sdpsettings('solver','+gurobi'),{A,B,C,D},x);

Now run a loop where we solve 10 problems. This takes 4 seconds, of which 3.5 is spent inside Gurobi, i.e., overhead 10-15%

for i = 1:10;
    A0 = randn(m,n);
    B0 = randn(m,1);
    C0 = randn(n,n)/10;
    D0 = randn(n,n)/10;
    [xoptimal,~,~,~,P]=P{{A0,B0,C0,D0}};
end

EDIT: I saw now that only A and B changes. That brings down the overhead in the optimizer approach to 5%.

EDIT: The next release compiles the parametrized problem in 3 seconds. Just mail me if you need a pre-release.