I am attempting to solve an entropy maximization problem numerically using CVXPY. Even though the problem passes the DCP checks, I keep getting an infinite result, and problem statuses like infeasible and unbounded_inaccurate (depending on the values I choose for the parameters). I manage to solve the same problem using a nonlinear constrained optimizer like Alabama in R. Below is a small example that reproduces the problem.
import cvxpy as cvx
vals = array([ 750., 770., 790., 810., 830., 850., 870., 890.,
910., 930., 950., 970., 990., 1010., 1030., 1050.,
1070., 1090., 1110., 1130.])
n = size(freq)
z = cvx.Variable(2,n)
a = cvx.Parameter(sign="positive", value=989.)
b = cvx.Parameter(sign="positive", value=.1)
d = cvx.Parameter(sign="positive", value=10.)
obj = cvx.Maximize(cvx.sum_entries(cvx.entr(z)))
cons = []
cons += [ z >= 0., cvx.sum_entries(z) == 1, cvx.sum_entries(z, axis=0) * vals == a ]
for i in range(n):
cons += [ cvx.logistic(b*(vals[i] - a - d)) * z[1,i] == cvx.exp(b*(vals[i] - a - d)) * (z[0,i] + z[1,i]) ]
prob = cvx.Problem(obj, cons)
prob.solve(solver=cvx.SCS)
I don't quite understand why CVX is having a harder time to solve this problem than other algorithms that are not designed for convex programming problems. Have I overlooked something in the way I have written the constraints?
* EDIT *
Having received no answer so far, I will also ask this question on the CVXPY Google group. I will update this thread accordingly.
I managed to solve my problem. The solution was to store the numeric value of the logit distribution using Numpy functions and then use its components in the constraints:
qre = np.exp(b.value*(vals - a - d.value))/(1.+np.exp(b.value*(vals - a - d.value)))
...
cons += [ qre[i] * (z[0,i]+z[1,i]) == z[1,i] ]