The I-Projection, Pythagorean Identity, and Iterative Proportional Fitting (Blahut-Arimoto)

Intuition Beginner

You have two probability distributions. Call one the "target" $q$ and the other a "candidate" $p$. The KL divergence $D(p\Vert q)$ measures how different $p$ is from $q$. Now suppose you want to find the distribution $p$ that is closest to $q$ (in the KL sense) but belongs to a constrained set. For example, $p$ might need to match certain marginal distributions.

This is the I-projection: project $q$ onto a set of distributions by finding the member of the set that minimises $D(p\Vert q)$. Think of it as finding the closest point to $q$ inside a fenced region, where "closest" means least KL divergence.

A beautiful property of I-projection is the Pythagorean identity. If $p^{\ast }$ is the I-projection of $q$ onto a convex set $\mathcal{P}$, then for any $p\in \mathcal{P}$:

$$D(p\Vert q)=D(p\Vert p^{\ast })+D(p^{\ast }\Vert q).$$

This looks like the Pythagorean theorem: the squared distance from $p$ to $q$ decomposes into the distance from $p$ to the projection $p^{\ast }$ plus the distance from $p^{\ast }$ to $q$. The I-projection splits the divergence into orthogonal components.

The Blahut-Arimoto algorithm uses this geometry to compute channel capacity. It alternates between two steps: (1) find the best output distribution given the current input distribution, and (2) find the best input distribution given the current output distribution. Each step is an I-projection, and the Pythagorean identity guarantees that the mutual information increases monotonically. The algorithm converges to the capacity-achieving input distribution.

Visual Beginner

Concept	Analogy	Mathematical object
KL divergence $D(p|q)$	Squared Euclidean distance	Measure of dissimilarity
I-projection	Orthogonal projection onto a set	$\arg {\min }_{p\in \mathcal{P}}D(p|q)$
Pythagorean identity	$a^2+b^2=c^2$	$D(p|q)=D(p|p^{\ast })+D(p^{\ast }|q)$
Blahut-Arimoto	Alternating projection (von Neumann)	Iterative I-projections converge to capacity

Figure: a Venn-diagram-style picture. A point $q$ sits outside a convex set $\mathcal{P}$. The I-projection $p^{\ast }$ is the closest point in $\mathcal{P}$ to $q$ (under KL divergence). For any other point $p\in \mathcal{P}$, the three points form a right triangle: $D(p\Vert q)$ is the hypotenuse, $D(p\Vert p^{\ast })$ is one leg (inside $\mathcal{P}$), and $D(p^{\ast }\Vert q)$ is the other leg (the projection direction).

Worked example Beginner

A binary symmetric channel flips each bit with probability 0.1. What input distribution maximises the mutual information?

By symmetry, the uniform input $p(0)=p(1)=0.5$ should be optimal. The Blahut-Arimoto algorithm confirms this. Starting from any input distribution, it alternates:

Step A. Given input $p(x)$, compute the output $q(y)$ by mixing over all inputs with their channel transition probabilities.

Step B. Update the input: $p_{\text{new}}(x)\propto p(x)\cdot \exp (I(x;Y))$ where $I(x;Y)$ measures the information that sending symbol $x$ provides about the output.

Starting from $p(0)=0.3$, $p(1)=0.7$:

The output distribution is $q(0)=0.3\times 0.9+0.7\times 0.1=0.34$, $q(1)=0.66$. After one iteration, the algorithm pushes $p$ toward uniform. After a few iterations, $p(0)\approx p(1)\approx 0.5$, and the mutual information converges to $C=1-H_b(0.1)=1-0.469=0.531$ bits/use.

Check your understanding Beginner

Formal definition Intermediate+

Let $\mathcal{P}$ be a convex set of probability distributions on a finite alphabet $\mathcal{X}$, and let $q$ be a distribution on $\mathcal{X}$ with full support.

Definition (I-projection). The I-projection of $q$ onto $\mathcal{P}$ is:

