Blahut–Arimoto algorithm

For the case of channel capacity, the algorithm was independently invented by Suguru Arimoto^[1] and Richard Blahut.^[2] In addition, Blahut's treatment gives algorithms for computing rate distortion and generalized capacity with input contraints (i.e. the capacity-cost function, analogous to rate-distortion). These algorithms are most applicable to the case of arbitrary finite alphabet sources. Much work has been done to extend it to more general problem instances.^[3][4] Recently, a version of the algorithm that accounts for continuous and multivariate outputs was proposed with applications in cellular signaling.^[5] There exists also a version of Blahut–Arimoto algorithm for directed information.^[6]

A discrete memoryless channel (DMC) can be specified using two random variables ${\displaystyle X,Y}$ with alphabet ${\displaystyle {\mathcal {X}},{\mathcal {Y}}}$, and a channel law as a conditional probability distribution ${\displaystyle p(y|x)}$. The channel capacity, defined as ${\displaystyle C:=\sup _{p\sim X}I(X;Y)}$, indicates the maximum efficiency that a channel can communicate, in the unit of bit per use.^[7] Now if we denote the cardinality ${\displaystyle |{\mathcal {X}}|=n,|{\mathcal {Y}}|=m}$, then ${\displaystyle p_{Y|X}}$ is a ${\displaystyle n\times m}$ matrix, which we denote the ${\displaystyle i^{th}}$ row, ${\displaystyle j^{th}}$ column entry by ${\displaystyle w_{ij}}$. For the case of channel capacity, the algorithm was independently invented by Suguru Arimoto^[8] and Richard Blahut.^[9] They both found the following expression for the capacity of a DMC with channel law:

${\displaystyle C=\max _{\mathbf {p} }\max _{Q}\sum _{i=1}^{n}\sum _{j=1}^{m}p_{i}w_{ij}\log \left({\dfrac {Q_{ji}}{p_{i}}}\right)}$

where ${\displaystyle \mathbf {p} }$ and ${\displaystyle Q}$ are maximized over the following requirements:

• ${\displaystyle \mathbf {p} }$ is a probability distribution on ${\displaystyle X}$, That is, if we write ${\displaystyle \mathbf {p} }$ as ${\displaystyle (p_{1},p_{2}...,p_{n}),\sum _{i=1}^{n}p_{i}=1}$
• ${\displaystyle Q=(q_{ji})}$ is a ${\displaystyle m\times n}$ matrix that behaves like a transition matrix from ${\displaystyle Y}$ to ${\displaystyle X}$ with respect to the channel law. That is, For all ${\displaystyle 1\leq i\leq n,1\leq j\leq m}$:
  • ${\displaystyle q_{ji}\geq 0,q_{ji}=0\Leftrightarrow w_{ij}=0}$
  • Every row sums up to 1, i.e. ${\displaystyle \sum _{i=1}^{n}q_{ji}=1}$.

Then upon picking a random probability distribution ${\displaystyle \mathbf {p} ^{0}:=(p_{1}^{0},p_{2}^{0},...p_{n}^{0})}$ on ${\displaystyle X}$, we can generate a sequence ${\displaystyle (\mathbf {p} ^{0},Q^{0},\mathbf {p} ^{1},Q^{1}...)}$ iteratively as follows:

${\displaystyle (q_{ji}^{t}):={\dfrac {p_{i}^{t}w_{ij}}{\sum _{k=1}^{n}p_{k}^{t}w_{kj}}}}$

${\displaystyle p_{k}^{t+1}:={\dfrac {\prod _{j=1}^{m}(q_{jk}^{t})^{w_{kj}}}{\sum _{i=1}^{n}\prod _{j=1}^{m}(q_{ji}^{t})^{w_{ij}}}}}$

For ${\displaystyle t=0,1,2...}$.

Then, using the theory of optimization, specifically coordinate descent, Yeung^[10] showed that the sequence indeed converges to the required maximum. That is,

${\displaystyle \lim _{t\to \infty }\sum _{i=1}^{n}\sum _{j=1}^{m}p_{i}^{t}w_{ij}\log \left({\dfrac {Q_{ji}^{t}}{p_{i}^{t}}}\right)=C}$.

So given a channel law ${\displaystyle p(y|x)}$, the capacity can be numerically estimated up to arbitrary precision.

Suppose we have a source ${\displaystyle X}$ with probability ${\displaystyle p(x)}$ of any given symbol. We wish to find an encoding ${\displaystyle p({\hat {x}}|x)}$ that generates a compressed signal ${\displaystyle {\hat {X}}}$ from the original signal while minimizing the expected distortion ${\displaystyle \langle d(x,{\hat {x}})\rangle }$, where the expectation is taken over the joint probability of ${\displaystyle X}$ and ${\displaystyle {\hat {X}}}$. We can find an encoding that minimizes the rate-distortion functional locally by repeating the following iteration until convergence:

${\displaystyle p_{t+1}({\hat {x}}|x)={\frac {p_{t}({\hat {x}})\exp(-\beta d(x,{\hat {x}}))}{\sum _{\hat {x}}p_{t}({\hat {x}})\exp(-\beta d(x,{\hat {x}}))}}}$

${\displaystyle p_{t+1}({\hat {x}})=\sum _{x}p(x)p_{t+1}({\hat {x}}|x)}$

where ${\displaystyle \beta }$ is a parameter related to the slope in the rate-distortion curve that we are targeting and thus is related to how much we favor compression versus distortion (higher ${\displaystyle \beta }$ means less compression).