
\section{Physical example of KL divergence as energy cost} \label{app:physical_systems}

\begin{figure}[!t]
    \centering
   
    \includegraphics[width=3in]{figures/KL_div_work_v4.pdf}
   
   
    \caption{(top) A gas molecule is found underneath the $i$-th piston with 
    probability proportional to volume $V_i$. (bottom) Inserting impermeable 
    partitions between the pistons, isothermally compressing to new volumes 
    $V'_i$, and then removing the partitions incurs expected work on the system 
    proportional to the KL divergence between the initial and final probability 
    distributions of the molecule's position.}
    \label{fig:KL_div_work}
\end{figure}



We offer an example of the KL divergence as the cost of sampling from a 
target distribution $p$ given a ``base'' distribution $q$ (see discussion in 
section \ref{sec:KL_div_cost}).

Figure \ref{fig:KL_div_work} presents a slight generalization of the Szilard 
information engine \cite{Szilard29}: molecules of an ideal gas inhabit the 
space formed by movable pistons indexed by $i \in \mathcal{X}$.  Let $V_i$ denote 
the volume beneath the $i$-th piston and $V \triangleq \sum_{i \in \mathcal{X}} V_i$ be 
the total volume.  A molecule of gas is equally likely to be found anywhere 
within the space beneath the pistons, corresponding to probability distribution 
$p \triangleq (V_i/V)_{i \in \mathcal{X}}$ on the pistons.  Now imagine we add 
impermeable partitions between the pistons (vertical dashed lines) and move the 
pistons to new positions (green lines) at constant temperature (perhaps the 
bottom of the box is in thermal contact with a heat reservoir).  The partitions 
prevent mixing between different pistons during compression; we remove them 
afterwards.  Let $V_i'$ be the new volume beneath the $i$-th piston, $V' \triangleq 
\sum_{i \in \mathcal{X}} V'_i$ be the new total volume, and $q \triangleq (V'_i/V)_{i 
\in \mathcal{X}}$ be the new piston probability distribution after this deformation.  

