example of relative entropy optimisation using exponential cone via matlab

225 views
Skip to first unread message

Ronan MT Fleming

unread,
Feb 11, 2021, 5:28:14 PM2/11/21
to mosek
Hi All,

does anyone have example code implementing

min  (d.*x)'*(log(x./y) + c)
s.t. l <= A[x;y] <= u

where d,c,A,l,u are data and x,y are variables, using the exponential cone solver in mosek 9.2?

I was previously using mskenopt for problems of the form

min  x'*(log(x) + c)
s.t.  Ax <= b

I have reviewed the general introduction to exponential cone optimisation
and there is one example with a matlab interface to a different exponential cone optimisation problem, but it is hard to adapt that to the first problem above.

An example with relative entropy optimisation would really help.

Regards,

Ronan

Utkarsh Detha

unread,
Feb 11, 2021, 6:24:39 PM2/11/21
to mosek
Hi,

You can reformulate the problem as follows:

min   d*t + d*c*x
s.t.   t >= x*log(x/y)
         l <= A[x;y] <= u

which is equivalent to:

min   d*t + d*c*x
s.t.   (y, x, -t) \in K_{exp}
         l <= A[x;y] <= u

Such a problem could be formulated using the Affine conic constraints, as shown in the following code:

function ceo1()

clear prob;

[r, res] = mosekopt('symbcon');
% Specify the non-conic part of the problem.

% x, y and auxiliary variable t
prob.c = [d*c 0 d];

% A = [1 1 0] for illustration
prob.a = sparse([1 1 0]);
prob.blc = l;
prob.buc = u;
prob.blx = [-inf -inf -inf];
prob.bux = [ inf inf inf];

% Specify the cones.
% The exponential cone: Fx + g \in K_{exp}
% where g is zero.
FE = sparse([0 1 0; 1 0 0; 0 0 -1]);
gE = [0 0 0];
cE = [res.symbcon.MSK_CT_PEXP 3];

prob.f = [FE];
prob.g = [gE];
prob.cones = [cE];

[r,res]=mosekopt('minimize',prob);

Please note: I have not tried running this code, so their might be some unchecked mistake. However, I hope it conveys the basic idea of the implementation that you can implement.

Hope this helps.

Regards,
Utkarsh

Ronan MT Fleming

unread,
Feb 16, 2021, 2:55:18 AM2/16/21
to mosek
Hi Utkarsh,

with some slight modifications, I have been able to solve an exponential cone problem following your suggestion. In your suggestion, the exponential cone variables were unbounded, but I had to bound them from below otherwise the problem seemed to be unbounded. The matrix of affine conic constraints (prob.F) is somewhat more complicated than a univariate example but it seems that permuting the rows works well to enforce those constraints. However, when I test for the optimality criteria, not all residuals are below the tolerance. In particular, I when I try to test if I have recovered an optimal solution to the original problem I posed, by checking the value of
 || t - x*log(x/y) ||_inf
I do not get something small as I had expected. Perhaps you can help pinpoint where I have strayed from the correct approach? 

Regards,

Ronan