$$p^{\ast }=\arg \underset{p\in \mathcal{P}}{\min }D(p\Vert q)=\arg \underset{p\in \mathcal{P}}{\min }\sum _{x}p(x)\log \frac{p(x)}{q(x)}.$$

Since $D(p\Vert q)$ is strictly convex in $p$ (for fixed $q$ with full support) and $\mathcal{P}$ is convex, $p^{\ast }$ exists and is unique.

Theorem (Pythagorean identity for I-projection). If $\mathcal{P}$ is a linear family (a set of distributions satisfying linear constraints of the form ${\sum }_{x}p(x)a_i(x)=b_i$), then for any $p\in \mathcal{P}$:

$$D(p\Vert q)=D(p\Vert p^{\ast })+D(p^{\ast }\Vert q).$$

Proof. Since $\mathcal{P}$ is defined by linear constraints, the optimal $p^{\ast }$ satisfies the KKT conditions. The Lagrangian is:

$$\mathcal{L}(p,\lambda )=D(p\Vert q)+\sum _i{\lambda }_i\left(b_i-\sum _{x}p(x)a_i(x)\right).$$

Setting the derivative to zero gives $p^{\ast }(x)=q(x)\exp \left(-1-{\sum }_i{\lambda }_i{a}_i(x)\right)/Z$, i.e., $p^{\ast }$ has the exponential family form $p^{\ast }(x)=q(x)e^{\theta \cdot a(x)}/Z(\theta )$.

For any $p\in \mathcal{P}$:

$$D(p\Vert q)-D(p^{\ast }\Vert q)=\sum _x(p(x)-p^{\ast }(x))\log \frac{p^{\ast }(x)}{q(x)}=\sum _x(p(x)-p^{\ast }(x))\left(\theta \cdot a(x)-\log Z\right).$$

Since both $p$ and $p^{\ast }$ satisfy the same linear constraints, ${\sum }_x(p(x)-p^{\ast }(x))a_i(x)=0$ for all $i$. Therefore:

$$D(p\Vert q)-D(p^{\ast }\Vert q)=-\log Z\sum _x(p(x)-p^{\ast }(x))=0.$$

This gives $D(p\Vert q)=D(p^{\ast }\Vert q)+{\sum }_x(p(x)-p^{\ast }(x))\log (p^{\ast }(x)/q(x))$. The remaining term equals $D(p\Vert p^{\ast })$ by direct computation: $D(p\Vert p^{\ast })={\sum }_{x}p(x)\log (p(x)/p^{\ast }(x))$, and using $p^{\ast }(x)=q(x)e^{\theta \cdot a(x)}/Z$ with the linear constraints gives the result. $\square$

The Blahut-Arimoto algorithm

Definition (Blahut-Arimoto for channel capacity). For a discrete memoryless channel $p(y|x)$, the Blahut-Arimoto algorithm computes:

$$C=\underset{p(x)}{\max }I(X;Y)=\underset{p(x)}{\max }\sum _{x,y}p(x)p(y|x)\log \frac{p(y|x)}{q(y)}$$

by iterating:

1. Output update: $q_k(y)={\sum }_x{p}_k(x)p(y|x)$
2. Input update: $p_{k+1}(x)=\frac{p_k(x)\exp (I_k(x))}{{\sum }_{x^{\prime }}{p}_k(x^{\prime })\exp (I_k(x^{\prime }))}$

where $I_k(x)=D(p(y|x)\Vert q_k(y))={\sum }_{y}p(y|x)\log \frac{p(y|x)}{q_k(y)}$.

The capacity estimate at step $k$ is $C_k={\sum }_x{p}_k(x)I_k(x)$, and $C_k\nearrow C$ monotonically.

Connection to iterative proximal methods

The Blahut-Arimoto algorithm can be viewed as an alternating proximal step in the geometry of KL divergence 44.06.01. Each update minimises:

$$p_{k+1}=\arg \underset{p}{\min }\left[-I(X;Y)(p)+D(p\Vert p_k)\right],$$

