Finding a mixture of 1st and 0'th order Markov models that is closest to an empirical distribution I am interested in finding the distribution "$p^*$" closest to an empirical distribution $\hat{p}$ where $p^*$ is a mixture of first and zeroth order Markov models. That is, I want to find $$
p^* = \arg\min_p \sum_{i,j} D\left(\,\hat{p}(j|i)\, \| \, p(j|i) \,\right)
$$
subject to the following constraints


*

*$p(j|i) = \alpha \cdot r(j) + (1 - \alpha) \cdot q(j|i)$

*$r(j) \ge 0$ for all $j$

*$\sum_j r(j) = 1$

*$q(j|i) \ge 0$ for all $i$ and $j$

*$\sum_j q(j|i) = 1$ for all $i$


where $\alpha$ is a mixture parameter in $[0,1]$ that is given and fixed. (The model is not identifiable if $\alpha$ is also a parameter to be learned.)
I know I can use EM here -- the missing data is the indicator variable for which Markov chain was invoked at a given time-step. But I am more interested in taking a direct optimization approach, particularly if this is a convex problem.
How might one derive the updates for for optimizing this model? Any assistance in this endeavor is greatly appreciated.
 A: You can solve this problem by fitting only the more complex model, then transforming the resulting parameter estimates $q^*$ to recover $q$ and $r$. No iterative optimization needed!
Please forgjve me if jn thjs answer J have mjxed up $j$ and $i$. Jt's iust very djffjcult for me to keep them strajght. 
Suppose your dataset consists of one observation of length $T$: $\{X_t\}_{n=0}^T$. (Below, I'll extend my answer to a case with multiple series.) I understand your problem to be that at each time $t$, there is a latent variable $z_t$ that is Bernoulli with known parameter $\alpha$. If $z_t=0$, then $x_{t}$ is drawn from $r$, and if $z_t=1$, then $x_{t}$ is drawn from $q(\cdot|x_{t-1})$. If $z$ is shared for all series, that's a very different situation, which this answer does not address. So, the contribution of epoch $t$ to the likelihood is this. 
$$P(x_{t}|x_{t-1})
 = \sum_z P(z_t, x_t |x_{t-1} )
 = \sum_z P(z_t)P( x_t | x_{t-1}, z_t) 
 = (1-\alpha) r(x_{t}) + \alpha q(x_{t}|x_{t-1})$$
This is much easier to look at if one defines $q^*(i|j)\equiv (1-\alpha) r(j) + \alpha q(j|i)$: it becomes just $q^*(x_{t}|x_{t-1})$. The complete-data likelihood (conditioning on $x_0$) is $\prod_{t=1}^{T} q^*(x_{t}|x_{t-1})$, and you could use the usual maximum likelihood estimates: $\hat q^*(j|i) \equiv \frac{N_{i\rightarrow j }}{N_{j}}$, where $N_{j}$ is the number of occurrences of $i$ and $N_{i\rightarrow j }$ is the number of $i$ to $j$ transitions. You can also add pseudocounts if you want to "smooth" things out; those could be interpreted as MAP estimates under suitable priors on $q^*$. (How that relates to priors on $q$ and $r$ would be interesting to investigate...) 
Define the estimates for $r$ and $q$ (I'll call them $\hat r$ and $\hat q$) to be any numbers satisfying $\alpha\hat q(j|i) + (1-\alpha)\hat r(j) = \hat q^*(j|i)$. If the state space has cardinality $J$, you now face a linear system with $J^2 + J$ unknowns ($\hat q$'s and $\hat r$'s). The number of equations is tricky to count because some of the sum-to-one constraints are redundant, but my calculation is:


*

*J^2 equations from the definition I gave

*J equations from the sum-to-one constraints on $\hat q$

*1 equation from the sum-to-one constraints on $\hat r$, which I am going to ignore because I suspect it's redundant.


If you vectorize $\hat q,\hat r$ as follows (dropping the hats, sorry):
$$v = [q(1|1), q(1|2),... q(1|J), j(1), q(2|1), ... j(2), ... q(J|1), ... r(J)]$$ 
and write the RHS as 
$$w = [q^*(1|1), q^*(1|2), ...q^*(1|J), 1, ..., q^*(J|1), q^*(J|2), ...q^*(J|J), 1]$$, 
then your system is $Av = w$, with $w$ known and $A$ known. $A$ is square and has a block diagonal structure with blocks of size $J+1$. Each block has the form 
\begin{bmatrix}
    \alpha & 0 & 0 & 1-\alpha \\
    0 & \alpha & 0 & 1-\alpha \\
    0 & 0 & \alpha & 1-\alpha \\
    1 & 1 & 1 & 0             \\
\end{bmatrix}
. 
This is invertible, as can be seen by subtracting, from the bottom row, $\alpha^{-1}$ times each of the top three, to yield an upper triangular block with a nonzero diagonal. This implies that the parameters $q$ and $r$ are identifiable, and fittingly, it only holds when $\alpha$ is not zero or 1. (Of course, identifying $r$ and $q$ is impossible if there is one of them we don't sample from.) 
Actually, those row operations also encode an algorithm to compute $r$ and $q$, which maybe I could have doped out more directly if I weren't such a matrix nut. The procedure is:
For all $j$:


*

*$a(j|i) \leftarrow \alpha^{-1}q^*(j|i)$ for all $i$  (rescale top $J$ rows)

*$b(j) \leftarrow 1 - \sum_i a(j|i)$ (subtract top J rows from bottom)

*$b(j) \leftarrow b(j) (\frac{J}{\alpha})^{-1}$ (rescale coeff in bottom row to $(1-\alpha)$ 

*$b(j) \leftarrow b(j)\alpha/J$ (rescale coeff to $(1-\alpha)$ in bottom row

*$a(j|i) \leftarrow a(j|i) - b(j)$ (wipe out top J of right-hand column)

*$q(j|i) \leftarrow a(j|i)/\alpha$ (convert most of diagonal to 1's)

*$r(j) \leftarrow b(j|i)/(1-\alpha)$ (convert bottom of diagonal to 1)


Let me know if you don't end up with positive probabilities that sum to 1. I suspect that if you specify too small of an $\alpha$, you'll have a bad time with this. You may want to start with a big $\alpha$ and inch it up slowly until something goes wrong. That way you can measure the degree to which your system can be simply described. 