%         https://docs.mosek.com/modeling-cookbook/expo.html
%         min  (d.*x)'*(log(x./y) + c)
%         s.t. l <= A[x;y] <= u
%        
%         where d,c,A,l,u are data and x,y are variables, is equivalent to
%        
%         min   d*t + d*c*x
%         s.t.   t >= x*log(x/y)
%         l <= A[x;y] <= u
%        
%         which is equivalent to:
%        
%         min   d*t + d*c*x
%         s.t.   (y, x, -t) \in K_{exp}
%         l <= A[x;y] <= u
%        
%         Such a problem could be formulated using the Affine conic constraints, as shown in the following code:


       
        % Specify the non-conic part of the problem.
        [mint,nint] = size(N);
       
        On = sparse(nint,nint);
        Omn = sparse(mint,nint);
        Om1 = sparse(mint,1);
        On1 = sparse(nint,1);
        In = speye(nint);
        I1n = ones(1,nint);
        In1 = ones(nint,1);
        e   = ones(nint,1);

        % vf, vr, s and auxiliary variables tf & tr
        d = ones(nint,1);
        prob.c = [cf; cr; 0; d; d];


       % x = [vf;vr]
       % y = s;
        % t = [tf;tr];

        %        vf    vr    s     tf      tr
        A  = [   N     -N   Om1   Omn     Omn;
                     In    -In   On1        On       On;
                  -I1n   -I1n     1      O1n     O1n];
         
        prob.a = A;
        prob.blc = [b;l;0];
        prob.buc = [b;u;0];
        if 0
            prob.blx = zeros(4*nint+1,1);
        else
            prob.blx = -inf*ones(4*nint+1,1); %unbounded problem
            prob.blx = [zeros(2*nint+1,1);-inf*ones(2*nint,1)]; %unbounded problem
            prob.blx = [zeros(2*nint+1,1);-10*ones(2*nint,1)]; %unbounded problem
        end
        prob.bux =  inf*ones(4*nint+1,1);
       
        % Specify conic part of the problem
        % https://docs.mosek.com/9.2/toolbox/data-types.html#cones
       
        % Each cone is a primal exponential cone.
        % x3 <= x2 * log(x1 /x2), x1, x2 >= 0  <=> [x1;x2;x3] \in K_{exp}
        
        % The cones are specified using two index lists cones.subptr and cones.sub and list of cone-type identifiers cones.type.
        % For affine conic constraints Fx+g \in K_exp, where K = K_1 * ... * Ks, cones is a list consisting
        % of s concatenated cone descriptions.
        [~, res] = mosekopt('symbcon');
        if 1
            prob.cones(1:2:4*nint) = res.symbcon.MSK_CT_PEXP;
            prob.cones(2:2:4*nint) = 3;
        else
            %*** Error(1301): Wrong number of variables in the cone. Exactly 3 members are required for cones of type MSK_CT_PEXP.
            prob.cones(1:2*nint,1) = res.symbcon.MSK_CT_PEXP;
            prob.cones(1:2*nint,2) = 3;
        end
       
        %Conic problem with affine conic constraints
        %https://docs.mosek.com/9.2/toolbox/data-types.html#equation-doc-notation-conic
        %f (double[][]) – The matrix of affine conic constraints. It must be a sparse matrix.
        F =      [In On  On1  On  On;
                  On In  On1  On  On;
                  On On  In1  On  On;
                  On On  In1  On  On;
                  On On  On1 -In  On;
                  On On  On1  On -In];
       
        %permute the rows of F to form (y, x, -t) triples for each exponential cone
        ind(1:3:(6*nint)) = (2*nint+1:4*nint)';
        ind(2:3:(6*nint)) = (1:2*nint)';
        ind(3:3:(6*nint)) = (4*nint+1:6*nint)';
       
        prob.f = F(ind,:);     
        %g (double[]) – The constant term of affine conic constraints. If not present or g==[] it is assumed g=0
        if 1
            prob.g = zeros(size(prob.f,1),1);
        else
            %*** Error(1200): prob.g must be a dense vector with length(prob.g) == size(prob.f,1)
            prob.g = [];
        end
       
        [r,res]=mosekopt('minimize',prob);
       
       
        if strcmp(res.sol.itr.prosta,'PRIMAL_AND_DUAL_FEASIBLE') &&  strcmp(res.sol.itr.solsta,'OPTIMAL')
            stat=0;
        else
            fprintf('%s\t%s\t%s\n','fixedPointSolverSub:', res.sol.itr.prosta,res.sol.itr.solsta);
            stat=1;
        end

        % Primal variables
        x = res.sol.itr.xx;
        % Dual variables corresponding to constraints
        y = res.sol.itr.y;
        % Dual variables corresponding to bounds on variables.
        z = res.sol.itr.slx - res.sol.itr.sux;
        % Dual variables of affine conic constraints
        w = res.sol.itr.doty;
       
        fprintf('%s\n','Optimality conditions')
        fprintf('%6.2g\t%s\n',norm(prob.c - prob.a'*y - z - prob.f'*w,inf), '|| c - A''*y - z - F''w ||_inf');
        fprintf('%6.2g\t%s\n',norm(x(2*nint+2:4*nint+1) - x(1:2*nint).*log(x(1:2*nint)/x(2*nint+1)),inf), '|| t - x*log(x/y) ||_inf');


gives

Optimality conditions
3.7e-08    || c - A'*y - z - F'w ||_inf
   7.2    || t - x*log(x/y)

Michal Adamaszek

