I have a convex optimization problem

\begin{align*} \max_{P\, =\, (p_{ij})} &\sum_{ij} p_{ij} \big(a_{ij}-\log(p_{ij})\big)\\ \text{st}\quad &p_{ij}\in\mathcal{P}_{n, m} \; \text{and} \sum_{i=1}^n p_{ij} (b_i-1) = 0 \quad \forall j\\ \text{where,}\; &\mathcal{P}_{n, m} := \{(p_{ij}) \in \mathbb{R}_+^{n \times m}:\; \sum_{j} p_{ij} = 1,\forall i\} \end{align*}

This problem is related to A convex optimization problem. For that problem it was suggested to use proximal operator and it worked well but for this one the proximal map doesn't have a closed-form. How to go about this?

A bit of search shows that there are these generalized proximal mappings where one replaces squared Euclidean distance with other distance-like functions. The keywords here seem to be "entropic proximal mapping" and "mirror descent". Are there versions of splitting algorithms such as Douglas-Rachford for these generalized proximal operators?

  • $\begingroup$ Am I correct that we interpret $p_{ij} \log(p_{ij})$ to be equal to $0$ if $p_{ij} = 0$? $\endgroup$
    – littleO
    Commented Dec 8, 2016 at 3:12
  • $\begingroup$ @littleO indeed $\endgroup$
    – ie86
    Commented Dec 8, 2016 at 3:59
  • $\begingroup$ By the way, how large are $m$ and $n$ typically? $\endgroup$
    – littleO
    Commented Dec 8, 2016 at 19:03
  • $\begingroup$ @littleO usually m << n, n could be thousands, m is say ~10 $\endgroup$
    – ie86
    Commented Dec 8, 2016 at 19:28

1 Answer 1


(This is my second attempt at a solution. My initial attempt can be found below.)

Your optimization problem can be expressed as $$ \operatorname{minimize}_{p \in \mathbb R^{n \times m}} \quad F(p) + G(p) $$ where $$ F(p) = \sum_{i,j} p_{ij}(\log(p_{ij}) - a_{ij}) + I_\Omega(p) $$ and $I_\Omega$ is the indicator function of the set $$ \Omega = \{ p \in \mathbb R^{n \times m} \mid p_{ij} \geq 0 \, \forall \, i,j\} $$ (so $I_\Omega$ enforces nonnegativity on $p$) and $G$ is the indicator function of the set $$ S = \{ p \in \mathbb R^{n \times m} \mid \sum_j p_{ij} = 1 \, \forall \, i \quad \text{and}\quad \sum_{i=1}^n p_{ij}(b_i - 1) = 0 \, \forall j \}. $$ Evaluating the prox-operator of $G$ requires projecting onto $S$, which is a linear algebra problem with a standard closed-form solution.

We do not have a closed-form expression for the prox-operator of $F$, but $F$ is still a very simple function (for example, it is fully separable). Can we not evaluate the prox-operator of $F$ numerically to high precision very efficiently?

Evaluating the prox-operator of $F$ reduces to evaluating the prox-operator of the function $w:\mathbb R \to \mathbb R \cup \{ \infty \}$ defined by $$ w(x) = \begin{cases} x \log(x) & \quad \text{if } x \geq 0, \\ \infty & \quad \text{otherwise.} \end{cases} $$ (We interpret $x \log(x) = 0$ when $x = 0$.) To evaluate the prox-operator of $w$ at $\hat x$ (with parameter $t > 0$), we must solve \begin{align} \operatorname{minimize}_x & \quad x \log(x) + \frac{1}{2t} (x - \hat x)^2 \\ \text{subject to} & \quad x \geq 0. \end{align} If we visualize the graph of $ v(x) = x \log(x) + \frac{1}{2t} (x - \hat x)^2$, we see that $v$ is initially decreasing (as we move away from the origin, to the right) but then eventually increases to $+\infty$. So, $v$ has a minimizer in the interval $(0,\infty)$ which can be found by setting $v'(x) = 0$. This yields: $$ 1 + \log(x) + \frac{1}{t}(x - \hat x) = 0 \iff x + t \log(x) = \hat x - t. $$ The function $x + t \log(x)$ is strictly increasing, and it ranges from $-\infty$ to $+\infty$ as $x$ ranges from $0$ to $+\infty$. Thus, there is a single value of $x$ for which $x + t \log(x) = \hat x - t$. I think this value of $x$ can be found to high accuracy with Newton's method (and I think it won't take too many iterations). This allows us to evaluate the prox-operator of $F$ efficiently.

Alternatively (and even better), we can express the solution to $$ \frac{x}{t} + \log(x) = c $$ in terms of the Lambert W-function. Note that \begin{align} \frac{x}{t} + \log(x) = c & \iff \log(e^{x/t}) + \log(x) = \log(e^c) \\ &\iff x e^{x/t} = e^c \\ &\iff (x/t) e^{x/t} = \frac{e^c}{t}\\ &\iff x/t = W(e^c/t) \\ &\iff x = t W(e^c/t). \end{align} In Matlab, the Lambert W-function can be evaluated using lambertw.

Now that we have seen how to evaluate the prox-operators of both $F$ and $G$ efficiently, we are able to minimize $F(p) + G(p)$ using the Douglas-Rachford method (which solves convex problems of exactly this form).

Here is my initial attempt to solve the problem, which had a certain flaw (or issue, at least) which is discussed at the end.

The Chambolle-Pock algorithm solves convex optimization problems in the canonical form (or "graph form") minimize $f(x) + g(Cx)$, where the convex functions $f$ and $g$ are "simple" (meaning they have easy prox-operators) and $C$ is a matrix or a linear transformation. The Chambolle-Pock algorithm was extended by Condat (in this paper) to solve convex optimization problems of the form $$ \tag{$\spadesuit$} \operatorname{minimize}_x \quad h(x) + f(x) + g(Cx) $$ where $f,g$ and $C$ are as above and the convex function $h$ is differentiable with a Lipschitz continuous gradient. Your problem has the form $(\spadesuit)$ where $C$ is the identity mapping and \begin{align} h(p) = -\sum_{ij} \, p_{ij}(a_{ij} - \log(p_{ij})) \\ \end{align} and $f$ is the indicator function of $$ \{ p \in \mathbb R^{n \times m}_+ \mid \sum_j p_{ij} = 1 \, \forall i \} $$ and $g$ is the indicator function of $$ \{ p \in \mathbb R^{m \times n} \mid \sum_{i=1}^n p_{ij}(b_i - 1) = 0 \, \forall j \}. $$

NOTE: A problem (or at least an issue) with this approach is that $h$ does not have a Lipschitz continuous gradient and cannot be extended to a differentiable function $\mathbb R^{n \times m}$. That might mean this approach doesn't work.

I'd be interested to see what @MichaelGrant says about this problem.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.