solve a positive definite linear system with $K$ variables. Together with the quadratic convergence rate of Newton’s method, the complexity of updating ${\boldsymbol{\varTheta}}$ is $\mathcal{O}(K^3\log\log\frac{1}{\varepsilon})$, where $\varepsilon$ is the desired accuracy for the ${\boldsymbol{\varTheta}}$ update. Noticing that the complexity of a naive implementation of Newton’s method would be $\mathcal{O}(K^6\log\log\frac{1}{\varepsilon})$, the difference is big for moderately large $K$. The in-line implementation of this tailored Newton’s method <span style="font-variant:small-caps;">ThetaUpdate</span> and the detailed derivation can be found in the supplementary material.

$\lambda>0$ initialize ${\boldsymbol{M}}$ using [@huang2016nips] initialize ${\boldsymbol{\varTheta}}\leftarrow \frac{1}{K(K+1)}({\boldsymbol{I}}+{\boldsymbol{\mathit{1}}}\!{\boldsymbol{\mathit{1}}}^{{\!\top\!}})$ $\widetilde{{\boldsymbol{\varOmega}}}\leftarrow{\boldsymbol{\varOmega}}\big/{\boldsymbol{M}}{\boldsymbol{\varTheta}}{\boldsymbol{M}}^{{{\!\top\!}}}$ $\widetilde{{\boldsymbol{M}}}\leftarrow {\boldsymbol{M}}\ast \left( \widetilde{{\boldsymbol{\varOmega}}}{\boldsymbol{M}}{\boldsymbol{\varTheta}}^{{{\!\top\!}}} + \widetilde{{\boldsymbol{\varOmega}}}^{{\!\top\!}}{\boldsymbol{M}}{\boldsymbol{\varTheta}}\right)$ $\widetilde{{\boldsymbol{\varTheta}}}\leftarrow{\boldsymbol{M}}^{{{\!\top\!}}}\widetilde{{\boldsymbol{\varOmega}}}{\boldsymbol{M}}$ $\widetilde{{\boldsymbol{M}}}\leftarrow \widetilde{{\boldsymbol{M}}} \operatorname{Diag}({\boldsymbol{\mathit{1}}}^{{\!\top\!}}\widetilde{{\boldsymbol{M}}})^{-1}$ $\widetilde{{\boldsymbol{\varTheta}}}\leftarrow$ <span style="font-variant:small-caps;">ThetaUpdate</span> $({\boldsymbol{M}},{\boldsymbol{\varTheta}})\leftarrow$ Amijo line search between $({\boldsymbol{M}},{\boldsymbol{\varTheta}})$ and $(\widetilde{{\boldsymbol{M}}},\widetilde{{\boldsymbol{\varTheta}}})$ ${\boldsymbol{M}}$ and ${\boldsymbol{\varTheta}}$

The entire proposed algorithm to solve Problem  is summarized in Algorithm \[alg:hmm\_id\]. Notice that there is an additional line-search step to ensure decrease of the loss function. The constraint set of is convex, so the line-search step will not incur infeasibility. Computationally, we find that any operation that involves $\varPi_{n\ell kj}^r$ can be carried out succinctly by defining the intermediate matrix $\widetilde{{\boldsymbol{\varOmega}}}={\boldsymbol{\varOmega}}/{\boldsymbol{M}}{\boldsymbol{\varTheta}}{\boldsymbol{M}}^{{{\!\top\!}}}$, where “$/$” denotes element-wise division between two matrices of the same size. The per-iteration complexity of Algorithm \[alg:hmm\_id\] is completely dominated by the operations that involve computing 