unread,
Feb 16, 2021, 3:34:23 AM2/16/21
to mosek
If you do mosekopt('write problem.opf',prob) then you will get a human readable file, where you can inspect if the problem is really what you intended to model, preferable on small data. This should help you make sure if you do not have some modeling error. Note that the affine conic constraints will appear via slacks (see "Remark" in https://docs.mosek.com/9.2/toolbox/prob-def-affine-conic.html)

If on the other hand you think it is a numerical issue then please prepare reproducible data and send to support or post everything here, as you prefer.

Michal

Ronan MT Fleming

unread,
Feb 16, 2021, 8:38:18 AM2/16/21
to mosek
Hi Michal,

I have tried to strip it down to a very simple example, but I still do not get what I expect. It may be my lack of experience with conic programming.

I try to solve the problem
min f*log(f/s) + r*log(r/s)
s.t.  -f + r = 0.1
        -f -r + s = 0
        0 <= f,r,s

Which I understand to be the following exponential cone problem

min tf + tr
s.t.  -f + r = 0.1
        -f -r + s = 0
        0 <= f,r,s,tf,tr
        [s;f;-tf] \in K_exp
        [s;r;-tr] \in K_exp

Which in mosek format, detailed below from  mosekopt('write(problem.opf)',prob), is:

min tf + tr
s.t.  -f + r = 0.1
        -f -r + s = 0
        0 <= f,r,s,tf,tr
           - s + yf = 0
           - f + xf = 0
            tf + mtf = 0
            - s + yr = 0
            - r + xr = 0
            tr + mtr = 0
        [yf;xf;mtf] \in K_exp
        [yr;xr;mtr] \in K_exp

When I check the dual optimality conditions, I see that

        % Primal variables
        x = res.sol.itr.xx;
        % Dual variables corresponding to constraints
        y = res.sol.itr.y;
        % Dual variables corresponding to bounds on variables.
        z = res.sol.itr.slx - res.sol.itr.sux;
        % Dual variables of affine conic constraints
        w = res.sol.itr.doty;
        fprintf('%6.2g\t%s\n',norm(prob.c - prob.a'*y - z - prob.f'*w,inf), '|| c - A''*y - z - F''w ||_inf');
gives 
8.7e-09    || c - A'*y - z - F'w ||_inf
as expected.

However, as tf replaced f.* log(f / s) and tr replaced r.* log(r / s) in the transition to an exponential cone problem I had expected the following to be almost zero, but they are not:
  0.51    || tf - f.* log(f / s) ||_inf
  0.54    || tr - r.* log(r / s) ||_inf

I am at a loss as to why this is not the case. Any suggestion why not would be most appreciated.

Regards,

Ronan


[comment]
   Written by MOSEK version 9.2.37
   Date 16-02-21
   Time 12:08:20
[/comment]

[hints]
  [hint NUMVAR] 11 [/hint]
  [hint NUMCON] 10 [/hint]
  [hint NUMANZ] 21 [/hint]
  [hint NUMQNZ] 0 [/hint]
  [hint NUMCONE] 2 [/hint]
[/hints]

[variables disallow_new_variables]
  f r s tf tr yf xf mtf yr xr mtr
[/variables]

[objective minimize]
   tf + tr
[/objective]

[constraints]
  [con 'c000000000_']  - f + r = -1e-1 [/con]
  [con 'c000000001_']  f - r = 1e-1 [/con]
  [con 'c000000002_']  f - r [/con]
  [con 'c000000003_']  - f - r + s = 0 [/con]
  [con 'c000000004_']  - s + yf = 0 [/con]
  [con 'c000000005_']  - f + xf = 0 [/con]
  [con 'c000000006_']  tf + mtf = 0 [/con]
  [con 'c000000007_']  - s + yr = 0 [/con]
  [con 'c000000008_']  - r + xr = 0 [/con]
  [con 'c000000009_']  tr + mtr = 0 [/con]
[/constraints]

[bounds]
  [b] 0e+00      <= f,r,s,tf,tr [/b]
  [b]               yf,xf,mtf,yr,xr,mtr free [/b]
  [cone pexp 'k000000005_'] yf, xf, mtf [/cone]
  [cone pexp 'k000000009_'] yr, xr, mtr [/cone]
[/bounds]

Ronan MT Fleming

unread,
Feb 16, 2021, 8:39:51 AM2/16/21
to mosek
Hi Michal,

thanks for your suggestion. I have tried to strip it down to a very simple example. However, I still cannot see what is wrong. It may be my lack of experience of conic programming, so possibly it is a simple mistake.

I try to solve this problem:

min  f*log(f/s) + r*log(r/s)
s.t.  -f + r  = 0.1 : y
       - f - r  + s = 0 : p
        0<= f, r, s : q

With dual optimality condition
log(f/s) + 1 + y  + q - q = 0
log(r/s) + 1 - y  + q - q = 0
Giving
log(r/f) = 2y
as desired for modelling purposes in biochemistry.

In terms of an exponential cone problem, I understand the optimisation problem to be:

min  tf + tr
s.t.  -f + r  = 0.1
       - f - r  + s = 0
       0 <= f, r, s, tf, tr
       [s; f; - tf] \in K_{exp}
       [s; r; - tr] \in K_{exp}

And in mosek format is:

min  tf + tr
s.t.  -f + r  = 0.1
       - f - r  + s = 0
       0 <= f, r, s, tf, tr
         - s + yf = 0
         - f + xf = 0
       tf + mtf = 0
         - s + yr = 0
         - r + xr = 0
        tr + mtr = 0
       [yf; xf; mtf] \in K_{exp}
       [yr; xr; mtr] \in K_{exp}

When I solve the problem, specified below using mosekopt('write(problem.opf)',prob), I can see that part of the dual optimality conditions are satisfied:

        % Primal variables
        x = res.sol.itr.xx;
        % Dual variables corresponding to constraints
        y = res.sol.itr.y;
        % Dual variables corresponding to bounds on variables.
        z = res.sol.itr.slx - res.sol.itr.sux;
        % Dual variables of affine conic constraints
        w = res.sol.itr.doty;
        
fprintf('%6.2g\t%s\n',norm(prob.c - prob.a'*y - z - prob.f'*w,inf), '|| c - A''*y - z - F''w ||_inf');
8.7e-09 =    || c - A'*y - z - F'w ||_inf
     0   =  || -y + z ||_inf

But I had anticipated that, because tf is a replacement for f*log(f/s), at optimality, the following should be zero, but are not:
  0.51    || tf - f.* log(f / s) ||_inf
  0.54    || tr - r.* log(r / s) ||_inf

Any suggestions to remedy the situation are most welcome.

Michal Adamaszek

unread,
Feb 16, 2021, 9:19:43 AM2/16/21
to mosek
The condition tf,tr>=0 is not needed. In fact, your entropy tends to be negative, and therefore further constraining tf>=0 forces this constraint to be binding and the inequality with log to be non-binding. So in the solution, as you can see, you get tf=0.

Without these bounds Mosek correctly identifies that the problem is unbounded (objective tends to -inf as f goes to +inf as one can see from the plot)

Ronan MT Fleming

unread,
Feb 17, 2021, 12:42:08 AM2/17/21
to mosek
Hi Michal,

thanks a lot for that. By box constraining the value of y, the problem is bounded and I could verify (some of) the optimality conditions involving the exponential cone.

What I am missing, and I cannot seem to find it on the mosek website, is a matlab expression of the optimality conditions involving the exponential cone that does not involve an abstract representation of an exponential cone.

According to the definition of a problem with an affine conic constraint here

I currently use the following code to test dual optimality criteria:

        [r,res]=mosekopt('minimize',prob);
               
        x = res.sol.itr.xx;
        y = res.sol.itr.y;

        z = res.sol.itr.slx - res.sol.itr.sux;
        w = res.sol.itr.doty;
       
            %
            fprintf('%6.2g\t%s\n',norm(prob.c - prob.a'*y - z - prob.f'*w,inf), '|| c - A''*y - z - F''w ||_inf');

            fprintf('%6.2g\t%s\n',norm(-y + res.sol.itr.slc - res.sol.itr.suc,inf), '|| -y + z ||_inf');

            fprintf('%6.2g\t%s\n',norm(x(2*n+2:4*n+1) - x(1:2*n).*log(x(1:2*n)/x(2*n+1)),inf), '|| t - x*log(x/y) ||_inf');

            fprintf('%6.2g\t%s\n',norm(F(:,1:2*n)'*w + log(x(1:2*n)/x(2*n+1)) + 1,inf), '|| F''w - log(x/y) + 1||_inf');

However, I have had to reverse engineer the latter two lines of code to test part of the optimality criteria and I am missing the rest. Also, I am not sure if they are entirely correct, because the values seem to deviate from zero for problems that are more numerically challenging.

I would like to be able to write an expression in matlab that allows me to recover the PFEAS and  DFEAS as reported in the iteration log (example below). I understand that the actual problem being solved is (slightly different from the "vanilla" version (https://docs.mosek.com/9.2/toolbox/prob-def-affine-conic.html). However, I do not see in the res.sol.itr structure where I might obtain τ and κ, the two additional scalar variables.

Once I can express in matlab exactly what the primal and dual optimality criteria are of the actual problem being solved I can use that to evaluate when numerical issues are starting. This interface work is to bring mosek into the core of a highly used simulation tool in systems biology (https://opencobra.github.io/cobratoolbox/stable/). All sorts of biochemical models might be thrown at the solver so I need to be able to alert the unsuspecting end user when there may be numerical issues.

Regards,

Ronan


----------------
MOSEK Version 9.2.37 (Build date: 2021-2-8 10:03:08)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: Linux/64-X86

Problem
  Name                   :                
  Objective sense        : min            
  Type                   : CONIC (conic optimization problem)
  Constraints            : 63626          
  Cones                  : 16576          
  Scalar variables       : 84686          
  Matrix variables       : 0              
  Integer variables      : 0              

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 2                 time                   : 0.00           
Lin. dep.  - tries                  : 1                 time                   : 0.01           
Lin. dep.  - number                 : 97             
Presolve terminated. Time: 0.11   
Problem
  Name                   :                
  Objective sense        : min            
  Type                   : CONIC (conic optimization problem)
  Constraints            : 63626          
  Cones                  : 16576          
  Scalar variables       : 84686          
  Matrix variables       : 0              
  Integer variables      : 0              

Optimizer  - threads                : 4              
Optimizer  - solved problem         : the primal     
Optimizer  - Constraints            : 26598
Optimizer  - Cones                  : 16576
Optimizer  - Scalar variables       : 56041             conic                  : 49728          
Optimizer  - Semi-definite variables: 0                 scalarized             : 0              
Factor     - setup time             : 0.13              dense det. time        : 0.00           
Factor     - ML order time          : 0.02              GP order time          : 0.00           
Factor     - nonzeros before factor : 2.06e+05          after factor           : 2.70e+05       
Factor     - dense dim.             : 3                 flops                  : 1.17e+07       
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME 
...
39  5.0e-08  1.9e-11  2.6e-10  1.00e+00   -9.715680642e+03  -9.715680709e+03  1.5e-11  3.07 
40  5.0e-08  1.9e-11  2.6e-10  1.00e+00   -9.715680644e+03  -9.715680711e+03  1.5e-11  3.14 
41  2.6e-08  4.7e-12  3.1e-11  1.00e+00   -9.715703865e+03  -9.715703881e+03  3.6e-12  3.26 
Optimizer terminated. Time: 3.39   


Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -9.7157038651e+03   nrm: 1e+03    Viol.  con: 2e-05    var: 0e+00    cones: 0e+00 
  Dual.    obj: -9.7157038810e+03   nrm: 1e+03    Viol.  con: 6e-13    var: 1e-09    cones: 0e+00 
Optimizer summary
  Optimizer                 -                        time: 3.39   
    Interior-point          - iterations : 41        time: 3.29   

Erling D. Andersen

unread,
Feb 18, 2021, 2:01:42 AM2/18/21
to mosek
Recovering PFEAS and DFEAS make no sense in my opinion and is impossible in any way. Since these are some internal
algorithmic numbers that largely only make sense for algorithm.

I assume what you want is to verify optimality of the solution Mosek reports.

If you have your own optimality conditions you will somehow have to figure out the relationship between those and the one Mosek employs.
You have to transform the Mosek solution to be solution to your optimality conditions. That may require a bit work.

In general it is easier to work we conic duality as we do in Mosek. Conic duality is a generalization of the usual LP duality. Indeed we state the dual problem at


so you should verify that the dual solution Mosek reports is a feasible solution to that problem. It should not be too hard for you.
Next you verify that the primal solution is feasible to the primal primal. Finally, if primal objective value is almost identical to dual objective value then you have proved optimality.

That is in fact what is going at

Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -9.7157038651e+03   nrm: 1e+03    Viol.  con: 2e-05    var: 0e+00    cones: 0e+00 
  Dual.    obj: -9.7157038810e+03   nrm: 1e+03    Viol.  con: 6e-13    var: 1e-09    cones: 0e+00 

It shows that the primal and dual solutions are nearly feasible and about 8 figures in the objective value is correct.
So looks like we computed a pretty decent solution.

We are looking forward to see a lot of Mosek usage in relation to Cobra toolbox. And maybe some sales too maybe :-).
If you/we can generate some test cases using the Cobra toolbox, then we would interested in that. 

Hope this helps.

PS. For those who read this and may want to learn more about conic optimization then consult

Reply all
Reply to author
Forward
0 new messages