which is a proximal iteration with KL divergence as the Bregman distance. This connects to mirror descent and other first-order methods in convex optimisation.

Key theorem with proof Intermediate+

Theorem (Blahut-Arimoto convergence). The Blahut-Arimoto algorithm converges to the channel capacity $C$ from any initial distribution $p_0$ with full support. The iterates satisfy:

$$C_{k+1}\ge C_k\text{and}\underset{k\to \infty }{\lim }{C}_k=C.$$

Proof. Define the double maximisation:

$$C=\underset{p(x)}{\max }\underset{q(y)}{\max }\sum _{x,y}p(x)p(y|x)\log \frac{p(y|x)}{q(y)}.$$

The inner maximisation over $q$ (for fixed $p$) is achieved by $q^{\ast }(y)={\sum }_{x}p(x)p(y|x)$. The outer maximisation over $p$ (for fixed $q$) is achieved by $p^{\ast }(x)\propto p(x)\exp (D(p(y|x)\Vert q(y)))$.

Each step of Blahut-Arimoto performs one of these maximisations, so the objective increases monotonically. By the alternating maximisation lemma (a consequence of convexity), the iterates converge to the global maximum. $\square$

Bridge. The I-projection and Pythagorean identity build toward the iterative algorithms in 44.06.01 where proximal methods and mirror descent use the same Bregman-divergence structure; this appears again in 46.02.05 where the rate-distortion function is computed by the dual Blahut-Arimoto algorithm; the foundational reason the Pythagorean identity holds is the linear constraint structure of the projection set, which is dual to the exponential family form of the optimal distribution; the central insight is that alternating I-projections converge monotonically, and putting these together gives a unified computational framework for all information-theoretic optimisation problems.

Exercises Intermediate+

Lean formalization Intermediate+

Mathlib has no formalisation of I-projection, the Pythagorean identity for KL divergence, or the Blahut-Arimoto algorithm. The KL divergence $D(p\Vert q)$ is not defined as a functional on distributions in Mathlib. The optimisation over conditional distributions required for the algorithm is beyond Mathlib's current capabilities. The convergence proof requires the alternating maximisation lemma, which is a result in convex optimisation not present in Mathlib. The exponential family form of the I-projection requires solving the KKT conditions for constrained optimisation over probability distributions, which Mathlib's optimisation library does not support. This unit ships without a lean_module.

Advanced results Master

The information geometry of I-projection

Imre Csiszar's 1975 paper established that the space of probability distributions equipped with KL divergence has a differential-geometric structure. The I-projection is the analogue of orthogonal projection in this geometry. The tangent space at a distribution $p$ consists of score functions $f(x)$ with ${\mathbb{E}}_p[f]=0$, and the KL divergence induces a Riemannian metric given by the Fisher information matrix $g_{ij}={\mathbb{E}}_p[({\partial }_i\log p)({\partial }_j\log p)]$.

The Pythagorean identity is a manifestation of the flatness of the exponential family manifold. Linear families (defined by moment constraints) and exponential families (defined by exponential tilting) are dual flat submanifolds in this geometry. The I-projection onto a linear family is the unique point where the exponential family "orthogonal" to the linear family intersects it.

Iterative proportional fitting (IPF)

The iterative proportional fitting procedure, also known as the RAS algorithm or matrix scaling, is a special case of alternating I-projections. Given a target matrix of row and column marginals, IPF alternates between scaling rows to match row marginals and scaling columns to match column marginals. Each scaling step is an I-projection onto the set of matrices with the given marginal.

IPF converges by the same alternating projection principle as Blahut-Arimoto. The convergence rate depends on the condition number of the marginal constraints. The algorithm is used in contingency table analysis, transportation planning, and the calibration of input-output matrices in economics.

The alternating minimisation lemma

The general principle underlying both Blahut-Arimoto and IPF is the alternating minimisation lemma: if $f(p,q)$ is convex in each argument separately and the feasible sets are convex, then the sequence generated by alternating minimisation converges to a stationary point. For the specific structure of KL divergence, the stationary point is the global optimum.