What is the work used to perform this deformation?  The work to isothermally 
compress an ideal gas is $k_\text{B} \tau \log(V_i/V'_i)$, where $k_\text{B}$ 
is Boltzmann's constant and $\tau$ is the temperature.  A gas molecule 
occupies the space beneath the $i$-th piston with probability $V_i/V$ before 
compression, so the expected work to move the pistons is
\begin{equation}\label{gaspistons}
    \sum_{i \in \mathcal{X}} \frac{V_i}{V} k_\text{B} \tau 
    \log\left(\frac{V_i}{V'_i}\right) = k_\text{B} \tau (D(p||q) + \log(V/V'))
\end{equation}
per molecule of gas.  If the new pistons positions are such that the total 
volume is unchanged, then $V' = V$ and the work is proportional to the KL 
divergence $D(p||q)$. We can imagine a sequence of such gas boxes and 
deformations, where the pre-deformation volumes $(V_{t,i})_{i \in \mathcal{X}}$ at 
time $t$ are determined by drawing a single molecule from volumes 
$(V_{t-1,i})_{i \in \mathcal{X}}$ at time $t-1$, forming a Markov chain with KL 
divergence control cost.

\section{Proof of Theorem 1} \label{app:proof}

\subsection{Part 1}
We wish to solve the following problem:
\begin{equation}
    P^* \triangleq \argmin_{P: \pi(P)= \pi^*} D(P || Q) \label{eq:def_optimization_problem}
\end{equation}

For ease of manipulation, we work with the empirical joint transition probability distribution $N = (N[i,j])_{i,j\in\mathcal{X}}$
\begin{equation}
	N[i,j](x_0^n) \triangleq \frac{1}{n}|\{t : (x_{t-1},x_t) = (i,j)\}|.
\end{equation}

We can solve (\ref{eq:def_optimization_problem}) by setting up the Lagrangian: 
\begin{align}
    \Lambda \triangleq & D(P||Q) + \alpha \left(\sum_{i,j \in \mathcal{X}} N_P[i,j] - 1\right) \nonumber \\
    &+ \sum_{i \in \mathcal{X}} \lambda_i \left(\sum_{j \in \mathcal{X}} N_P[i,j] - \pi^*[i]\right) \\
    & + \sum_{j \in \mathcal{X}} \eta_j \left(\sum_{i \in \mathcal{X}} N_P[i,j] - \pi^*[j]\right) \label{eq:def_Lagrangian}
\end{align}
where $\alpha$ is a Lagrange multiplier enforcing normalization of the joint transition probability distribution under $P$, $N_P$, and $(\lambda_i)_{i \in \mathcal{X}}$, $(\eta_i)_{i \in \mathcal{X}}$ are Lagrange multipliers enforcing the stationary distribution condition $\pi(P) = \pi^*$.  Our solution is a stationary point of the Lagrangian with respect to $N_P$:
\begin{align}
&\left.\frac{\partial \Lambda}{\partial N_P[i,j]} \right|_{N_P = N^*} = \left. \frac{\partial D(P||Q)}{\partial N_P[i,j]} + \alpha + \lambda_i + \eta_j \right|_{N_P = N^*} = 0 \nonumber \\
&\forall i,j \label{eq:Lagrangian_partial_derivative} 
\end{align}
Since $\pi_P = \pi^*$, then $N^*$, $A^*$ satisfy:
\begin{align}
    \left.\frac{\partial D(P||Q)}{\partial N_P[i,j]}\right|_{N_P = N^*} &= \left.\frac{\partial}{\partial N_P[i,j]}\left(D(N_P||N_Q) - D(\pi^*||\pi)\right)\right|_{N_P = N^*} \\
    &= \left.\frac{\partial}{\partial N_P[i,j]}D(N_P||N_Q)\right|_{N_P = N^*} \\
    &= \left.1 + \log\left(\frac{N_P[i,j]}{N_Q[i,j]}\right)\right|_{N_P = N^*} \label{eq:Lagrangian_partial_derivative_expression}
\end{align}
Now using (\ref{eq:def_Lagrangian}), (\ref{eq:Lagrangian_partial_derivative}), and (\ref{eq:Lagrangian_partial_derivative_expression}) we find:
\begin{align}
    N^*[i,j] = N_Q[i,j]e^{-1-\alpha-\lambda_i-\eta_j}
\end{align}
Equivalently
\begin{equation}
    N^* = \text{diag}((e^{-\lambda_i})_{i \in \mathcal{X}}) \ N_Q \ \text{diag}((e^{-\eta_i})_{i \in \mathcal{X}}) e^{-1-\alpha}
\end{equation}
Now using the condition $\pi(P) = \pi^*$:
\begin{align}
    \pi^*[i] &= \sum_{j \in \mathcal{X}} N^*[i,j] \\
    &= e^{-1-\alpha-\lambda_i}\sum_{j \in \mathcal{X}} N_Q[i,j] e^{-\eta_j}
\end{align}
Therefore
\begin{equation}
    \lambda_i = \log\left(\sum_{j \in \mathcal{X}} N_Q[i,j] e^{-\eta_j}\right) - \log(\pi^*[i]) - 1 - \alpha \label{eq:recursive_lambda}
\end{equation}
Analogously
\begin{equation}
    \eta_i = \log\left(\sum_{j \in \mathcal{X}} N_Q[j,i] e^{-\lambda_j}\right) - \log(\pi^*[i]) - 1 - \alpha \label{eq:recursive_eta}
\end{equation}
Now using the normalization condition
\begin{equation}
    1 = \sum_{i,j} N^*[i,j]
\end{equation}
we find
\begin{equation}
    \alpha = \log\left(\sum_{i,j \in \mathcal{X}} N_Q[i,j]e^{-\lambda_i - \eta_j}\right) - 1 \label{eq:recursive_alpha}
\end{equation}

\subsection{Part 2}
In the continuous case, we work with rate matrices instead of probability transition matrices. We wish to solve the following problem:

\begin{equation}
\begin{aligned}
& \underset{\B{P}: \pi{(\B{P})} = \pi^*}{\text{min}}
& & D(\B{P}||\B{Q}) \\
\end{aligned}
\end{equation}
The Lagrangian is $D(\B{P}||\B{Q})+\sum_i \lambda_i \sum_j \B{P}_{ij} + \sum_j \eta_j \sum_i \pi_i \B{P}_{ij}.$ Differentiating it we get the conditions 
\begin{equation}
      \B{P}_{ij}^* = \begin{cases} \label{eq:rec} 
	\B{Q}_{ij} e^{\eta_i-\eta_j} & i \neq j \\
	-\sum_{j:j\neq i} \B{Q}_{ij} e^{\eta_i-\eta_j} & i = j \\
   \end{cases}
\end{equation}


\section{Discussion}\label{sec:discussion}

This work derives the minimum power required to maintain a target stationary 
distribution given uncontrolled Markov chain dynamics in both discrete and 
continuous time.  We relate KL divergence-like penalties from control theory 
\cite{Todorov06,Todorov09} to the power used to control a Markov process, using 
muscular molecular motors as a guiding example.  The problem of minimizing a KL 
divergence subject to a constrained stationary distribution is familiar from 
large deviations theory 
\cite{DemboZeitouniBook,BaldiPiccioni98,CsiszarCoverChoi87}; the novelty of our 
work is in relating these large deviations results to the thermodynamics of 
``holding'' a distribution, and in computing the minimum-cost chains in some 
important examples: the birth-and-death process in continuous time (a toy model 
for a muscular molecular motor) and two-state chains in discrete time.


To the best of our knowledge, this is the first time a lower bound on average 
power consumption has been studied in detail for control of the stationary 
distribution.  \cite{HorowitzEngland17} study a related quantity, the minimum 
entropy production rate associated with adding edges (allowing control to 
increase but not decrease transition rates) to a continuous time Markov chain, 
but their notion of cost has the interpretation of energy efficiency, as 
opposed to ours, which is to be interpreted as total energy usage.  
Unsurprisingly, different notions of cost lead to different optimal controlled 
chains: the optimal controlled chain of \cite{HorowitzEngland17} depends on the 
underlying uncontrolled chain only in the requirement that it be much faster, 
while our minimum-cost controlled chain is a function of the uncontrolled 
chain; this function is easy to compute (\ref{eq:reversible_cont_time}) in the 
case of a continuous time, reversible uncontrolled chain, an important case in 
modeling biological processes.


\section{Introduction} \label{sec:introduction}

Let $(X_1,X_2,\ldots)$ be a sequence of random variables forming a first-order 
Markov chain on a finite set $\mathcal{X}$ with
transition probability 
matrix $Q$, where $Q_{ij} = \mathbb{P}_Q(X_{t+1}=j|X_t=i)$
for all $t$ and $i,j\in\mathcal{X}$
We think of $Q$ as the ``uncontrolled'' or ``base'' chain.
If $Q$ is irreducible and aperiodic, then there exists a unique stationary 
distribution $\pi$, viewed as a row vector:
\begin{equation}
   
    \pi Q = \pi
   
   
   
\end{equation}
Let the initial state $X_1$ have distribution $\mu$.  Then, denoting $X_1^T 
\triangleq (X_1, \ldots X_T)$, the distribution of $X_1^T$ is:
\begin{align}
    p_{Q,\mu}(x_1^T) &\triangleq \mathbb{P}_{Q,\mu}(X_1^T = x_1^T) = \mu_{x_1} 
    \prod_{t=2}^T Q_{x_{t-1},x_t} \label{eq:def_MC_distribution}
\end{align}




Let $\pi^*$ be some probability distribution on $\mathcal{X}$ viewed as a row 
vector.  We study the nearest Markov chain transition matrix $P$ to $Q$ having 
$\pi^*$ as its stationary distribution:
\begin{align}
   
   
    P^* &= \argmin_{P: \pi^* P = \pi^*} D(P||Q) 
    \label{eq:optimization_problem_in_introduction}
\end{align}
where $D(P||Q)$ is the Kullback-Leibler (KL) divergence rate between Markov 
chains with transition matrices $P$ and $Q$ 
\cite{CsiszarCoverChoi87}~\footnote{This expression is independent of $\mu = 
\mathbb{P}(X_1)$ for aperiodic, irreducible $P$.}:
\begin{equation}
   
   
   
   
   
    D(P||Q) \triangleq \lim_{T \rightarrow \infty} \frac{1}{T} \mathbb{E}_P 
    \log\left(\frac{p_P(X_0^T)}{p_Q(X_0^T)}\right) 
    \label{eq:def_KL_div_MC_discrete_time_in_introduction}
\end{equation}
where we note that $D(P||Q)$ is independent of the initial distribution $\mu$.

We think of $P$ as the ``controlled'' or ``driven'' chain and $D(P||Q)$ as the 
cost of control per unit time -- the power.
We further consider the analogous question for the case of a continuous-time 
Markov chain $(X_t)_{t \in [0,\infty)}$.


This setting is inspired by the following thought experiment due to Feynman 
\cite{FeynmanLectures}: a person holds a heavy bag above the floor for an hour 
and gets tired.  The net work done on the object is zero\footnote{Since work is 
the product of force and distance displaced, the latter of which is zero.}, so 
why does she get tired?  A table could hold the same bag indefinitely without 
an energy source, and so could the person if she were frozen solid, 
the bag hanging on her stiff, lifeless limb.
The latter observation implicates the microscopic dynamics of muscles as key to 
this question.
A toy model for the motion of striated muscle fibers -- see 
\cite{Huxley57,Qian00} for an extended discussion -- is of a random walk in a 
periodic energy potential\footnote{The myosin protein joining the fibers 
together is doing the random walking, the energy potential having 
periodically-spaced minima corresponding to 
discrete steps along the fiber.}.  Attachment of the heavy bag pulls on the 
muscle fibers, biasing the random walk in the direction of gravity, tilting the 
energy potential.  The person must use chemical energy\footnote{The hydrolysis 
of ATP molecules.} to de-bias the random walk in such a way that the bag is 
held at the desired height above the floor.  
If the person is frozen solid,
then the underlying random walk stops\footnote{Or slows down a lot, as lowering 
the temperature reduces the transition rate between the potential's energy 
minima.}, and so chemical energy is no longer used to hold the bag.
We ask: what is a lower bound on the power to hold the bag?

We recast the above story as optimization problem 
(\ref{eq:optimization_problem_in_introduction}).
The state space $\mathcal{X}$ corresponds to the possible configurations of the 
system (the position of myosin along a fiber and its internal state).
The uncontrolled Markov chain $Q$ corresponds to the underlying fluctuations of 
the myosin molecule along a filament and the controlled chain $P$ corresponds 
to chemically driving the system.  

The control goal is macroscopic: the net force the person exerts on the bag is 
the sum of the forces due to each microscopic subsystem (myosin protein).  We 
do not get to control each subsystem separately, but can control them all in 
the same way, so that each subsystem corresponds to a trajectory drawn from 
Markov chain $P$ independently of other subsystems.  This notion of macroscopic 
control is reflected in the KL divergence cost 
(\ref{eq:def_KL_div_MC_discrete_time_in_introduction}), which is stated in 
terms of the probability distribution (\ref{eq:def_MC_distribution}) over 
microscopic trajectories, rather than in terms of a single trajectory.  

Our choice of the Kullback-Leibler divergence as the control cost function is 
motivated by this quantity's appearance in thermodynamics as proportional to 
the free energy difference from the equilibrium distribution over trajectories, 
which in turn lower bounds the work to prepare a non-equilibrium distribution 
over trajectories.
KL divergence control (the microscopic per-trajectory setting) was introduced 
in the reinforcement learning literature by Todorov \cite{Todorov06,Todorov09} 
and has connections to data compression; we discuss this and other motivations 
for our work in section \ref{sec:KL_div_cost}. The problem of maintaining a 
target nonequilibrium distribution has been studied recently by 
\cite{HorowitzZhouEngland17,HorowitzEngland17}, using a different notion of 
cost.  We discuss our work's relation to prior works in section \ref{subsec: 
comparison}.

