Stochastic Generative Automaton Optimization

Let $S = \lbrace s_1,\dots,s_N\rbrace$ denote a set of states and $Y$ an output alphabet with $\vert Y\vert=K$. From a state $s_i$ you transition to state $s_j$ via an edge labeled with $y\in Y$ and a probability $p_{y,ij}\in [0,1]$. This "stochastic automaton" that looks a bit like a Mealy machine has the following properties:

1. $\forall i.\sum_{y,j} p_{y,ij}=1$ (for all $s_i\in S$ the probabilities of outgoing transitions sum up to $1$)
2. for every $s\in S$ and $y\in Y$ there is at most one outgoing transition edge in $s$ labeled $y$

Starting from a given initial state $s^0$ you randomly take an outgoing transition edge with the probability given by its label. Every state transition $(s^t,s^{t+1})$ produces an output $y^{t+1}\in Y$. This is similar to a HMM, but the output that is produced is fixed for the transition edge taken, not randomly produced by the state itself. Property 2 ensures that, given an output sequence $(y^1,\dots,y^T)$ and the initial state $s^0$, you can reconstruct the state sequence $(s^0,s^1,\dots,s^T)$ unambiguously.

My problem: Given an initial state distribution $\pi^0$, what are the optimal $p_{y,ij}$ so that all $y\in Y$ appear equally often on average when running the automaton indefinitely? Note that the graph structure of the automaton, i.e. the state transition edges and their $y$-labels, is given and fixed. The only variables are all the non-zero $p_{y,ij}$.

My thoughts: My initial thought was that this is some kind of Markov Chain problem. Let $P_y=[p_{y,ij}] \in [0,1]^{N\times N}$ be the matrix of all transition probabilities that produce output $y$. Note that $P_y$ is a sparse matrix with at most one non-zero entry per row, due to property 2. Then $P=\sum_y P_y$ is a transition matrix. Let $\pi = (\pi^\top P)^\top$ be the stationary distribution of this matrix, then the probability of output $y$ being produced is $p_y:=\vert\pi^\top P_y\vert=\sum_i \pi_i \sum_j p_{y,ij}$. Now this is an optimization problem where I want all the $p_y$ to be the same, so I have to minimize the variance $$L(\theta)=\sum_y (\frac{1}{K}-p_y)^2 \longrightarrow \mathrm{min}$$ with $\theta=(P_y)_{y\in Y}$.

This is quite hard, though, because $p_y$ depends on the stationary distribution of a Markov Chain that changes whenever I change $p_{y,ij}$, so I cannot simply take the derivative $\partial L(\theta)/\partial p_{y,ij}$ for applying. Also I cannot solve this directly using gradient descent anyway because changing $p_{y,ij}$ would break property 1, therefore I'd have to solve this using unnormalized probabilities $\tilde{p}_{y,ij}$ as target parameters for optimization, which makes everything even more complex.

My questions:

• Am I thinking too complicated and is there a simpler model for my problem?
• If not, how can I approach to solve the optimization problem? Can you think of another way to solve this without using Markov chains and stationary distributions?

Thank you in advance! - smuecke