The connection to proximal methods 44.06.01 is that each step of Blahut-Arimoto can be written as:

$$p_{k+1}=\arg \underset{p}{\min }\left\{f(p)+\frac{1}{\eta }D(p\Vert p_k)\right\},$$

which is a mirror descent update with step size $\eta$ and KL divergence as the Bregman distance. The convergence rate is $O(1/k)$ for the capacity estimate, matching the rate of other first-order methods.

Synthesis. The I-projection and Pythagorean identity form the geometric backbone of information theory; the central insight is that KL divergence induces a flat geometry on the probability simplex, and putting these together we see that Blahut-Arimoto, IPF, and mirror descent are all instances of alternating I-projections in this geometry; this builds toward the computation of channel capacity 46.03.01 and rate-distortion 46.02.05 where these algorithms are deployed; the foundational reason for convergence is the convexity of KL divergence in its first argument combined with the linear structure of the constraint sets, and this generalises to any Bregman-divergence-based alternating optimisation.

Full proof set Master

Proposition (Existence and uniqueness of I-projection). For a convex set $\mathcal{P}$ of distributions and a distribution $q$ with full support, the I-projection $p^{\ast }=\arg {\min }_{p\in \mathcal{P}}D(p\Vert q)$ exists and is unique.

Proof. The function $p\mapsto D(p\Vert q)$ is strictly convex on the simplex (since its Hessian is the Fisher information, which is positive definite). The feasible set $\mathcal{P}$ is convex by assumption. A strictly convex function on a convex set has at most one minimiser.

For existence, consider a minimising sequence $p_n\in \mathcal{P}$ with $D(p_n\Vert q)\to {\inf }_{p\in \mathcal{P}}D(p\Vert q)$. Since $\mathcal{P}$ is a subset of the compact simplex, a subsequence converges to some $p^{\ast }$. By continuity of $D(p\Vert q)$ in $p$ (for fixed $q$ with full support), $D(p^{\ast }\Vert q)=\inf D(p\Vert q)$. By strict convexity, $p^{\ast }$ is unique. $\square$

Proposition (Dual Blahut-Arimoto for rate-distortion). The rate-distortion function $R(D)$ can be computed by an alternating minimisation: initialise $q_0(\hat{x})$, then iterate $p_{k+1}(\hat{x}|x)\propto q_k(\hat{x})e^{s\cdot d(x,\hat{x})}$ and $q_{k+1}(\hat{x})={\sum }_{x}p(x)p_{k+1}(\hat{x}|x)$. The sequence $R_k=I(X;\hat{X}{)}_k$ converges to $R(D)$ where $s$ is the Lagrange multiplier corresponding to $D$.

Proof. The Lagrangian formulation gives $R(D)={\max }_{s\le 0}{\min }_{p(\hat{x}|x)}[I(X;\hat{X})-s(\mathbb{E}[d]-D)]$. For fixed $s$, the inner minimisation is over $p(\hat{x}|x)$, which decomposes into independent problems for each $x$:

$$\underset{p(\hat{x}|x)}{\min }\left[\sum _{\hat{x}}p(\hat{x}|x)\log \frac{p(\hat{x}|x)}{q(\hat{x})}-s\sum _{\hat{x}}p(\hat{x}|x)d(x,\hat{x})\right].$$

Setting the derivative to zero gives $p^{\ast }(\hat{x}|x)\propto q(\hat{x})e^{s\cdot d(x,\hat{x})}$. Given $p(\hat{x}|x)$, the optimal $q$ is $q^{\ast }(\hat{x})={\sum }_{x}p(x)p(\hat{x}|x)$ (the marginal). Alternating these two steps converges by the alternating minimisation lemma. $\square$

Connections Master