Minimizing the KL divergence with respect to the first argument --  computing 
the I-projection -- 
connects this setting to large deviations theory.  Indeed, most of the 
minimum-cost controlled chains we compute first appeared in the computation of 
rate functions for large deviations of the empirical measure of both discrete- 
and continuous-time Markov chains 
\cite{DemboZeitouniBook,BaldiPiccioni98,CsiszarCoverChoi87}.  The novelty of 
our work lies in relating these results to the minimum-power control setting -- 
showing the minimized KL divergence to be the minimum power to ``hold'' a 
nonequilibrium distribution; in aggregating related problem statements -- in 
continuous and discrete time and for reversible base chains; in computing some 
of these minimizations more explicitly than previously available; and in 
computing these minimizations for a few common examples like the 
birth-and-death chain and the two-state chain.

This work is organized as follows.  Section \ref{sec:KL_div_cost} motivates the 
use of KL divergence as energy cost by drawing from information theory and 
optimization settings and contains definitions of this cost function in 
discrete and continuous time. Section \ref{sec:minimum_power} shows how this 
energy cost of holding a given target distribution may be analytically 
minimized.  Section \ref{sec:examples} contains several examples that apply our 
theory to calculate the minimum-power controlled chain, including a 
birth-and-death chain which serves as a toy model of the muscular fiber, 
addressing the motivating question of Feynman. We conclude with a summary and 
outlook in section \ref{sec:discussion}.

We release code for computing the minimum-cost chains in this work at 
\href{https://github.com/dmitrip/controlledMC}{https://github.com/dmitrip/controlledMC}.

\section{Kullback-Leibler divergence rate as the cost of control of Markov 
chains} \label{sec:KL_div_cost}

We motivate the KL divergence rate between Markov chains in both discrete and 
continuous time as the cost function lower bounding the power in the 
bag-holding thought experiment.  We present a thermodynamics perspective in 
subsection \ref{subsec:thermo} and an equivalent perspective due to Todorov 
\cite{Todorov06,Todorov09} of a Markov decision process with log ratio cost 
function in subsection \ref{subsec:log_loss_action_cost}.  We summarize known 
expressions for the KL divergence between Markov chains in discrete and 
continuous time in subsection \ref{sec:MC_discrete_time} and 
\ref{sec:MC_continuous_time} respectively. Finally, \ref{subsec: comparison} 
places this work in context with related work. 

\subsection{KL divergence in thermodynamics}\label{subsec:thermo}

We summarize briefly the appearance of the KL divergence in measuring work in 
statistical mechanics. 
Below, let $D(p||q) = \sum_{i \in \mathcal{X}} p(i) \log(p(i)/q(i))$ denote the 
Kullback-Leibler (KL) divergence between distributions $p$ and $q$ on finite 
set $\mathcal{X}$ and let $H(p) = \sum_{i \in \mathcal{X}} p(i) \log(1/p(i))$ denote the 
entropy of distribution $p$. 
Let \begin{equation}
    q = \left(\frac{1}{Z}e^{-\beta U(i)}\right)_{i \in \mathcal{X}} 
   
    \label{eq:Boltzmann_distribution}
\end{equation}
be the Boltzmann distribution on $\mathcal{X}$, where $U(\cdot)$
is the energy function (sometimes called the energy potential or internal 
energy), $\beta$ is the inverse temperature, and $Z = \sum_{i \in \mathcal{X}} e^{-\beta U(i)}$ is the partition function.
Denote the free energy of distribution $p$ by $F(p)$ 
\cite{MezardMontanariBook}:
\begin{align}
   
   
   
   
    F(p) &\triangleq \mathbb{E}_p(U(X)) - \frac{1}{\beta} H(p) \\
   
    &= \frac{1}{\beta} \mathbb{E}_p \log 
    \left(\frac{p_X}{q_X}\right) - \frac{1}{\beta} \log(Z)\\
   
    &= \frac{1}{\beta} D(p||q) - \frac{1}{\beta} \log(Z)
\end{align}
where the expectation $\mathbb{E}_p(\cdot)$ is over random variable $X$ with 
distribution $p$.  Then $F(\cdot)$ is minimized at equilibrium $p=q$, so that 
$F(q) = -\frac{1}{\beta}\log(Z)$.

In thermodynamics \cite{KittelKroemerBook,LandauLifshitz} the work $W$ to 
prepare the ``controlled'' distribution $p$ starting from the ``base'', 
equilibrium distribution $q$ (also known as the work on the system) is at least 
the free energy difference:
\begin{equation}
   
   
   
    W \geq F(p) - F(q) = \frac{1}{\beta} D(p||q) 
    \label{eq:work_lower_bound_free_energy}
\end{equation}
As is customary in this setting, there is in the background a notion of a 
stochastic process transforming initial states into final states, and the work 
$W$ in (\ref{eq:work_lower_bound_free_energy}) is an average over realizations 
of this process.  In Appendix (\ref{app:physical_systems}) we provide a 
physical example in the spirit of the Szilard's engine thought experiment in 
thermodynamics, for which the KL divergence does emerge as the work done. 
In the bag-holding thought experiment, there is a large collection of 
independent myosin systems, and the total work is the sum of the works on each 
system.  We imagine the number of microscopic systems to be large enough that
fluctuations about the average work $W$ per subsystem are small, so it is this 
average work that's our object of study.

The KL divergence cost is familiar in data compression, where the ``energy'' of 
symbol $X$ drawn from distribution $p$ is $U(X) = -\log(p_X)$.  Any compression 
scheme must use at least $\mathbb{E}_p(U(X)) = H(p)$ bits  to encode a sample from 
$p$ \cite{Shannon1948,CoverThomasBook} on average over draws from distribution 
$p$.  If we use a compression scheme that instead uses $-\log(q_X)$ bits to 
encode symbol $X$ -- a mismatched code -- then we would pay $D(p||q)$ extra 
bits per symbol on average.  Section \ref{sec:MC_discrete_time} contains 
analogous remarks for compressing samples drawn from Markov chain 
distributions.

\subsection{Markov chains in discrete time} \label{sec:MC_discrete_time}

In the Markov chain control setting, we apply the preceding picture with 
alphabet $\mathcal{X}^{T}$ (trajectories of length $T$) and with Markov chain 
distributions on $\mathcal{X}^T$ with a desired marginal distribution $\pi^*$.  We 
consider the continuous time setting in section \ref{sec:MC_continuous_time}.

A discrete time Markov chain distribution $p_{Q,\mu}^{(T)}$ on the set 
$\mathcal{X}^{T}$ is the Boltzmann distribution with energy function 
$U_{Q,\mu}^{(T)}$ parametrized by the stochastic transition matrix $Q = 
(Q_{ij})_{i,j\in\mathcal{X}}$ and the initial distribution vector $\mu = (\mu)_{i 
\in \mathcal{X}}$ (obtained by taking the logarithm of the rightmost quantity in 
(\ref{eq:def_MC_distribution})):  
\begin{equation}
   
    U_{Q,\mu}^{(T)}(x_1^T) = -\log(\mu_{x_1}) - \sum_{t=2}^{T} 
    \log(Q_{x_{t-1},x_t})
   
   
   
   
\end{equation}
where $x_1^T \triangleq (x_1, \ldots, x_T)$.  Then $p_{Q,\mu}^{(T)}(x_1^T) = e^{-U_{Q,\mu}^{(T)}(x_1^T)}$, where $Z = \beta = 
1$ (\ref{eq:Boltzmann_distribution}).  Given another transition matrix $P$ and 
initial distribution $\nu$, the work to prepare Markov chain distribution 
$p_{P,\nu}^{(T)}$ starting from the Boltzmann distribution $p_{Q,\mu}^{(T)}$ is 
lower bounded by the free energy difference $F(p_{P,\nu}^{(T)}) - 
F(p_{Q,\mu}^{(T)})$ (\ref{eq:work_lower_bound_free_energy}).  In the limit $T 
\rightarrow \infty$, the work per time step -- the power -- is lower bounded by
\begin{align}
    \lim_{T \rightarrow \infty} \frac{1}{T} W(p_{P,\nu}^{(T)}) &\geq \lim_{T \rightarrow 
    \infty} \frac{1}{T} \left(F(p_{P,\nu}^{(T)}) - F(p_{Q,\mu}^{(T)})\right) 
    \label{eq:power_lower_bound_free_energy_mc} \\
    &= \lim_{T \rightarrow \infty} \frac{1}{T} D(p_{P,\nu}^{(T)}||p_{Q,\mu}^{(T)}) \\
   
   
    &= \lim_{T \rightarrow \infty} \frac{1}{T} \mathbb{E}_{p_{P,\nu}^{(T)}} 
    \log\left(\frac{p_{P,\nu}^{(T)}(X_1,\ldots,X_T)}{p_{Q,\nu}^{(T)}(X_1,\ldots,X_T)}\right) 
    \label{eq:KL_div_MC_discrete_time_expression} \\
    &= \sum_{i \in \mathcal{X}} \pi(P)_i \sum_{j \in \mathcal{X}} 
    P_{ij} \log\left(\frac{P_{ij}}{Q_{ij}}\right) \\
    &\triangleq D(P||Q) \label{eq:def_KL_div_MC_discrete_time_1}
\end{align}
where $\pi(P)$ is the stationary distribution of transition matrix $P$ and the 
last equality defines the KL divergence rate \cite{CsiszarCoverChoi87} between 
Markov chains with transition matrices $P$ and $Q$.

In data compression, a sample $X_1^T \sim p_{Q,\mu}^{(T)}$ can be compressed on 
average to at least $T H(Q)$ bits \cite{CoverThomasBook}, where $H(Q)$ is the 
entropy rate of Markov chain with transition matrix $Q$:
\begin{equation}
    H(Q) = -\sum_{i \in \mathcal{X}} \pi(Q)_i \sum_{j \in \mathcal{X}} Q_{ij} 
    \log(Q_{ij})
\end{equation}
Encoding samples from distribution $p_{Q,\mu}^{(T)}$ with respect to a 
mismatched code based on distribution $p_{P,\nu}^{(T)}$ incurs an average cost 
per unit time of at least $D(P||Q)$ extra bits.

Maximizing the lower bound on power (\ref{eq:power_lower_bound_free_energy_mc}) 
over transition matrices $P$ with target stationary distribution $\pi^*$ yields 
the optimization problem (\ref{eq:optimization_problem_in_introduction}):
\begin{equation}
    P^* = \argmin_{P: \pi^* P = \pi^*} D(P||Q) 
    \label{eq:optimization_problem_in_MC}
\end{equation}

\subsection{Log loss action cost} \label{subsec:log_loss_action_cost}

Another path to the same optimization problem 
(\ref{eq:optimization_problem_in_MC}) (KL divergence as a lower bound on the 
work to maintain a nonequilibrium distribution) is in terms of a Markov 
decision process with log loss action cost, a setting introduced by Todorov 
\cite{Todorov06,Todorov09}.  Let $Q$ be the uncontrolled chain and let $P$ be 
the controlled chain.  Let $c(i, j, P, Q)$ be the \textit{microscopic cost} 
paid when a transition is made from $X_t = i$ to $X_{t+1} = j$ when the 
controller chooses transition probability matrix $P$.  KL divergence control 
amounts to using the log likelihood ratio:
\begin{equation}
   
    c(i, j, P, Q) \triangleq \log\left(\frac{P_{ij}}{Q_{ij}}\right)
\end{equation}
If we view the rows of $P$ and $Q$ as Boltzmann distributions with different 
energy potentials -- that is, if we choose energy functions $E_i(\cdot)$, 
$E'_i(\cdot)$ such that $P_{ij} = e^{-E_i(j)}$ and $Q_{ij} = e^{-E'_i(j)}$ -- 
then the microscopic cost is the difference in energies: $c(i,j,P,Q) = E'_i(j) 
- E_i(j)$.


