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)
Ad = np.random.randn(p, n)
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
        sparse.hstack([sparse.csc_matrix((p, n)), -Ad]),
        # 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