• 46.04.01 — Hypothesis testing uses KL divergence as the error exponent; the I-projection of the null onto the alternative's typical set determines the optimal test.
• 46.01.02 — KL divergence and mutual information are the basic quantities that the I-projection minimises.
• 46.03.01 — Channel capacity is computed by Blahut-Arimoto; the capacity-achieving distribution is the fixed point of the I-projection iteration.
• 44.06.01 — Proximal and subgradient methods are the optimisation-theoretic counterparts of Blahut-Arimoto; mirror descent with KL divergence is equivalent.
• 46.02.05 — The rate-distortion function is computed by the dual Blahut-Arimoto algorithm; the same I-projection geometry governs both.

Historical & philosophical context Master

Imre Csiszar introduced the I-projection in his 1975 paper "I-Divergence Geometry of Probability Distributions and Identification Problems" (Periodica Mathematica Hungarica 6(2-3):149-161). Csiszar recognised that the KL divergence endows the probability simplex with a differential-geometric structure, and the I-projection plays the role of orthogonal projection in this geometry. The Pythagorean identity is the information-geometric analogue of the Euclidean Pythagorean theorem.

Richard Blahut and Suguru Arimoto independently discovered the alternating maximisation algorithm for channel capacity in 1972. Blahut's paper ("Computation of Channel Capacity and Rate-Distortion Functions," IEEE Trans. IT 18(4):460-473) presented the algorithm in full generality for both channel capacity and rate-distortion. Arimoto's paper ("An Algorithm for Computing the Capacity of Arbitrary Discrete Memoryless Channels," IEEE Trans. IT 18(1):14-20) focused on channel capacity and proved the monotonic convergence property.

The iterative proportional fitting procedure predates both: it was introduced by Deming and Stephan in 1940 ("On a Least Squares Adjustment of a Sampled Frequency Table When the Expected Marginal Totals Are Known," Annals of Mathematical Statistics 11(4):427-444) for contingency table adjustment. The connection to I-projection and KL divergence was made by Csiszar in 1975, who showed that IPF is alternating I-projection onto marginal constraint sets.

The information geometry perspective was developed by Amari and Nagaoka in their 2000 book Methods of Information Geometry (AMS/Oxford), which established the dual flat structure of the exponential family manifold and unified the geometric treatment of I-projection, Pythagorean identities, and alternating projections.

Bibliography Master

@article{csiszar1975,
  author  = {Csisz{\'a}r, I.},
  title   = {I-Divergence Geometry of Probability Distributions and Identification Problems},
  journal = {Periodica Mathematica Hungarica},
  volume  = {6},
  number  = {2-3},
  pages   = {149--161},
  year    = {1975},
}
@article{blahut1972,
  author  = {Blahut, R. E.},
  title   = {Computation of Channel Capacity and Rate-Distortion Functions},
  journal = {IEEE Transactions on Information Theory},
  volume  = {18},
  number  = {4},
  pages   = {460--473},
  year    = {1972},
}
@article{arimoto1972,
  author  = {Arimoto, S.},
  title   = {An Algorithm for Computing the Capacity of Arbitrary Discrete Memoryless Channels},
  journal = {IEEE Transactions on Information Theory},
  volume  = {18},
  number  = {1},
  pages   = {14--20},
  year    = {1972},
}
@article{deming-stephan1940,
  author  = {Deming, W. E. and Stephan, F. F.},
  title   = {On a Least Squares Adjustment of a Sampled Frequency Table When the Expected Marginal Totals Are Known},
  journal = {Annals of Mathematical Statistics},
  volume  = {11},
  number  = {4},
  pages   = {427--444},
  year    = {1940},
}
@book{amari-nagaoka2000,
  author    = {Amari, S. and Nagaoka, H.},
  title     = {Methods of Information Geometry},
  publisher = {AMS/Oxford University Press},
  year      = {2000},
}
@book{cover-thomas2006,
  author    = {Cover, T. M. and Thomas, J. A.},
  title     = {Elements of Information Theory},
  edition   = {2nd},
  publisher = {Wiley},
  year      = {2006},
}