If $X_t \sim \mu$, then let the \textit{cost} $D_\mu(P||Q)$
be the expected microscopic cost (the average cost paid per microscopic 
system):
\begin{align}
   
   
    D_\mu(P||Q) &\triangleq \mathbb{E}_{\mu, P}(c(X_t, X_{t+1}, P, Q)) \\
   
   
    &= \sum_{i \in \mathcal{X}} \mu_i \sum_{j \in \mathcal{X}} P_{ij}\log 
    \left(\frac{P_{ij}}{Q_{ij}}\right) \label{eq:cost_disc_time_KL_div}
\end{align}
Thus 
$D_\mu(P||Q)$ is a $\mu$-weighted KL divergence between the rows of transition 
matrices $P$ and $Q$.  We are interested in macroscopic control of $\mu_t$, 
rather than microscopic control of $X_t$, so our setup differs from the setting 
introduced by \cite{Todorov06,Todorov09}: we average the control cost over 
$\mu_t$, so there is no randomness in our setting.  Finally for irreducible, 
aperiodic transition matrix $P$, we have the identity
\begin{equation}
   
   
   
    D_{\pi(P)}(P||Q) = D(P||Q) \label{eq:def_KL_div_MC_discrete_time}
\end{equation}
where $\pi(P)$ is the stationary distribution of $P$.  Minimizing the cost 
$D_{\pi^*}(P||Q)$ with respect to $P$ such that $\pi^* P = \pi^*$ is 
optimization problem (\ref{eq:optimization_problem_in_MC}).

\subsection{Markov chains in continuous time} \label{sec:MC_continuous_time}

The setup of section \ref{sec:MC_discrete_time} has a natural counterpart for 
continuous-time Markov chains.  Let $\bar{Q} = (\bar{Q}_{ij})_{i,j \in \mathcal{X}}$ 
denote the transition rate matrix of the uncontrolled continuous-time Markov 
chain $(X_t)_{t \in \mathbb{R}_{\geq 0}}$, where henceforth the overbar notation 
corresponds to rate matrices.  Let $\bar{P}$ be the controlled rate matrix and 
let $X_t \sim \mu_t$.  Then
\begin{equation}
    \frac{d}{dt} \mu_t = \left(\frac{d}{dt} \mathbb{P}_{\bar{P}}(X_t = i)\right)_{i 
    \in \mathcal{X}} = \mu_t \bar{P} \Rightarrow \mu_t = \mu_0 e^{\bar{P} t} 
    \label{eq:controlled_state_continuous_time}
\end{equation}
where $e^{\bar{P}t}$ denotes the matrix exponential.  Note that every rate 
matrix $\bar{P}$ satisfies $\bar{P}_{ij} \geq 0$ for $i \neq j$ and 
$\bar{P}_{ii} = -\sum_{j\in\mathcal{X}:j\neq i} \bar{P}_{ij} \leq 0$, so the row 
sums of $\bar{P}$ are 0.  Conditioned on $X_t = i$, the time until the next 
jump is exponentially-distributed with a mean of $-1/\bar{P}_{ii}$, and the 
probability to jump to $j$ is proportional to $\bar{P}_{ij}$ for $i \neq j$.

