Finite-alphabet Blahut–Arimoto rate-distortion algorithm

For the exact-real recurrence, write $[m]=\{1,\dots ,m\}$ and $[n]=\{1,\dots ,n\}$ for positive integers $m,n$. Let $p_1,\dots ,p_m>0$ be source probabilities summing to one, let $d_{ij}\ge 0$ be finite distortions for $(i,j)\in [m]\times [n]$, let $\beta \ge 0$ be finite, and let $q_0$ be any full-support reproduction pmf. For $t=0,1,\dots$ set

$$Z_i(q_t)=\sum _{j=1}^n{q}_t(j)e^{-\beta d_{ij}},W_t(j\mid i)=\frac{q_t(j)e^{-\beta d_{ij}}}{Z_i(q_t)},q_{t+1}(j)=\sum _{i=1}^m{p}_i{W}_t(j\mid i)\mathrm{.}$$

All denominators are positive, every row of $W_t$ sums to one, and $q_{t+1}$ has full support. Define the nats objective $L_{\beta }(W)=I_{\mathrm{n}\mathrm{a}\mathrm{t}}(X;\widehat{X})+\beta \mathbb{E}d(X,\widehat{X})$, the current channel value $L_t=L_{\beta }(W_t)$, and the a-posteriori factor

$${\gamma }_t:=\underset{1\le j\le n}{\max }\sum _{i=1}^m{p}_i\frac{e^{-\beta d_{ij}}}{Z_i(q_t)}\ge 1.$$

For $\tau >0$, the output is the first processed iteration index $t$ for which a rigorous upper interval endpoint $u_t\ge \ln {\gamma }_t$ satisfies $u_t\le \tau$, together with $W_t$, $q_{t+1}$, $L_t$, and interval enclosures of these values and $\ln {\gamma }_t$. Compute the objective by $L_t={\sum }_{i=1}^m{\sum }_{j=1}^n{p}_i{W}_t(j\mid i)\ln \frac{W_t(j\mid i)}{q_{t+1}(j)}+\beta {\sum }_{i=1}^m{\sum }_{j=1}^n{p}_i{W}_t(j\mid i)d_{ij}\mathrm{.}$ The exact-real recurrence defines all values being enclosed.

The certified interval procedure is defined only when every $p_i,d_{ij},\beta ,q_0(j)$ and $\tau$ is rational, with the positivity, normalization, finiteness and full-support conditions above and $\tau >0$. Execute the following loops. For $t=0,1,2,\dots$, try precisions $r=2,3,4,\dots$ in order. At precision $r$, use the dyadic format $(p,e_{\min },e_{\max })=(r^2,-r,r)$ and truncation index $N=r$ in Certified range evaluation for elementary real and complex functions. Round each rational input outward to an interval in this format. Build a finite typed expression DAG for the recurrence from $q_0$ through $q_t$, and request enclosures for every $W_t(j\mid i)$, $q_{t+1}(j)$, $L_t$, and $s_{j,t}:={\sum }_i{p}_i{e}^{-\beta d_{ij}}\mathrm{/}{Z}_i(q_t)$. If rounding any input overflows or the evaluator returns a non-finite status, increase $r$ and repeat this precision loop. Otherwise form the interval $G_t=[{\max }_j{\underset{‾}{s}}_{j,t},{\max }_j{\bar{s}}_{j,t}]$ from the rational endpoints of the $s_{j,t}$ enclosures, and evaluate $\ln G_t$ with the same format and truncation index. If this evaluation is non-finite, increase $r$ and repeat. Continue increasing $r$ until all requested output enclosures, including that of $\ln {\gamma }_t$, have width less than $1\mathrm{/}(t+1)$. Let $u_t$ be the upper endpoint of the resulting $\ln G_t$ enclosure. If $u_t\le \tau$, return the stated output and enclosures; otherwise end the precision loop and proceed to $t+1$. For arbitrary exact-real inputs the recurrence remains mathematically defined, but this interval implementation is specified only under the rational representability conditions. No bit-complexity bound is asserted.