You are using Lagrangian multipliers incorrectly. For this model I recommend not using Lagrangian multipliers but simply the reparametrization of $K-1$ parameters $\boldsymbol{\theta} = \left(\theta_1, \cdots, \theta_{K-1}\right)$. Any inference for $\theta_K$ can be obtained by the invariance property of the MLEs since it is a function of these $K-1$ parameters.

So your gradient is correct, ignoring the last term; however, you must replace $\theta_K$ with $1 - \sum_{i=1}^{K-1} \theta_i$. Now you will have a linear system of $K-1$ equations and parameters, which can be solved. Let $\boldsymbol{1}_{K-1} \in \mathbb{R}^{K-1}$ denote a vector of 1's. You should arrive at a formula that looks like \begin{eqnarray*} \left(N_K Diag\left(\frac{1}{N_1}, \cdots, \frac{1}{N_{K-1}}\right) +\boldsymbol{1}_{K-1}\boldsymbol{1}_{K-1}^{\prime}\right)\widehat{\boldsymbol{\theta}} &=& \boldsymbol{1}_{K-1} \\ \widehat{\boldsymbol{\theta}} &=& \begin{pmatrix} \frac{N_1}{N_k} \\ \vdots \\ \frac{N_{K-1}}{N_K} \end{pmatrix}, \end{eqnarray*} where the second equality follows from the Sherman-Morrison formula and matrix algebra.

If you want to use Lagrangian multipliers then one does not need to change the parameter space. Let $\boldsymbol{\theta}^{\ast} = \left(\theta_1, \cdots, \theta_K \right)$. We wish to optimize the likelihood in the reduced domain such that $\sum_{i=1}^K \theta_i = 1$. One way to code the Lagrange multiplier and restraint function is to add the term $\lambda(\sum_{i=1}^K \theta_i -1)$ to the log-likelihood. Therefore, the log-likelihood with the Lagrange multiplier is \begin{eqnarray*} l(\boldsymbol{\theta}^{\ast}, \lambda) = \sum_{i=1}^K N_i \log \left(\theta_i\right) + \lambda(\sum_{i=1}^K \theta_i -1). \end{eqnarray*}

You are using Lagrangian multipliers incorrectly. For this model I recommend not using Lagrangian multipliers but simply the reparametrization of $K-1$ parameters $\boldsymbol{\theta} = \left(\theta_1, \cdots, \theta_{K-1}\right)$. Any inference for $\theta_K$ can be obtained by the invariance property of the MLEs since it is a function of these $K-1$ parameters.

So your gradient is correct, ignoring the last term; however, you must replace $\theta_K$ with $1 - \sum_{i=1}^{K-1} \theta_i$. Now you will have a linear system of $K-1$ equations and parameters, which can be solved. Let $\boldsymbol{1}_{K-1} \in \mathbb{R}^{K-1}$ denote a vector of 1's. You should arrive at a formula that looks like \begin{eqnarray*} \left(N_K Diag\left(\frac{1}{N_1}, \cdots, \frac{1}{N_{K-1}}\right) +\boldsymbol{1}_{K-1}\boldsymbol{1}_{K-1}^{\prime}\right)\widehat{\boldsymbol{\theta}} &=& \boldsymbol{1}_{K-1} \\ \widehat{\boldsymbol{\theta}} &=& \begin{pmatrix} \frac{N_1}{N_k} \\ \vdots \\ \frac{N_{K-1}}{N_K} \end{pmatrix}, \end{eqnarray*} where the second equality follows from the Sherman-Morrison formula and matrix algebra.