The natural notion of KL divergence rate $D(\bar{P}||\bar{Q})$ between 
transition rate matrices $\bar{P}$ and $\bar{Q}$ is 
\cite{Kesidis93,BaldiPiccioni98,deLaFortelle00, BertiniFaggionatoGabrielli14} 
the limiting log likelihood ratio, analogous to 
(\ref{eq:KL_div_MC_discrete_time_expression}):
\begin{align}
    D(\bar{P}\|\bar{Q}) &\triangleq \lim_{T \rightarrow \infty} \frac{1}{T} 
    \mathbb{E}_{p_{\bar{P},\mu_0}^{(T)}} 
    \log\left(\frac{p_{\bar{P},\mu_0}^{(T)}((X_t)_{t \in 
    [0,T]})}{p_{\bar{Q},\nu_0}^{(T)}((X_t)_{t \in [0,T]})}\right) \\
    &= \sum_{i \in \mathcal{X}} \pi(\bar{P})_i \sum_{j \in \mathcal{X}: j \neq i} 
    \left(\bar{Q}_{ij} - \bar{P}_{ij} + \bar{P}_{ij} 
    \log\left(\frac{\bar{P}_{ij}}{\bar{Q}_{ij}}\right)\right) 
    \label{eq:KL_div_MC_continuous_time_expression}
\end{align}
where $p_{\bar{P},\mu_0}^{(T)}(X_0^T)$ denotes the likelihood under rate matrix 
$\bar{P}$ and initial distribution $\mu_0$, and where $\pi(\bar{P})$ is the 
stationary distribution of rate matrix $\bar{P}$~\footnote{That is, 
$\pi(\bar{P}) \bar{P} = (0)_{i \in \mathcal{X}}$.  Equivalently, $\pi(\bar{P}) 
e^{\bar{P} t} = \pi(\bar{P})$ for all $t \geq 0$.}.  The quantity in the second 
summation in (\ref{eq:KL_div_MC_continuous_time_expression}) is the KL 
divergence between two Poisson distributions with means $\bar{P}_{ij}$ and 
$\bar{Q}_{ij}$.

The optimization problem analogous to (\ref{eq:optimization_problem_in_MC}) in 
continuous time is:
\begin{equation}
    \bar{P}^* = \argmin_{\bar{P}: \pi(\bar{P}) = \pi^*} D(\bar{P}||\bar{Q}) 
    \label{eq:optimization_problem_in_MC_continuous_time}
\end{equation}

\subsection{Comparison to prior work} \label{subsec: comparison}

Recent work \cite{HorowitzZhouEngland17,HorowitzEngland17} considers the 
question of the minimum power used to maintain a nonequilibrium state.  Their 
setting uses a different notion of cost than we do and also makes some 
restrictions about the base and controlled chains (they work in continuous 
time, assume that the base chain $\bar{Q}$ is reversible, and only allow 
controlled chains $\bar{P}$ such that $\bar{P}_{ij} \geq \bar{Q}_{ij}$ for all 
$i \neq j$ -- this corresponds to the biochemical mechanism of adding transitions with non-negative rates).  \cite{HorowitzZhouEngland17,HorowitzEngland17} minimize the 
entropy production rate among all controlled chains with the desired target 
distribution $\pi^*$ and find that ``fast control is optimal'': there is in 
general no optimally controlled chain, but given any chain $\bar{P}$ that has 
the target distribution $\pi^*$, we can come arbitrarily close to the minimum 
entropy production bound by speeding up $\bar{P}$ arbitrarily 
much\footnote{That is, using $c \bar{P} = (c \bar{P}_{ij})_{i,j\in\mathcal{X}}$ as 
the controlled chain and letting $c \rightarrow \infty$. This corresponds to the 
statement in \cite{HorowitzZhouEngland17} than the ``added edges (should) 
operate much faster than the equilibrium transitions''.} (while incurring an 
arbitrarily large KL divergence cost $D(\bar{P}||\bar{Q})$ according to our 
metric).

The difference between the two notions of cost -- KL divergence rate in this 
work and entropy production rate in 
\cite{HorowitzZhouEngland17,HorowitzEngland17} is the difference between total 
energy used and the efficiency with which that energy is used as measured by 
entropy production rate.  The very fast controlled chain of 
\cite{HorowitzZhouEngland17,HorowitzEngland17} uses a lot of energy 
efficiently, while our chain $\bar{P}^*$ minimizes energy use by the 
controller, but not the efficiency.  A consequence of this difference is that 
our optimal controlled chain $\bar{P}^*$ 
(\ref{eq:optimization_problem_in_MC_continuous_time}) depends on the 
uncontrolled chain $\bar{Q}$ (see section \ref{sec:minimum_power}), while the 
very fast close-to-optimal chain of 
\cite{HorowitzZhouEngland17,HorowitzEngland17} does not, except in the 
requirement that it be much faster than $\bar{Q}$.

\cite{Todorov06,Todorov09} introduced the KL divergence control setting.  
\cite{Todorov06,Todorov09} uses ``microscopic'' control cost, assigning a cost 
to a trajectory rather than a distribution over trajectories.  Similarly, the 
control goal in \cite{Todorov06,Todorov09} is microscopic (to reach a certain 
subset of the state space $\mathcal{X}$), rather than macroscopic (to maintain a 
target distribution $\pi^*$ over $\mathcal{X}$).

\cite{Gopalkrishnan16} considers the problem of erasing a bit of information 
encoded in the stationary distribution of a two-state continuous time Markov 
chain and uses the KL divergence to measure the cost of control, as does our 
work.  Whereas our control goal is to hold a target distribution $\pi^*$ and 
minimize the cost per unit time, \cite{Gopalkrishnan16}'s control goal is to 
have (in the notation of section \ref{sec:MC_continuous_time}) $\mu_T = \pi^*$ 
by a fixed time $T$ and to minimize the total cost used to achieve this.  
Consequently, \cite{Gopalkrishnan16} uses a time-varying controlled chain, 
while ours is constant in time.

\section{Minimum-power controlled chains} \label{sec:minimum_power}

In this section we minimize the power used to hold a desired nonequilibrium 
stationary distribution $\pi^*$. 

\begin{theorem} (Minimum-power chain) \label{thm:I_proj_stat_distrib} Let 
    $\pi^*$ be a distribution on finite set $\mathcal{X}$ with $\pi^*_i > 0$ for all 
    $i \in \mathcal{X}$.  Let $\pi(P)$ and $\pi(\bar{P})$ denote the stationary 
    distributions of discrete- and continuous-time chains $P$ and $\bar{P}$, 
    respectively.
    \begin{enumerate}
        \item (Discrete time)
   
            Let $Q$ be an irreducible, aperiodic transition probability matrix 
            (the uncontrolled discrete-time chain), let $\pi = \pi(Q)$, and let 
            $P^*$ denote the minimum-power controlled chain with the desired 
            stationary distribution $\pi^*$:
    \begin{equation}
        P^* = \argmin_{P: \pi(P) = \pi^*} D(P||Q) 
        \label{eq:def_I_proj_stat_distrib_disc_time}
    \end{equation}
    where the minimum is over all transition probability matrices with the 
            desired stationary distribution and where $D(P||Q)$ is as defined 
            in (\ref{eq:def_KL_div_MC_discrete_time_1}).  Then $P^*$ exists, is 
            unique, and satisfies for all $i,j\in\mathcal{X}$:
    \begin{equation}
        P^*_{ij} = Q_{ij} e^{\lambda_i + \eta_j} 
        \label{eq:I_proj_stat_distrib_dic_time_element}
    \end{equation}
    where $(\lambda_i)_{i \in \mathcal{X}}, (\eta_i)_{i \in \mathcal{X}}$ are real-valued 
            constants satisfying the recursive relations:
    \begin{align}
        \lambda_i &= - \log\bigg(\sum_{j \in \mathcal{X}} Q_{ij} e^{\eta_j}\bigg) 
        \label{eq:eta_disc_time} \\
        \eta_i &= \log(\pi^*_i) - \log\bigg(\sum_{j \in \mathcal{X}} \pi_j Q_{ji} 
        e^{\lambda_j}\bigg) \label{eq:lambda_disc_time}
    \end{align}
\item (Continuous time)
    Let $\bar{Q}$ be a transition rate matrix (the uncontrolled continuous-time 
            chain) with $e^{\bar{Q}}$ irreducible and aperiodic, and let 
            $\bar{P}^*$ denote the minimum-power controlled chain with the 
            desired stationary distribution $\pi^*$:
    \begin{equation}\label{eq:minpowerchain}
        \bar{P}^* = \argmin_{\bar{P}: \pi(\bar{P}) = \pi^*} D(\bar{P}||\bar{Q})
    \end{equation}
    where the minimum is over all transition rate matrices with the desired 
            stationary distribution and where $D(\bar{P}||\bar{Q})$ is as 
            defined in (\ref{eq:KL_div_MC_continuous_time_expression}).  Then 
            $\bar{P}^*$ exists, is unique, and
    satisfies:
    \begin{equation}
        \bar{P}^*_{ij} = \begin{cases}
            \bar{Q}_{ij} e^{\lambda_i - \lambda_j} &: i \neq j \\
            -\sum_{j \in \mathcal{X}: j \neq i} \bar{P}^*_{ij} &: i = j 
            \label{eq:def_I_proj_stat_distrib_cont_time}
        \end{cases}
    \end{equation}
    where $(\lambda_i)_{i \in \mathcal{X}}$ are real-valued constants satisfying the 
            recursive relations:
    \begin{align}
        \lambda_i = \frac{1}{2} \log\left(\frac{\sum_{j \in \mathcal{X}: j \neq i} 
        \pi^*_j \bar{Q}_{ji} e^{\lambda_j}}{\sum_{j \in \mathcal{X}: j \neq i} 
        \pi^*_i \bar{Q}_{ij} e^{-\lambda_j}}\right) \label{eq:lambda_cont_time}
    \end{align}
    \end{enumerate}
\end{theorem}

Existence and uniqueness of $P^*$ (\ref{eq:def_I_proj_stat_distrib_disc_time}) 
follow as a special case of Lemma 1 of \cite{CsiszarCoverChoi87}.  Existence 
and uniqueness of $\bar{P}^*$ (\ref{eq:def_I_proj_stat_distrib_cont_time}) were 
shown in \cite{BaldiPiccioni98}.  We prove expressions 
(\ref{eq:I_proj_stat_distrib_dic_time_element}) and 
(\ref{eq:def_I_proj_stat_distrib_cont_time}) for $P^*_{ij}$ and 
$\bar{P}^*_{ij}$
by setting up a Lagrange multiplier optimization problem, where $\alpha, 
(\lambda_i)_{i \in \mathcal{X}}, (\eta_i)_{i \in \mathcal{X}}$ are Lagrange multipliers.  
See Appendix \ref{app:proof} for proof details.


The recursive relations (\ref{eq:eta_disc_time}), (\ref{eq:lambda_disc_time}), 
and (\ref{eq:lambda_cont_time}) enable an iterative computation of the chains 
$P^*$ (\ref{eq:def_I_proj_stat_distrib_disc_time})
and $\bar{P}^*$ (\ref{eq:def_I_proj_stat_distrib_cont_time}).  In the 
continuous-time case, for example, we initialize $(\lambda_i^{(0)})_{i \in 
\mathcal{X}}$ to some value and then use relation (\ref{eq:lambda_cont_time}) to 
compute $\lambda^{(t)}_i$ as a function of $(\lambda^{(t-1)}_i)_{i \in \mathcal{X}}$ 
at the $t$-th iteration until numerical convergence.  The example in section 
\ref{sec:examples} is computed in this way.

The chains $P^*$ and $\bar{P}^*$ are the I-projections of the chains $P$ and 
$\bar{P}$ on the set of discrete- and continuous-time Markov chains, 
respectively, with a fixed stationary distribution.  The discrete-time case is 
presented in \cite{DemboZeitouniBook,CsiszarCoverChoi87} and the 
continuous-time case in \cite{BaldiPiccioni98}, where $D(P^*||Q)$ and 
$D(\bar{P}^*||\bar{Q})$ arise as large deviations rate functions for the 
empirical marginal distribution.
To our knowledge Theorem \ref{thm:I_proj_stat_distrib} presents the most 
explicit characterization of the I-projection in terms of the Lagrange 
multipliers.

We next specialize our results to the case of reversible uncontrolled Markov 
chains, a case important in equilibrium thermodynamics.  
Let $Q^\text{r} = (Q^\text{r}_{ij})_{i,j\in\mathcal{X}}$ denote the \textit{time-reverse} of a 
transition probability matrix $Q$.  That is:
\begin{equation}
    Q^\text{r}_{ij} = \frac{\pi_j}{\pi_i} Q_{ji} \label{eq:time_rev_chain_disc_time}
\end{equation}
where $\pi = \pi(Q) = \pi(Q^\text{r})$ is the stationary distribution.  Then 
$Q^\text{r}_{ij} = \lim_{t \rightarrow \infty} \mathbb{P}_Q(X_t=j|X_{t+1}=i)$; if $X_0 \sim 
\pi$, then $Q^\text{r}_{ij} = \mathbb{P}_Q(X_t=j|X_{t+1}=i)$ for all $t$.   Analogously, 
in continuous time, the time-reverse $\bar{Q}^\text{r}$ of a transition rate matrix 
$\bar{Q}$ satisfies $\bar{Q}_{ij} = \frac{\pi_j}{\pi_i} \bar{Q}_{ji}$ for all 
$i,j\in\mathcal{X}$.  A chain $Q$ is \textit{reversible} if $Q = Q^\text{r}$ 
(analogously, $\bar{Q} = \bar{Q}^\text{r}$ in continuous time).  

\begin{theorem} \label{thm:I_proj_stat_distrib_reversible_base_chain} 
    (Reversible uncontrolled chain) Let notation be as in the statement of 
    Theorem \ref{thm:I_proj_stat_distrib}.
   
   
   
   
   
    \begin{enumerate}
        \item If the uncontrolled transition probability matrix $Q$ is 
            reversible, then so is the minimum-power chain $P^*$ 
            (\ref{eq:def_I_proj_stat_distrib_disc_time}).  Analogously, in 
            continuous time, if the uncontrolled transition rate matrix 
            $\bar{Q}$ is reversible, then so is $\bar{P}^*$ 
            (\ref{eq:def_I_proj_stat_distrib_cont_time}).
        \item In discrete time, the Lagrange multipliers 
            (\ref{eq:eta_disc_time}), (\ref{eq:lambda_disc_time}) satisfy 
            $\eta_i = \lambda_i + \log(\pi^*_i/\pi_i)$ for all $i \in \mathcal{X}$.
        \item In continuous time, $\bar{P}^*$ satisfies for $i \neq j$:
            \begin{equation}
                \bar{P}^*_{ij} = \bar{Q}_{ij} \sqrt{\frac{\pi_i \pi^*_j}{\pi_j 
                \pi^*_i}} \label{eq:reversible_cont_time}
            \end{equation}
            and
            \begin{equation}
                D(\bar{P}^*||\bar{Q}) = \frac{1}{2}\sum_{i,j\in\mathcal{X}:i\neq j} 
                \left(\sqrt{\pi^*_i \bar{Q}_{ij}} - \sqrt{\pi^*_j 
                \bar{Q}_{ji}}\right)^2 \label{eq:KL_div_rev_cont_time}
            \end{equation}
    \end{enumerate}
