Hello everyone, I have a question in optimization

Please find the attached question

[John D'Errico]

I suppose this is not a homework assignment. However, all three models were identical as you wrote them, so your question is totally confused. In any case, use GA.

This is the modet appropriate tool, since you use 1 norms - that makes the problem a non-differentiable one. As well, the presence of an infinity norm in the constraints again will cause problems for other optimizers. Just use GA.

Since models 2 and 3 as written are identical to the first, it is impossible to know more.

[Bruno Luong]

In this pdf: r_m c_m cannot be vectors since they are the bounds of the norms, which must be scalars.

[Ahmed Galal]

Sorry, they are not identical only the constraints are identical. The objective function is L1, L2, and Linf, not only L1. I did a mistake

Answer 1

[Bruno Luong]

The constraint (1c)

norm(R11∗*w+r00, Inf) >= rm

can be transformed as a union of 42 halfplanes

R11(i,:)*w+r00(i,:) >= rm 

or

R11(i,:)*w+r00(i,:) <= -rm

for i=1,2,...21.

I would then suggest to solve 42 sub linear-programing problems by replacing the (1c) with one of those conditions. The sub problem can be solved with intlinprog, then we just take the argmin of those 42 problem solutions.

It must be more preditable and robust than GA.

[Matt J]

@Ahmed,

Once you replace (1c) with the half-plane constraint mentioned by Bruno and convert all the other constraints to matrix form, the resulting sub-problems can be solved with lsqlin, minl1lin, and fminimax (for the L2,L1, and L-inf norms respectively), as we discussed in an earlier thread.

The matrix form for the L1-constraint (1d) can be obtained as follows,

N=21;

clear A;
[A{1:N}]=ndgrid([-1,1]/sqrt(N));
A=reshape( cat(N+1,A{:}) ,[],N );

b=cm*abs(A(:,1))-A*c00;
A=A*C11;  %Constraint (1d) is the same as A*w<=b

However, in your case, this is a problem with more than 2 million linear constraints. It will be interesting to see if any of the Optimization Toolbox solvers will handle that many.

[Bruno Luong]

Here is a code for l^2-norm using QUADPROG for a toy example.

(this might be the last fish you might get from me).

Up to you to do l1 or linfinity norm (hint uing LINPROG or other tools we suggested you in other threads) .

EDIT: some bugs are found and corrected

N1=eye(2);
N0=[2; 4];
nh=-1000;
H1=[1 0; 0 1];
u1=[-0.7; 10];
R11=eye(2);
r00=[0; 0];
rm=2;
C11=eye(2);
c00=zeros(2,1);
cm=3;

% quadratic form and linear form where the sum is to be minimized
n = size(N1,2);
H = N1'*N1;
f = N0'*N1;
%
A = [-N1;  % constraint (a)
      H1]; % constraint (b)
b = [N0 - nh;
     u1];
% m slack variables to force L1-norm constraint (d)
m = size(C11,1);
As = [ C11 -eye(m);
      -C11 -eye(m)];
bs= [-c00;
     +c00];
Ad = [zeros(1,n) ones(1,m)];
bd = cm;
A = [A zeros(size(A,1),m);
     As;
     Ad];
b = [b;
     bs;
     bd];  
H(end+m,end+m) = 0;
f(end+m) = 0; % slack variables matches vector bounds, probably can be set the m-last elements to 1s

% Loop for constraint (c), decomposing
fmin = Inf;     
for i=1:size(R11,1)
    for s=[-1 1]     
        Ack = [-s*R11(i,:), zeros(1,m)];
        bck = s*r00(i)-rm;
        Ak = [A; Ack];
        bk = [b; bck];
        [xk, fk, flag] = quadprog(H, f, Ak, bk);
        if flag == -3
            fmin = -Inf;
            break
        elseif flag >= 0 && fk < fmin
            xmin = xk;
            fmin = fk;
        end
    end
    if flag < 0 && flag ~= -2 % not fesasible
        break
    end
end
if isfinite(fmin)
    w = xmin(1:n);
    obj = norm(N1*w + N0);
    tol = 1e-6;
    constraints = [all(N1*w + N0        >= nh-tol), ...	(a)
                   all(H1*w - u1        <= 0+tol),  ...	(b)
                   norm(R11*w+r00,Inf)  >= rm-tol,  ...	(c)
                   norm(C11*w+c00,1)    <= cm+tol   ...	(d)
                   ];
    fprintf('constraints satisfied = %s\n', mat2str(constraints));
    fprintf('solution = %s\n', mat2str(w));
    fprintf('minnorm = %g\n', obj);
    if ~all(constraints)
        fprintf('constraints not satisfied (bug?)\n')
    end
elseif fmin == -Inf
    fprintf('problem unbounded\n')
else
    error('No solution found')
end

[Bruno Luong]

@Matt: "However, in your case, this is a problem with more than 2 million linear constraints. It will be interesting to see if any of the Optimization Toolbox solvers will handle that many."

See my code for usage of slack variables to limit the matrix of l1-norm constraint.

[Bruno Luong]

You can carry on pretty straigh forward with other inputs Aeq, beq, lb, ub of quadprog. If you need the syntax and argument here it is.

Remember, I though I never giver you a fish again. However you are free to ask question on my code though.

[Ahmed Galal]

Thank you very much Mr.Bruno

This is the result when I used your code for my data, i donot where is is mistake

Error using quadprog (line 285)

The number of rows in A must be the same as the number of elements of b.

Error in Neighbouring_criterion_burnL1L2LINF (line 787)

[xk, fk, flag] = quadprog(HH, f1, Ak, bk,[],[],lb,ub,w0);

[Bruno Luong]

In my formulation I add m (=21) slack variables to the tail of the main w decision parameters. The quadprog formulation works on vector of length 56:

x = [w; s]; with size(w,1)=35 and m=size(s,1)=21.

If you want to impose lower bounds for w, lbw of length(=35) you need to add 21 -inf to the tail, and similarly for the upper bounds

lbx = [lbw; -inf(m,1)]
ubx = [ubw; inf(m,1)]
[xk, fk, flag] = quadprog(HH, f1, Ak, bk,[],[],lbx,ubx);

You might need also to make some change of w0 (I have no idea what you use there).

[Matt J]

I like the use of slack variables, but I do not see the need for quadprog. Once you have all the constraints set up linearly, you can use lsqlin for the L2 case. You can also use minl1lin for the L1 case and fminimax for L_inf as already discussed.

[rn1,cn1]=size(N1)
C=[N1,zeros(rn1,m)]; C(rn1+m,end)=0;  %pad for slack
d=[-N0;zeros(m,1)];

 x=lsqlin(C,d,Ak, bk,[],[],lbx,ubx); 
w=x(1:cn1);  %L2

 x=minl1lin(C,d,Ak, bk,[],[],lbx,ubx);  
w=x(1:cn1); %L1

 options = optimoptions('fminimax','AbsoluteMaxObjectiveCount',1);
x=fminimax(@(x)C*x-d, rand(m+cn1,1),Ak, bk,[],[],lbx,ubx);  
w=x(1:cn1);  %L-inf

[Bruno Luong]

What's wrong with quadprog (compared to lsqlin)? Is lsqlin preferable for semi-convex problem?

Reading the doc I have impression they have the same engine under the hood, just the interface is different.

[Matt J]

I was under the impression that it is numerically better to evaluate gradients as,

C.'*(C*x-d)  %lsqlin

instead of,

(C.'*C)*x-C.'*d  %quadprog

As you get close to a solution point where C*x-d=0, the former version would do the subtraction on numbers of smaller magnitude and conserve precision.

[Bruno Luong]

In LSQLIN doc page, the algorithm part at the bottom, it just says it uses quadprog active-set and interior-point (default method). Not clear about Trust-Region-Reflective Algorithm, they both link to the same paper.

If you edit the code of lsqlin.m and follow interior-point branch or active-set branch, you see that computes H and f exactly like in my code; so no the gradient does not use parethesis like your.

        % compute for later use in algorithm
        f = -C'*d;
           ...
            H = C'*C;
            [X, ~, exitflag, output, lambda] = ...
                ipqpdense(H, f, ...);

But OK if you feel safer and more comfortable with lsqlin, why not.

[Matt J]

so no the gradient does not use parethesis like your.

OK, but they should do it that way, no? They missed the opportunity to take advantage of the special form of the problem.

[Bruno Luong]

Does it really matter for accuracy? At the cost of 2 matrix-vector product per gradient evaluation instead of one?