Numerical scheme for convex optimization

Question

Given $(e_n)_{-N\le n\le N}\in\mathbb R^{2N+1}$ and $-1<x<1$, solve

\begin{eqnarray} &&\max_{(q_n)_{-N\le n\le N}\in\mathbb R^{2N+1}_+}~ \sum_{n=-N}^N (e_n-\log(q_n))q_n \\ \mbox{s.t.} &&\quad \sum_{n=-N}^N q_n=1 \mbox{ and } \sum_{n=-N}^N nq_n/N=x. \end{eqnarray}

My question is whether we have (very efficient) packages (C++ or Python) to solve numerically this convex optimization problem? E.g. processing time less than 0.1 second with $N=100$. The motivation comes from my previous question, and I need to solve 40,000 optimization problems of this type to construct a discrete Markov martingale.

Any answers and comments are highly appreciated!

My thoughts : Introduce the Lagrangian multipliers and define

$$L(q_{-N},\ldots, q_N,\lambda_1,\lambda_2):=\sum_{n=-N}^N (e_n-\log(q_n))q_n+\lambda_1(\sum_{n=-N}^N q_n-1)+\lambda_2(\sum_{n=-N}^N nq_n/N-x).$$

Setting $\partial_{q_n}L=e_n-\log(q_n)-1+\lambda_1+n\lambda_2/N=0$, $\partial_{\lambda_1}L=\sum_{n=-N}^Nq_n-1=0$ and $\partial_{\lambda_2}L=\sum_{n=-N}^Nnq_n/N-x=0$, one has $q_n=\exp(e_n-1+\lambda_1+n\lambda_2/N)$ and thus

$$\sum_{n=-N}^N \exp(e_n-1+\lambda_1+n\lambda_2/N)=1\quad \mbox{and}\quad \sum_{n=-N}^N \exp(e_n-1+\lambda_1+n\lambda_2/N)n/N=x.$$

We obtain finally

$$\sum_{n=-N}^N \exp(e_n+n\lambda_2/N)(n/N-x)=0.$$

Let $z:=\exp(\lambda_2/N)$, one has

$$\sum_{n=-N}^N \exp(e_n)z^{n}(n/N-x)=0=\sum_{n=-N}^N \exp(e_n)z^{n+N}(n/N-x).$$

Define $A_k:=\exp(e_{k-N})(k/N-1-x)$ for $0\le k\le 2N+1$. Then solving the above opitmization problem is equivalent to search for the strictly positive root of the polynomial $\sum_{k=0}^{2N+1}A_kz^k=0$. Is there any quick solver for this?

Answer 1

You can solve this problem very quickly. From $\partial_{q_n} L=0$, we have $$q_n(\lambda_1,\lambda_2) = e^{e_n - 1 + \lambda_1 + \lambda_2 n/N}.$$

We know that $q_n$ has to sum to 1 (the first constraint), so we can eliminate $\lambda_1$:

$$q_n(\lambda_1,\lambda_2) = \frac{e^{e_n - 1 + \lambda_1 + \lambda_2 n/N}}{\sum_{n'} e^{e_{n'} - 1 + \lambda_1 + \lambda_2 n'/N}}= \frac{e^{e_n + \lambda_2 n/N}}{\sum_{n'} e^{e_{n'} + \lambda_2 n'/N}}\equiv q_n(\lambda_2).$$

This determines the solution vector $\{q_n\}_{n}$ up to a single coefficient $\lambda_2$, which will be determined by your second constraint. You can find $\lambda_2$ using a 1-dimensional root-finding algorithm, i.e., satisfy $f(\lambda_2)=0$ where $$f(\lambda_2) := \sum_n q_n(\lambda_2) n/N - x$$ You can use for instance Newton's method, which should take milliseconds.