\end{theorem}
\begin{proof}
    1) We can check that if $Q = Q^\text{r}$, then $D(P||Q) = D(P^\text{r}||Q)$ for all 
    $P$.  Suppose that $P^*$ is not reversible.  Let $P^{\leftrightarrow} 
    \triangleq (P^*+P^{*\text{r}})/2$.  Then $P^{\leftrightarrow}$ is reversible 
    and $\pi(P^{\leftrightarrow}) = \pi(P^*) = \pi^*$.  Since $D(P||Q)$ is 
    strictly convex in $P$, we have $D(P^{\leftrightarrow}||Q) < (D(P^*||Q) + 
    D(P^{*\text{r}}||Q))/2 = D(P^*||Q)$, contradicting the optimality of $P^*$.  
    Another proof: suppose that $P^*$ is not reversible, then $D(P^*||Q) = 
    D(P^{*\text{r}}||Q)$, contradicting the uniqueness of $P^*$ established in 
    Theorem \ref{thm:I_proj_stat_distrib}.  Therefore $P^*$ is reversible.  An 
    analogous argument proves $\bar{P}^*$ is reversible.  2) follows by using 
    time reversal (\ref{eq:time_rev_chain_disc_time}) twice along with the 
    reversibility of $Q$ and $P^*$, established in part 1):
    \begin{align}
        \pi^*_i Q_{ij}e^{\lambda_i + \eta_j} &= \pi^*_i P^*_{ij} = \pi^*_j 
        P^*_{ji} = \pi^*_j Q_{ji} e^{\lambda_j + \eta_i} \\
        &= \pi^*_j \frac{\pi_i}{\pi_j} Q_{ij} e^{\lambda_j + \eta_i}
    \end{align}
    and collecting $i$- and $j$-dependent terms to separate sides of the 
    equality to conclude that $\frac{\pi^*_i}{\pi_i} e^{\lambda_i - \eta_i} = 
    a$ for some constant $a$ for all $i \in \mathcal{X}$.  Choosing $a=1$ yields the 
    result.

    3) Expression (\ref{eq:reversible_cont_time}) for $\bar{P}^*$ is derived in 
    \cite{BaldiPiccioni98} and $D(\bar{P}^*||Q)$ 
    (\ref{eq:KL_div_rev_cont_time}) is derived as a large deviations rate 
    function in \cite{BertiniFaggionatoGabrielli14}.
\end{proof}

Part 2) of Theorem \ref{thm:I_proj_stat_distrib_reversible_base_chain} lets us 
simplify computation of $P^*$ (\ref{eq:def_I_proj_stat_distrib_disc_time}) 
somewhat when the uncontrolled chain $Q$ is reversible.

\section{Examples} \label{sec:examples}


We conclude with numerical examples of minimum-power controlled Markov chains 
with a target stationary distribution.  The first example is a two-state chain 
in discrete time and the second example is a birth-and-death chain in 
continuous and discrete time -- a toy model of the muscle fiber thought 
experiment in the introduction (section \ref{sec:introduction}).

\subsection{Two-state chain in discrete time}

Let $Q$ be a two-state discrete-time Markov chain on set $\mathcal{X} = \{1,2\}$ and 
let $\pi^* = (\pi^*_1, 1-\pi^*_1)$ be our desired nonequilibrium distribution 
with $\pi^*_1 \in (0,1)$.  All two-state chains are reversible, so we apply 
Theorem \ref{thm:I_proj_stat_distrib_reversible_base_chain} part 2) to compute 
the minimum-power controlled chain with stationary distribution $\pi^*$.  A 
computation shows the $2 \times 2$ minimum-power transition matrix $P^*$ 
(\ref{eq:def_I_proj_stat_distrib_disc_time}) has off-diagonal entries:
    \begin{equation}
        P^*_{ij} = \frac{1}{\pi^*_i} \cdot \frac{1-\sqrt{1-4 \pi^*_1 
        (1-\pi^*_1) s}}{2 s} \ : i \neq j
    \end{equation}
    where the second factor is independent of $i,j$ and
    \begin{equation}
        s \triangleq \frac{1}{Q_{1,2}} + \frac{1}{Q_{2,1}} - \frac{1}{Q_{1,2} 
        Q_{2,1}}
    \end{equation}
       
    The diagonal terms of $P^*$ are such that the row sums are $1$. 
   
   

\subsection{Birth-and-death chain}


We next present the example of the birth-and-death chain as a toy model of 
Feynman's muscle fiber thought experiment (see section \ref{sec:introduction}).  
For detailed models of molecular motors see \cite{Julicher} and \cite{Qian00}.

\subsubsection{Continuous time} \label{sec:cont_time_birth_and_death_chain}

\begin{figure}[!t]
    \centering
    \includegraphics[width=3in]{figures/birthanddeathchain.pdf}
    \caption{Uncontrolled chain $Q$ (a toy model for the slack arm pulled down 
    by gravity, with smaller-index states closer to the ground and rates $q_- > 
    q_+$): the birth-and-death chain with 5 states in continuous time.  The 
    states correspond to minima in the potential landscape experienced by a 
    single protein at different positions along the muscle fiber.}
    \label{fig:birth_and_death_chain}
\end{figure}

Let $\bar{Q}$ be the continuous-time birth-and-death chain on set 
$\{1,2,\ldots,K\}$ with parameters $q_+, q_- > 0$ depicted in Figure 
\ref{fig:birth_and_death_chain} for $K = 5$.  The chain transitions increments 
from state $i$ to $i + 1$ (resp. decrements to $i-1$) with rate $q_+$ (resp.  
$q_-$).  All other transitions, as well as decrementing from state $1$ and 
incrementing from state $K$, have rate $0$.  $\bar{Q}$ is reversible and its 
stationary distribution is, up to normalization, \cite{LevinPeresWilmerBook}:
\begin{equation}
    \pi_i \sim \left(\frac{q_+}{q_-}\right)^i 
    \label{eq:stat_distrib_birth_and_death_chain}
\end{equation}

Let our control objective be to maintain the target distribution $\pi^*$, a 
geometric distribution on $\{1,\ldots,K\}$:
\begin{equation}
    \pi^*_i \sim b^i \label{eq:target_distribution_birth_and_death_chain}
\end{equation}
where $b > 0$.  Then applying Theorems \ref{thm:I_proj_stat_distrib} and 
\ref{thm:I_proj_stat_distrib_reversible_base_chain} we find the minimum-power 
controlled chain $\bar{P}^*$ (\ref{eq:minpowerchain}) with stationary 
distribution $\pi^*$ to be another birth-and-death chain with increment and 
decrement rates $p_+^*, p_-^*$:
\begin{equation}
    p_{\pm}^* = (q_+ q_-)^{\frac{1}{2}} b^{\pm \frac{1}{2}} 
    \label{eq:birth_and_death_chain_cont_time_parameters}
\end{equation}
The cost per unit time of this birth-and-death chain $\bar{P}^*$ is 
(\ref{eq:KL_div_rev_cont_time}):
\begin{equation}
   
   
    D(\bar{P}^*||\bar{Q}^*) = \left(\sqrt{q_+} - \sqrt{q_- b}\right)^2 
    \left(\frac{1-b^{K-1}}{1-b^K}\right) 
    \label{eq:birth_and_death_chain_cont_time_cost}
\end{equation}
If $b < 1$, then as $K \rightarrow \infty$, we have $D(\bar{P}^*||\bar{Q}) \rightarrow 
\left(\sqrt{q_+} - \sqrt{q_- b}\right)^2$.

\begin{figure}[!t]
    \centering
   
    \includegraphics[width=3in]{figures/periodic_potential.pdf}
    \caption{A tilted periodic potential energy landscape.
   
    We consider the $i$-th local minimum of the potential as a discrete 
    position state with energy $U_i$ and energy difference $B = U_{i+1} - U_i$ 
    between adjacent states.  The activation energy to go left (resp. right) is 
    $A$ (resp. $A+B$).}
    \label{fig:energy_potential_snapshot}
\end{figure}

Recalling the motivating example of section \ref{sec:introduction}, we can 
think of the birth-and-death chain as a biased random walk, where the random 
walker tends to spend more time at small values of $i$ if $q_- > q_+$.  This is 
a toy model of a myosin protein moving on an actin filament in muscles -- a 
random walk in a tilted energy potential with periodically spaced minima 
corresponding to discrete steps along the fiber.  The energy potential is 
depicted in Figure \ref{fig:energy_potential_snapshot}: a transition from the 
$i+1$-th energy minimum to the $i$-th energy minimum must overcome activation 
energy $A$, and the reverse transition must overcome activation energy $B+A$, 
where $B$ is the energy difference between adjacent states.  In the bag-holding 
thought experiment, $B$ is the gravitational potential energy difference 
between adjacent states, with state $1$ being closest to the ground.

