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: