# Maximum Entropy

This example demonstrates an instance of using the exponential cone. In this problem we want find the maximum entropy point inside a convex polytope, ie, to solve

$\begin{split}\begin{array}{ll} \mbox{maximize} & -\sum_i^n x_i \log x_i \\ \mbox{subhect to} & {\bf 1}^T x = 1 \\ & Ax - b \geq 0 \end{array}\end{split}$

over variable $$x \in \mathbf{R}^{n}$$, where $$A \in \mathbf{R}^{m \times n}$$ and $$b \in \mathbf{R}^m$$ are data. The problem has the following equivalent form,

$\begin{split}\begin{array}{ll} \mbox{minimize} & -{\bf 1}^T t \\ \mbox{subject to} & {\bf 1}^T x = 1 \\ & Ax - b \geq 0 \\ & \begin{bmatrix} t_i \\ x_i \\ 1 \end{bmatrix} \in \mathcal{K}_\mathrm{exp}, \quad i=1,\ldots,n, \end{array}\end{split}$

over variables $$x \in \mathbf{R}^{n}$$, $$t \in \mathbf{R}^{n}$$ and where $$\mathcal{K}_\mathrm{exp} \subset \mathbf{R}^3$$ denotes the exponential cone.

Python code to solve this is below.

import scs
import numpy as np
from scipy import sparse

# Generate problem data
np.random.seed(1)

# Matrix size parameters
n = 50  # Number of variables
p = 20  # Number of constraints

# Generate random problem data
tmp = np.random.rand(n)
tmp /= np.sum(tmp)
bd = 0.5 * Ad.dot(tmp) + 0.01 * np.random.rand(p)

# Build the A, b rows corresponding to the exponential cone
A_exp = sparse.lil_matrix((3 * n, 2 * n))
b_exp = np.zeros(3 * n)
for i in range(n):
A_exp[i * 3, i] = -1  # t
A_exp[i * 3 + 1, i + n] = -1  # x
b_exp[i * 3 + 2] = 1

A = sparse.vstack(
[
# zero cone
sparse.hstack([sparse.csc_matrix((1, n)), np.ones((1, n))]),
# positive cone
# exponential cones
A_exp,
],
format="csc",
)
b = np.hstack([1, -bd, b_exp])
c = np.hstack([-np.ones(n), np.zeros(n)])

# SCS data
data = dict(A=A, b=b, c=c)
# ep is exponential cone (primal), with n triples
cone = dict(z=1, l=p, ep=n)

# Setup workspace
solver = scs.SCS(
data,
cone,
)
sol = solver.solve()
x_scs = sol["x"][-n:]

# Verify solution with CVXPY
try:
import cvxpy as cp
except ModuleNotFoundError:
print("This example requires CVXPY installed to run.")
raise

x = cp.Variable(shape=n)
obj = cp.Maximize(cp.sum(cp.entr(x)))
constraints = [cp.sum(x) == 1, Ad @ x >= bd]
prob = cp.Problem(obj, constraints)
prob.solve(solver=cp.ECOS)
x_cvxpy = x.value

print(f"CVXPY optimal value is:", prob.value)
print(f"Solution norm difference: {np.linalg.norm(x_scs - x_cvxpy, np.inf)}")


After following the python install instructions, we can run the code yielding output:

------------------------------------------------------------------
SCS v3.1.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 100, constraints m: 171
cones: 	  z: primal zero / dual free vars: 1
l: linear vars: 20
e: exp vars: 150, dual exp vars: 0
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct
nnz(A): 1150, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
0| 1.11e+01  6.95e-01  5.59e+02 -3.06e+02  1.00e-01  9.02e-04
225| 8.91e-06  1.30e-05  3.24e-05 -3.89e+00  5.45e+00  1.43e-01
------------------------------------------------------------------
status:  solved
timings: total: 1.44e-01s = setup: 5.53e-04s + solve: 1.43e-01s
lin-sys: 1.06e-03s, cones: 1.42e-01s, accel: 4.26e-05s
------------------------------------------------------------------
objective = -3.888803
------------------------------------------------------------------
CVXPY optimal value is: 3.888784623785257
Solution norm difference: 2.241074541050464e-06