In terms of these energies, the increment and decrement rates $q_+, q_-$ are:
\begin{align}
    q_- &= c e^{-\beta A} \\
    q_+ &= c e^{-\beta (A + B)}
\end{align}
for some constant $c$, where $\beta$ is the inverse temperature.  The 
stationary distribution of the uncontrolled chain $\bar{Q}$ is, up to 
normalization:
\begin{equation}
    \pi_i \sim \left(\frac{q_+}{q_-}\right)^i = e^{-\beta B i}
\end{equation}

Let's write the target nonequilibrium distribution as:
\begin{equation}
    \pi^*_i \sim b^i = e^{-\beta B^* i} 
    \label{eq:target_distribution_cont_time_birth_and_death_chain}
\end{equation}
with energy difference $B^* = (-1/\beta) \log(b)$ between adjacent states.


Since the optimal controlled chain $\bar{P}^*$ is another birth-and-death 
chain, we can write its parameters $p^*_+, p^*_-$ 
(\ref{eq:birth_and_death_chain_cont_time_parameters}) in terms of energies 
$A^*$ and $B^*$:
\begin{align}
    p^*_- &= c e^{-\beta A^*} \\
    p^*_+ &= c e^{-\beta (A^* + B^*)}
\end{align}
where, using (\ref{eq:birth_and_death_chain_cont_time_parameters}), we find the 
controlled activation energy $A^*$:
\begin{equation}
    A^* = A + \frac{1}{2}\left(B - B^*\right)
\end{equation}
The cost per unit time of this birth-and-death chain $\bar{P}^*$ is 
(\ref{eq:birth_and_death_chain_cont_time_cost}):
\begin{equation}
   
   
    D(\bar{P}^*||\bar{Q}^*) = c e^{-\beta A} \left(e^{-\frac{\beta}{2} B} - 
    e^{-\frac{\beta}{2} B^*}\right)^2 \left(\frac{1-e^{-\beta B^* 
    (K-1)}}{1-e^{-\beta B^* K}}\right)
\end{equation}
where if $B^* > 0$, then the last factor tends to $1$ as $K \rightarrow \infty$.

In the muscle fiber thought experiment (where $B$ is the gravitational 
potential energy difference between adjacent states and state $0$ is closest to 
the ground) if $B^* < B$, then target distribution $\pi^*$ 
(\ref{eq:target_distribution_cont_time_birth_and_death_chain}) corresponds to 
imposing a constant force upwards (away from state $1$) on the random-walking 
myosin protein.  The control objective is macroscopic: rather than control the 
microscopic trajectory of a single myosin protein, we imagine controlling a 
large collection of identical, independent myosin proteins in the same way by 
imposing the same controlled chain $\bar{P}^*$ on all myosins; the bag-holder's 
arm position is determined by an average over the positions of this collection 
of myosins.

\subsubsection{Discrete time} \label{sec:discrete_time_birth_and_death_chain}

Let $Q$ be the discrete-time birth-and-death chain on set $\mathcal{X} = 
\{1,\ldots,K\}$ with transition probability $q_+$ (resp. $q_-$) to increment 
(resp. decrement) the state from $i$ to $i+1$ (resp. $i-1$).  The stationary 
distribution $\pi$ of $Q$ is as in 
(\ref{eq:stat_distrib_birth_and_death_chain}),
the same as in the continuous time case with transition rates $q_+, q_-$.  Let 
our control objective be to maintain the target distribution $\pi^*_i \sim 
b^i$, a geometric distribution on $\{1,\ldots,K\}$ 
(\ref{eq:target_distribution_birth_and_death_chain}).

Then in contrast to the continuous time case of section 
\ref{sec:cont_time_birth_and_death_chain}, the minimum-power controlled chain 
$P^*$ (\ref{eq:def_I_proj_stat_distrib_disc_time}) is not in general a 
birth-and-death chain.  Consider this numerical example: let the target 
nonequilibrium stationary distribution be as in 
(\ref{eq:target_distribution_birth_and_death_chain}) with $b = \sqrt{q_-/q_+}$:
\begin{equation}
    \pi^*_i \sim \sqrt{\pi_{K-i}} \sim \left(\frac{q_+}{q_-}\right)^{-i/2} 
    \label{eq:target_distribution_birth_and_death_chain_discrete_time}
\end{equation}
$\pi^*$ (\ref{eq:target_distribution_birth_and_death_chain_discrete_time}) is 
biased the other way from $\pi_i$, assigning most of its mass to large values 
of $i$ if $q_+ < q_-$.  The square root in 
(\ref{eq:target_distribution_birth_and_death_chain_discrete_time}) makes 
$\pi^*$ look more uniform than $\pi$.

Let $q_- = 0.2$, $q_+ = 0.1$, so that the stationary distribution of $Q$ is 
$\pi_i \sim 2^{-i}$, and $\pi^*_i \sim 2^{i/2}$.  Figure 
\ref{fig:birth_and_death_chain_I_proj} shows the non-zero off-diagonal elements 
of $P^*_{ij}$ with three-digit precision; the increment and decrement 
probabilities vary with state $i$, so $P^*$ is not a birth-and-death chain.

\begin{figure}[!t]
    \centering
    \includegraphics[width=3in]{figures/birthanddeathchain_I_proj.pdf}
    \caption{The minimum-power discrete-time controlled chain $P^*$ 
    (\ref{eq:def_I_proj_stat_distrib_disc_time}) for the uncontrolled 
    birth-and-death chain $Q$ with parameters $q_- = 0.2$, $q_+ = 0.1$ and 
    target stationary distribution $\pi^*_i \sim 2^{i/2}$.  Transitions from a 
    state to itself are not shown.}
    \label{fig:birth_and_death_chain_I_proj}
\end{figure}

Finally, Figure \ref{fig:birth_and_death_chain_cost} depicts the time evolution 
of distribution $\mu_t = \mu_0 P^{*t}$ with $\mu_0 = \pi$, and the cost 
$D_{\mu_t}(P^*||Q)$ converging to the minimum power $D(P^*||Q) \approx 0.0315$ 
to maintain the nonequilibrium distribution $\pi^*$.

\begin{figure}[!t]
    \centering
    \includegraphics[width=3.8in]{figures/birthanddeathchaincost.pdf}
    \caption{(top) Time evolution of the $|\mathcal{X}| = 5$ components of $\mu_t = 
    \mu_0 (P^{*})^t$ with $\mu_0 = \pi$, where $P^*$ is the minimum cost chain 
    (\ref{eq:def_I_proj_stat_distrib_disc_time}), showing $\mu_t \rightarrow \pi^*$ as 
    $t \rightarrow \infty$.  (bottom) The cost (blue) $D_{\mu_t}(P^*||Q)$ 
    (\ref{eq:cost_disc_time_KL_div}) and the minimum power (red) $D(P^*||Q)$.}
    \label{fig:birth_and_death_chain_cost}
\end{figure}

Returning to the molecular motor picture, the discrete-time chain $P^*$ differs 
from its continuous-time cousin of section 
\ref{sec:cont_time_birth_and_death_chain} in that $P^*$ is not a 
birth-and-death chain.  Thus an optimal discrete-time control policy modifies 
the base birth-and-death chain $Q$ in a way that depends on state $i$, and so 
can't be thought of as corresponding to a constant, state-independent force 
upwards (away from state $1$) as in the continuous time case.


\section*{Acknowledgment}

The authors gratefully acknowledge Hideo Mabuchi for suggestions and insightful discussions.

\bibliographystyle{IEEEtran}
