we define $\bm{\varXi}^r=|\det{\boldsymbol{\varTheta}}^r|({\boldsymbol{\varTheta}}^r)^{-{{\!\top\!}}}$. Problem  decouples with respect to ${\boldsymbol{M}}$ and ${\boldsymbol{\varTheta}}$, so we can work out their updates individually.

The update of ${\boldsymbol{M}}$ admits a simple closed form solution, which can be derived via checking the KKT conditions. We denote the dual variable corresponding to ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{M}}={\boldsymbol{\mathit{1}}}^{{\!\top\!}}$ as ${\boldsymbol{\mu}}\in{\mathbf{R}}^K$. Setting the gradient of the Lagrangian with respect to $M_{nk}$ equal to zero, we have $$M_{nk} = 
\sum_{\ell=1}^N\sum_{j=1}^{K}\left(\varOmega_{n\ell}\varPi_{n\ell kj}^r+\varOmega_{\ell n}\varPi_{\ell njk}^r\right) \bigg/ \mu_k$$ and ${\boldsymbol{\mu}}$ should be chosen so that the constraint ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{M}}={\boldsymbol{\mathit{1}}}^{{\!\top\!}}$ is satisfied, which amounts to a simple re-scaling.

The update of ${\boldsymbol{\varTheta}}$ is not as simple as a closed form expression, but it can still be obtained very efficiently. Noticing that the nonnegativity constraint is implicitly implied by the individual $\log$ functions in the loss function, we propose to solve it using Newton’s method with equality constraints [@boyd2004convex §10.2]. Although Newton’s method requires solving a linear system of equations with $K^2$ number of variables in each iteration, there is special structure we can exploit to reduce the per-iteration complexity down to $\mathcal{O}(K^3)$: The Hessian of the loss function of  is diagonal, and the linear equality constraints are highly structured; using block elimination [@boyd2004convex §10.4.2], we ultimately only need to 