How to adapt a linear time Newton-Raphson numerical method for an optimisation problem with positivity constraints?

Question

I would appreciate some assistance in understanding how I can adapt a linear time Newton-Raphson root finding algorithm for unconstrained optimisation, to solve a problem where I introduce positivity constraints on the arguments I am optimising with respect to.

Query.

I am seeking tutorials, textbook or canonical paper references on how (or whether) I can use Newton-Raphson to solve the following generically specified constrained optimisation problem:

$$\underset{r_1, r_2}{\max} L(r_1, r_2) \quad \text{s.t.} \quad r_1 > 0, r_2 >0.\tag{1}$$

Where $L(r_1, r_2)$ is a function of $r_1$ and $r_2$ and other terms that can be treated as constants.

---

Further context and additional problem structure.

As I understand, one would normally solve the above for $r_1$ and $r_2$ explicitly using Karush-Kuhn-Tucker conditions, because it is an optimisation problem with inequality constraints.

On why Newton-Raphson is used, and why I seek to adapt a linear-time Newton-Raphson method. This is best seen from the perspective of the unconstrained optimisation case, for which I have derived and verified a solution.

1. The reason for the use of Newton-Raphson in the unconstrained case is that the arguments $r_1$ and $r_2$ are coupled in a nonlinear fashion (inside a digamma function). Meaning that in the unconstrained case, while we can compute partial derivatives with respect to each $r_1$ and $r_2$ and set to $0$ for first order conditions, we cannot separate and solve for $r_1$ and $r_2$ explicitly. Therefore in the unconstrained case, the iterative updates are

$$\mathbf{r}^{(t+1)} = \mathbf{r}^{(t)} - \mathbf{H}(\mathbf{r}^{(t)})^{-1} \mathbf{g}(\mathbf{r}^{(t)}),$$

with $\mathbf{r} = (r_1, r_2)^T$, $\mathbf{H}^{-1} \in \mathbb{R}^{2 \times 2}$ is the inverse of the Hessian, $\mathbf{g} \in \mathbb{R}^2$ is the gradient.

2. The reason why the unconstrained case admits a linear time Newton Raphson numerical method as a solution is due to the fact that additional structure in the Hessian can be used to speed up its inversion. That is, the Hessian has the structure

$$\mathbf{H} = \text{diag}(\mathbf{h}) + \mathbf{1}z\mathbf{1}^T,$$

where $\mathbf{h} \in \mathbb{R}^2$ and $z \in \mathbb{R}$.

For further context, $L(r_1, r_2)$ are those parts of an evidence lower bound that contain terms in $r_1$ and $r_2$. The arguments $r_1$ and $r_2$ which I am optimising with respect to are parameters of a (variational) Dirichlet distribution, hence the need to enforce the positivity constraints. And the constrained optimisation problem $(1)$ is a subroutine inside a variational E-step, whereby maximisation of $L(r_1, r_2)$ is with respect to the free variational parameters $r_1$ and $r_2$ as part of a co-ordinate ascent procedure.

The linear time Newton Raphson method I am looking to adapt to the constrained case comes from:

Blei, D., Ng, A., Jordan, M. (2003). Latent Dirichlet Allocation. Journal of Machine Learning Research 3 (2003) pp993-1022.

Answer 1

If there is actually a maximum (as opposed to the function increasing towards the constrained boundary) you can do reasonably well with Newton's method and step-halving. Compute the update; if it violates the constraints, repeatedly halve the step size until it doesn't violate the constraints any more.

This will be slow if there is no maximum (and it can potentially get into an infinite loop if the constraints are $\leq$ rather than $<$ and the maximum is on the boundary)

Answer 2

The class of methods devised to solve this generic problem are known in the optimisation literature as projected Newton methods. The following is extracted from Bertsekas and Gafni (1983), with some minor additions to the notation.

In their simplest formulation, these are methods which solve simply constrained problems, which take the form

$$\min{J(x)} \quad \text{subject to} \quad x \geq 0, \tag{4}$$

where $x = [x^1, \dots, x^n] \in \mathbb{R}^n$, $J: \mathbb{R}^n \rightarrow \mathbb{R}$ is a continuously differentiable function, and the vector inequality $x \geq 0$ is component-wise, i.e. $x^i \geq 0$ for all $i = 1, \dots, n$.

The above problem admits an iterative solution of the form

$$x_{k+1} = [x_k - \alpha_k D_k \nabla J(x_k)]^+, \tag{5}$$

where $a_k$ is a positive scalar step size, $D_k$ is a positive definite symmetric matrix which is diagonal with respect to some of the co-ordinates of $x$, and $[\cdot]^+$ denotes projection (with respect to the standard norm) on the positive orthant. For avoidance of doubt on this latter operation, Bertsekas (1982) defines $[\cdot]^+$ as

$$[z]^+ = \begin{bmatrix} \max \{0, z^1 \} \\ \vdots \\ \max \{0, z^n \} \\ \end{bmatrix},$$

where $z \in \mathbb{R}^n$. The precise details of how $\alpha_k$ and $D_k$ are chosen; as well as convergence properties (in particular superlinear convergence in the case where $D_k$ is chosen on the basis of the second derivatives of $J$) can be found in the aforementioned references. These are:

Bertsekas, D. (1982). Projected Newton methods for optimization problems with simple constraints. SIAM J. Control and Optimization, 20(2), 221-246. https://doi.org/10.1137/0320018

Bertsekas, D., Gafni, E. (1983). Projected Newton methods and optimization of multicommodity flows. IEEE Transactions on Automatic Control, 28(12), 1090-1096. https://doi.org/10.1109/TAC.1983.1103183