How to implement following function in CVX?

I have a function of form which I have observed as a convex optimization problem w.r.t {\tau_0,\tau_1}.

Knowing, rel_entr(x,y)=x.*log(x/y) and -rel_entr( x, x+y)= xlog( 1+y/x ).

How to convert my function, so that CVX accepts the function?
I have tried but couldn’t succeed. Help is highly appreciated.

I think you should be able to adapt my solution at Xlog( 1+ Y/(X+Y) ): DCP rules, and build-in functions in cvx to accommodate the constants in your function.

Yes. I also think that my equation can be transformed to the form that CVX accepts. I am still trying to transform to the suitable form. If anyone can figure out suitable transformation, please let me know. Thanks.

My linked answer provides the suitable form: it can be expressed using rel_entr. I think that should work, provided all the constants are nonnegative.

\tau_1 log\left( 1 + \frac{\mathbf {A_i} \tau_0 }{\mathbf{B_i} \tau_0 +C_i \tau_1 }\right) =\tau_1 log\left( 1 + \frac{ \tau_0 }{ \frac{\mathbf{B_i}}{\mathbf {A_i}} \tau_0 +\frac{C_i}{\mathbf {A_i}} \tau_1 }\right)
According to your solution at (Xlog( 1+ Y/(X+Y) ): DCP rules, and build-in functions in cvx) my function can be transformed as
(\tau_0+\tau_1) log\left( 1 + \frac{ \tau_0 }{ \frac{\mathbf{B_i}}{\mathbf {A_i}} \tau_0 +\frac{C_i}{\mathbf {A_i}} \tau_1 }\right) - \tau_0 \ log_2\left( 1 + \frac{ \tau_0 }{ \frac{\mathbf{B_i}}{\mathbf {A_i}} \tau_0 +\frac{C_i}{\mathbf {A_i}} \tau_1 }\right).

I noticed that the second term, - \tau_0 \ log\left( 1 + \frac{ \tau_0 }{ \frac{\mathbf{B_i}}{\mathbf {A_i}} \tau_0 +\frac{C_i}{\mathbf {A_i}} \tau_1 }\right), can be written as -\text{rel_entr}( \frac{C_i}{A_i} \tau_1 + \frac{B_i}{A_i} \tau_0, \frac{C_i}{A_i} \tau_1 + \frac{B_i}{A_i} \tau_0+ \tau_0) - \text{rel_entr}( \frac{C_i}{A_i} \tau_1 + \frac{B_i}{A_i} \tau_0+ \tau_0, \frac{C_i}{A_i} \tau_1 + \frac{B_i}{A_i} \tau_0) .

How can we change the first term? Help is highly appreciated. Thanks.

Mu solution works when A = B = c = 1. You just need to adapt it by looking at the derivation, to handle other values of A,B,c.

You will also make things easier for everyone if you show formulas as code, not images. Also, it’s easy to check your work by picking random values of the inputs and seeing if the numerical value of the reformulation matches the original.

I noticed that in CVX, I can write my function \tau_1 \ log \left( 1 + \frac{ \mathbf{A_i} \tau_0 }{ c_i \tau_1 + \mathbf {B_i} \tau_0 } \right) as

- rel_{-} \operatorname{entr} \left( \tau_1+ \frac{ \mathbf{B_i} }{c_i} \tau_0 , \tau_1+ \frac{ \mathbf {B_i} }{c_i} \tau_0 + \frac{ \mathbf{A_i}}{ c_i} \tau_0 \right) -rel_{-} \operatorname{entr}( \frac{ \mathbf {B_i} }{ \mathbf {A_i} } \tau_1 + \frac{ \mathbf {B_i}^2}{ \mathbf {A_i} c_i} \tau_0, \frac{ \mathbf {B_i}}{ \mathbf {A_i}} \tau_1 + \frac{ \mathbf {B_i}^2}{ \mathbf {A_i} c_i} \tau_0+\frac{ \mathbf {B_i}}{ c_i} \tau_0)
- rel_{-} \operatorname{entr}( \frac{ \mathbf {B_i}}{ \mathbf {A_i}} \tau_1 + \frac{ \mathbf {B_i}^2}{ \mathbf {A_i} c_i} \tau_0+ \frac{ \mathbf {B_i} }{ c_i}\tau_0, \frac{ \mathbf{B_i} }{ \mathbf{A_i}} \tau_1 + \frac{ \mathbf{B_i}^2}{ \mathbf{A_i} c_i} \tau_0) .

Then, I wrote the following code for my optimization problem.

