Problem with solution of maximum entropy problem (or request for documentation)

[gajomi]

I thought it would be nice to contribute a pair of examples for finding exponential family distributions (http://en.wikipedia.org/wiki/Exponential_family) with a given set of features and expectation values. The first example would consist of maximizing the relative entropy between the target distribution and a base measure under constraints n the expectation values. The second example would be the same thing, but where the maximization occurs on a much lower dimensional dual space (the Langrange multipliers that fix expectation values for a set of chosen features), which presumably is the easier problem. The point of the example would be to show that the two methods yield the same result.

Here is a bit of code that (at least initially), appeared to solve the first part of the problem for the case where the solution should be a Binomial distribution

using Gadfly
using Distributions
using Convex, SCS

#space and probalities
N = 32
X = 0:N
p = Variable(length(X))
probability_contraints = [0 <= p, p <= 1, sum(p) == 1];
#base measure
h = pdf(Binomial(N,1/2))
#feature expectation values
q = .3
μ =N*q
fs = [identity]
η = [μ]
C = hcat([map(f,X) for f in fs]...)'
feature_constraints = C * p == η
#definte optimization problem and solve
problem = maximize(entropy(p)+dot(p,log(h)),probability_contraints...,feature_constraints)
solve!(problem, SCSSolver(verbose=0))
#plot solution
plot(layer(x=X,y=p.value,Geom.point),
     layer(x=X,y=pdf(Binomial(N,q)),Geom.bar),Scale.x_discrete)

However, trying out a larger number of states (N=64) results in:

In both cases problem.solution.status is :Optimal. However, from the figure above it is clear that the solution is considerably far from the linear constraint. So, is there a bug (in convergence or constraint enforcement) at larger N values or do I just have the wrong idea about the range of final states that yield :Optimal convergence status? If the later is the case, is there some documentation one sensitivity parameters to ensure one is near the constraint subspace?

P.S.

This is my first time looking into Convex.jl and this whole thing about automatic convexity detection within such a rich modeling language seems very magical. Its great :)

[karanveerm]

Thanks for taking the time to write some examples! If you run SCS with verbose=1, you'll see a message which says:

Status: Solved/Inaccurate
Hit max_iters, solution may be inaccurate

So, we're hitting the maximum number of allowed iterations. If you instead use, say solve!(problem, SCSSolver(verbose=1, max_iters=15000)), you should see a graph that looks more like this:

I don't think Solved/Inaccurate should be interpreted as a status of :Optimal, since this is misleading. I'm going to open up another issue for that. Let us know if you run into other problems!

[mlubin]

@karanveerm, usually the :UserLimit status is used to indicate that the solver stopped because of the iteration limit.

[mlubin]

@gajomi, we're starting to collect example notebooks at https://github.com/JuliaOpt/juliaopt-notebooks, if you'd like to put together this example with a bit of explanation and LaTeX math, the contribution would be more than welcome.

[gajomi]

@karanveerm - Thanks. This fixes my problem. @mlubin - Yeah, I can put together a little notebook about all this. I should have some time next week.