How can I compute fast the minimum of a linear plus Kulback-Leibler on the unit simplex?

Question

Given $a, x^0 \in \mathbb{R}^n$ I wish to compute $$\min_{x \in \Delta_n} a^t x + \sum_{i=1}^n x_i\log(x_i/x^0_i) - x_i +x^0_i $$ where $\Delta_n$ is the unit simplex $\{x \in \mathbb{R}^n \mid \sum_{i=1}^n x_i \leq 1, x \geq 0\}$.

How can I solve this problem efficiently?

I tried using cvx directly, but this procedure takes a lot of time, and I need to do it a lot of times per iteration as it's a subproblem of my algorithm.

Since this is a very specific problem I thought about asking, is there a quick way of computing this? Another way to pose the problem or something? Is this already solved anywhere known?

for tldr'ers the question ends here, but I'm gonna post other thing I tried, but didn't work:

Another thing I tried was this: Since optimality conditions give us $$ 0 \in a_i + \log(x_i) - \log(x^0_i) + n_i$$ for every index $i \in \{1, \dots n\}$, where $n$ is a vector on the normal cone $N_{\Delta_n}(x)$, selecting a vector $n$ where all the components are equal a number say $\lambda$, this will give $$ x_i = \exp(-a_i - \lambda + \log(x^0_i)) = e^{-\lambda}e^{-a_i}x^0_i$$ and then imposing $\sum_{i=1}^n x_i = 1$, the number $e^{-\lambda}$ can be computed as $$e^{-\lambda} = \frac{1}{\sum_{i=1}^n e^{-a_i}x^0_i}$$ However this doesn't seem to be true for some reason, because when I come back to check the optimality conditions, it doesn't hold. Is this procedure wrong? Why? If it was right it would be awesome because it takes so much less time.

I'd appreciate any ideas

Answer 1

You are not saying what constitutes too long in your case or what the dimensions are, but with a standard nonlinear solver such as ipopt, a random problem with $n=1000$ is solved in 0.2 seconds on an old laptop. I guess you could reduce that by a factor of 10 if you manually code a solver for this particular problem.

Tested with this snippet of YALMIP code

n = 1000;
x0 = rand(n,1);
a = randn(n,1);
x = sdpvar(n,1);
p = entropy(sqrt(eps)+x)-sum(x.*log(x0)) + sum(-x + x0);
optimize([x>=0,sum(x)<=1],p,sdpsettings('solver','ipopt'))

$n=10000$ is solved in around 1 second.