for i=1:K
cvx_begin
variable tau0 nonnegative
tau1=1-tau0;
R1=-rel_entr(tau1+ ((B(i)./c(i)) .* tau0 ), tau1+((B(i)+A(i))./c(i)) .* tau0 );

       R2=-rel_entr( (B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0, ...
(B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0+ (B(i) ./ c(i)).*tau0) ;

R3=-rel_entr((B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0+ (B(i) ./ c(i)).*tau0,...
(B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0) ;

maximize (tau1.*(R1+R2+R3));
subject to
tau0<=1;
cvx_end
tau0_opt(i)=tau0;
end


Expressions R1, R2, R3, and R4 in my code are concave expressions. Hence, the objective function is a product of an affine and concave expression. Thus, when running this code, I am getting an error which is as follows.

Error using .* (line 173)
Disciplined convex programming error:
Cannot perform the operation: {real affine} .* {concave}

How to modify my code to handle this issue in CVX? Help required. Thanks in advance.

I’m confused as to what mathematical problem you are trying to solve.

Do you have the expression involving the product of tau and log(1+…) in terms of sum of -rel_entr() ? If so, why are you attempting to multiply the -rel_entr terms by tau again?

As a check, pick some random numbers for all the inputs and verify the original non-DCP compliant formulation produces the same result, within roundoff error, as the reformulation using -rel_entr.

K=4;
A =[ 24.0173 17.6637 78.0226 60.4402];
B =[ 156.1265 162.4801 102.1212 119.7036];

 cvx_begin
variable tau0 nonnegative
tau1=1-tau0;
for i=1:K
R1=-rel_entr(tau1+ ((B(i)./c(i)) .* tau0 ), tau1+((B(i)+A(i))./c(i)) .* tau0 );

R2=-rel_entr( (B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0, ...
(B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0+ (B(i) ./ c(i)).*tau0) ;

R3=-rel_entr((B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i)) ).*tau0+ (B(i) ./ c(i)).*tau0,...
(B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0) ;
end
maximize (R1+R2+R3);
subject to
tau0<=1;
cvx_end


CVX Warning:
Models involving “rel_entr” or other functions in the log, exp, and entropy
family are solved using an experimental successive approximation method.
This method is slower and less reliable than the method CVX employs for
other models. Please see the section of the user’s guide entitled
The successive approximation method
for more details about the approach, and for instructions on how to
suppress this warning message in the future.

Cones | Errors | Mov/Act | Centering Exp cone Poly cone | Status --------±--------------------------------±-------- 3/ 3 | 1.452e+00 1.422e-01 0.000e+00 | Failed 3/ 3 | 1.966e-01 3.158e-03 0.000e+00 | Failed 2/ 2 | 1.044e-02 6.073e-06 0.000e+00 | Failed 0/ 2 | 8.000e+00 1.399e+06 1.399e+06 | Failed 0/ 0 | 0.000e+00 0.000e+00 0.000e+00 | Failed 0/ 2 | 8.000e+00 1.417e+08 1.417e+08 | Failed

Status: Failed
Optimal value (cvx_optval): NaN

For my code , CVX failed to solve the problem. And I have noticed a suggestion for such issues. Hence I installed CVXQUAD.
As I noticed, rel_entr_quad(x,y)=x.*log(x./y). I replaced rel_entr(x,y) in the above program as rel_entr_quad(x,y). Then, again problem is not solved. And obtained as follows.

Calling SDPT3 4.0: 242 variables, 97 equality constraints For improved efficiency, SDPT3 is solving the dual problem.

num. of constraints = 97
dim. of sdp var = 144, num. of sdp blk = 72
dim. of linear var = 26

SDPT3: Infeasible path-following algorithms

number of iterations = 27 primal objective value = -9.20978219e+00 dual objective value = -9.20895652e+00 gap := trace(XZ) = 3.95e-03 relative gap = 2.03e-04 actual relative gap = -4.25e-05 rel. primal infeas (scaled problem) = 1.45e-07 rel. dual " " " = 1.63e-07 rel. primal infeas (unscaled problem) = 0.00e+00 rel. dual " " " = 0.00e+00 norm(X), norm(y), norm(Z) = 1.0e+03, 1.9e+04, 2.5e+04 norm(A), norm(b), norm© = 3.1e+01, 2.8e+01, 6.7e+00 Total CPU time (secs) = 0.60 CPU time per iteration = 0.02 termination code = -5 DIMACS: 2.0e-07 0.0e+00 5.5e-07 0.0e+00 -4.3e-05 2.0e-04

Status: Failed
Optimal value (cvx_optval): NaN

Anyway as per the suggestion here, no modifications are needed for my original code. Still I cannot understand why my problem is not solved.

Your program is not reproducible, because it is missing input data. Therefore, i can’t run it.

One thing you can try yourself is to specify sedumi as the solver.

Better yet, if you have access to Mosek 9.x, restore CVX’s version of exponential.m, make sure you are using CVX 2.2, then start a fresh MATLAB session. Sspecify Mosek as solver, and run it. That will utilize Mosek 9.x’s native exponential cone capability, which is the most robust and reliable of the solution options for this problem. It will also print our warnings and other useful diagnostic information which Mosek employees who read this forum can assess.

K=4;
A =[ 24.0173 17.6637 78.0226 60.4402];
B =[ 156.1265 162.4801 102.1212 119.7036];
c=1.0000e-08*ones(1,4);

cvx_begin
variable tau0 nonnegative
tau1=1-tau0;
for i=1:K
R1=-rel_entr(tau1+ ((B(i)./c(i)) .* tau0 ), tau1+((B(i)+A(i))./c(i)) .* tau0 );

R2=-rel_entr( (B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0, ...
(B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0+ (B(i) ./ c(i)).*tau0) ;

R3=-rel_entr((B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i)) ).*tau0+ (B(i) ./ c(i)).*tau0,...
(B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0) ;
end
maximize (R1+R2+R3);
subject to
tau0<=1;
cvx_solver mosek
cvx_end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mosek solver result is as follows.

Calling Mosek 9.1.9: 38 variables, 13 equality constraints For improved efficiency, Mosek is solving the dual problem.

MOSEK Version 9.1.9 (Build date: 2019-11-21 11:34:40)
Platform: Windows/64-X86

MOSEK warning 62: The A matrix contains a large value of -1.6e+10 in constraint ‘’ (12) at variable ‘’ (2).
MOSEK warning 62: The A matrix contains a large value of -1.0e+11 in constraint ‘’ (12) at variable ‘’ (5).
MOSEK warning 62: The A matrix contains a large value of 4.3e+10 in constraint ‘’ (12) at variable ‘’ (7).
MOSEK warning 62: The A matrix contains a large value of -1.2e+11 in constraint ‘’ (12) at variable ‘’ (8).
MOSEK warning 62: The A matrix contains a large value of 3.7e+10 in constraint ‘’ (12) at variable ‘’ (10).
MOSEK warning 62: The A matrix contains a large value of -1.6e+10 in constraint ‘’ (12) at variable ‘’ (11).
MOSEK warning 62: The A matrix contains a large value of -1.5e+11 in constraint ‘’ (12) at variable ‘’ (14).
MOSEK warning 62: The A matrix contains a large value of 6.1e+10 in constraint ‘’ (12) at variable ‘’ (16).
MOSEK warning 62: The A matrix contains a large value of -1.7e+11 in constraint ‘’ (12) at variable ‘’ (17).
MOSEK warning 62: The A matrix contains a large value of 5.5e+10 in constraint ‘’ (12) at variable ‘’ (19).
Warning number 62 is disabled.
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 13
Cones : 12
Scalar variables : 38
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 - tries : 1 time : 0.00
Lin. dep. - tries : 1 time : 0.01
Lin. dep. - number : 0
Presolve terminated. Time: 0.11
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 13
Cones : 12
Scalar variables : 38
Matrix variables : 0
Integer variables : 0

Optimizer - solved problem : the primal
Optimizer - Constraints : 1
Optimizer - Cones : 12
Optimizer - Scalar variables : 38 conic : 36
Optimizer - Semi-definite variables: 0 scalarized : 0
Factor - setup time : 0.00 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 1 after factor : 1
Factor - dense dim. : 0 flops : 5.30e+01
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 3.3e+00 2.7e+11 5.8e+01 0.00e+00 5.468532531e+01 -2.415306005e+00 1.0e+00 0.16
1 7.8e-01 6.5e+10 2.8e+01 -1.00e+00 3.623288422e+01 -1.764751023e+01 2.4e-01 0.39
2 2.4e-01 2.0e+10 1.6e+01 -1.00e+00 -2.023743658e+01 -6.457482339e+01 7.3e-02 0.39
3 2.6e-02 2.1e+09 5.1e+00 -1.00e+00 -7.557907627e+02 -6.856443070e+02 7.8e-03 0.41
4 8.8e-05 7.4e+06 3.0e-01 -1.00e+00 -2.408314867e+05 -2.035103403e+05 2.7e-05 0.41
5 1.5e-10 1.3e+01 3.3e-04 -1.00e+00 -9.698071853e+10 -8.193276523e+10 4.7e-11 0.42
Optimizer terminated. Time: 0.63

Interior-point solution summary
Problem status : DUAL_INFEASIBLE
Solution status : DUAL_INFEASIBLE_CER
Primal. obj: -6.4447685340e+00 nrm: 6e+00 Viol. con: 4e+01 var: 6e+00 cones: 2e-11
Optimizer summary
Optimizer - time: 0.63
Interior-point - iterations : 5 time: 0.42
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00

Status: Infeasible
Optimal value (cvx_optval): -Inf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
And sedumi solver result is as follows.

Cones | Errors | Mov/Act | Centering Exp cone Poly cone | Status --------±--------------------------------±-------- 0/ 1 | 8.000e+00 5.373e+06 5.373e+06 | Failed 0/ 1 | 8.000e+00 3.909e+08 3.909e+08 | Failed 0/ 1 | 8.000e+00 6.995e+07 6.995e+07 | Failed

Status: Failed
Optimal value (cvx_optval): NaN

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Both solver failed to solve.

I wrote “It <Mosek> will also print our warnings and other useful diagnostic information” Did you look at the warnings from Mosek about large coefficients? Those large numbers arise from dividing by 1e-8.

Mosek reports the problem is infeasible. But it clearly is feasible. Even with the objective removed, Mosek still reports the problem as infeasible.

if the for loop assigning R1, R2, R3 is also removed, then the (feasibility) problem is solved (determined to be feasible). Apparently, even though R1, R2, R3 aren’t used anywhere in the problem if the objective is removed, R1, R2, R3 “stuff” is somehow getting into the problem sent from CVX to Mosek, and then the large numbers in the model are “screwing up” Mosek, resulting in an incorrect infeasibility determination. I have to defer to CVX developer @mcg why that is the case. Apparently CVX is not “smart” enough to see that R1, R2, R3 serve no value when the objective is removed. CVX is apparently including constraints based on them in the problem sent to the solver, even though the constraints are irrrelevant when the objective is omitted.

Finally, i note that R1, R2, R3 are overwritten every time though the for loop. So the objective function is just based on their values the last time through the loop (when i = K). Do you want to delcare R1, R2, R3 as expression holders http://cvxr.com/cvx/doc/basics.html#assignment-and-expression-holders and use R1(i) = …, etc., and then perhaps sum them in the objective?

Yes, thanks for pointing out the issue with R1,R2, and R3. I have to declare them as expression holders and have to sum them in the objective.

And Mosek is not solving the problem. I hope someone can figure out the issue.

K=4;
A =[ 24.0173 17.6637 78.0226 60.4402];
B =[ 156.1265 162.4801 102.1212 119.7036];
c=1.0000e-08*ones(1,4);

cvx_begin
variable tau0 nonnegative
expressions R1(K) R2(K) R3(K)
tau1=1-tau0;
for i=1:K
R1(i)=-rel_entr(tau1+ ((B(i)./c(i)) .* tau0 ), tau1+((B(i)+A(i))./c(i)) .* tau0 );

R2(i)=-rel_entr( (B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0, ...
(B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0+ (B(i) ./ c(i)).*tau0) ;

R3(i)=-rel_entr((B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i)) ).*tau0+ (B(i) ./ c(i)).*tau0,...
(B(i) ./ A(i)).*tau1 + (B(i)^2 ./( A(i).*c(i) ) ).*tau0) ;
end
maximize (sum(R1)+sum(R2)+sum(R3));
subject to
tau0<=1;
cvx_solver mosek
cvx_end

Yes, the issue is huge coefficients in your model. See if you can rescale so as to avoid that.

Mosek didn’t issue the warnings because it was bored and just looking for something to say.

1 Like

By decreasing c from 10^{-2} (where it solves nicely) to 10^{-8} you see from the log output that the number of iterations increases, the dual solution norm reported by Mosek grows, violations grow, the entries in the linear constraints grow and the expected value of \tau_0 is closer and closer to 0. So it is not surprising that at some point you hit the numerical accuracy issues.

Note that your problems ultimately contains subexpressions such as 10^{10}\tau_0 in various places. It is hard to expect it to be nicely solvable.

2 Likes

Yes. Thanks for your valuable suggestions. I could rescale the value of c and could run the program properly.