If you want to use Lagrangian multipliers then one does not need to change the parameter space. Let $\boldsymbol{\theta}^{\ast} = \left(\theta_1, \cdots, \theta_K \right)$. We wish to optimize the likelihood in the reduced domain such that $\sum_{i=1}^K \theta_i = 1$. One way to code the Lagrange multiplier and restraint function is to subtract the term $\lambda(\sum_{i=1}^K \theta_i -1)$ from the log-likelihood. Therefore, the log-likelihood with the Lagrange multiplier is \begin{eqnarray*} l(\boldsymbol{\theta}^{\ast}, \lambda) &=& \sum_{i=1}^K N_i \log \left(\theta_i\right) - \lambda(\sum_{i=1}^K \theta_i -1) \end{eqnarray*} Of interest are the following partial derivatives: \begin{eqnarray*} \frac{\partial l(\boldsymbol{\theta}^{\ast}, \lambda)}{\partial \theta_j} &=& \frac{N_i}{\theta_i} - \lambda \quad \mbox{for} \quad j=1,\cdots,K\\ \frac{\partial l(\boldsymbol{\theta}^{\ast}, \lambda)}{\partial \lambda} &=& 1- \sum_{i=1}^K \theta_i \end{eqnarray*} Hence summing over all $K$ partial derivatives in the first equation, after equating them to 0, one obtains \begin{eqnarray*} \lambda \sum_{i=1}^K \theta_i &=& N\\ \lambda &=& N, \end{eqnarray*} where the second equation holds by setting $\frac{\partial l(\boldsymbol{\theta}^{\ast}, \lambda)}{\partial \lambda}$ to 0. Now plug-in $\lambda$ back to the first partial derivative equations to obtain the MLEs.

You are using Lagrangian multipliers incorrectly. For this model I recommend not using Lagrangian multipliers but simply the reparametrization of $K-1$ parameters $\boldsymbol{\theta} = \left(\theta_1, \cdots, \theta_{K-1}\right)$. Any inference for $\theta_K$ can be obtained by the invariance property of the MLEs since it is a function of these $K-1$ parameters.

So your gradient is correct, ignoring the last term; however, you must replace $\theta_K$ with $1 - \sum_{i=1}^{K-1} \theta_i$. Now you will have a linear system of $K-1$ equations and parameters, which can be solved. Let $\boldsymbol{1}_{K-1} \in \mathbb{R}^{K-1}$ denote a vector of 1's. You should arrive at a formula that looks like \begin{eqnarray*} \left(N_K Diag\left(\frac{1}{N_1}, \cdots, \frac{1}{N_{K-1}}\right) +\boldsymbol{1}_{K-1}\boldsymbol{1}_{K-1}^{\prime}\right)\widehat{\boldsymbol{\theta}} &=& \boldsymbol{1}_{K-1} \\ \widehat{\boldsymbol{\theta}} &=& \begin{pmatrix} \frac{N_1}{N_k} \\ \vdots \\ \frac{N_{K-1}}{N_K} \end{pmatrix}, \end{eqnarray*} where the second equality follows from the Sherman-Morrison formula and matrix algebra.

You are using Lagrangian multipliers incorrectly. For this model I recommend not using Lagrangian multipliers but simply the reparametrization of $K-1$ parameters $\boldsymbol{\theta} = \left(\theta_1, \cdots, \theta_{K-1}\right)$. Any inference for $\theta_K$ can be obtained by the invariance property of the MLEs since it is a function of these $K-1$ parameters.

So your gradient is correct, ignoring the last term; however, you must replace $\theta_K$ with $1 - \sum_{i=1}^{K-1} \theta_i$. Now you will have a linear system of $K-1$ equations and parameters, which can be solved. Let $\boldsymbol{1}_{K-1} \in \mathbb{R}^{K-1}$ denote a vector of 1's. You should arrive at a formula that looks like \begin{eqnarray*} \left(N_K Diag\left(\frac{1}{N_1}, \cdots, \frac{1}{N_{K-1}}\right) +\boldsymbol{1}_{K-1}\boldsymbol{1}_{K-1}^{\prime}\right)\widehat{\boldsymbol{\theta}} &=& \boldsymbol{1}_{K-1} \\ \widehat{\boldsymbol{\theta}} &=& \begin{pmatrix} \frac{N_1}{N_k} \\ \vdots \\ \frac{N_{K-1}}{N_K} \end{pmatrix}, \end{eqnarray*} where the second equality follows from the Sherman-Morrison formula and matrix algebra.

If you want to use Lagrangian multipliers then one does not need to change the parameter space. Let $\boldsymbol{\theta}^{\ast} = \left(\theta_1, \cdots, \theta_K \right)$. We wish to optimize the likelihood in the reduced domain such that $\sum_{i=1}^K \theta_i = 1$. One way to code the Lagrange multiplier and restraint function is to add the term $\lambda(\sum_{i=1}^K \theta_i -1)$ to the log-likelihood. Therefore, the log-likelihood with the Lagrange multiplier is \begin{eqnarray*} l(\boldsymbol{\theta}^{\ast}, \lambda) = \sum_{i=1}^K N_i \log \left(\theta_i\right) + \lambda(\sum_{i=1}^K \theta_i -1). \end{eqnarray*}