Iterative Proportional Fitting (IPF): Algorithm, Applications, and Implications

• Iterative Proportional Fitting (IPF) is an algorithm designed to adjust a nonnegative table or probability distribution to meet specified marginal, aggregate, or moment constraints, using repetitive multiplicative updates to converge to a constrained distribution.
• IPF has broad applications in matrix scaling, Bayesian networks, graphical models, and optimization problems, enabling adjustments to match pre-defined distributions and maintaining the stability where necessary.
• When applied to matrix scaling, IPF alternately rescales rows and columns while preserving the zeros in the initial matrix and replacing sparse feasibility checks using Hall-types conditions.

Iterative proportional fitting (IPF) is an algorithm for modifying a nonnegative table or probability distribution so that specified marginal, aggregate, or moment constraints are satisfied. Starting from a seed distribution, IPF repeatedly applies multiplicative adjustments to one constraint at a time. Each update is a Kullback–Leibler (KL) or I-projection onto the corresponding constraint set, and under appropriate feasibility and support conditions the sequence converges to a constrained distribution closest to the seed in KL divergence. In matrix form, IPF alternately rescales rows and columns; in probabilistic models, it fits prescribed marginals; in entropy-regularized optimal transport, it is also known as Sinkhorn’s algorithm or matrix scaling.

1. Mathematical formulation and information geometry

Let $Q_0(x)$ be a strictly positive initial distribution over a finite configuration space, and let $R_i(y_i)$ specify a target distribution for a subset $Y_i$ of variables. The feasible set associated with $R_i$ is

$\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$

The IPF objective is to find a distribution satisfying all constraints while minimizing the I-divergence

$I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$

For probability distributions, this is the Kullback–Leibler divergence. It is nonnegative and equals zero only when its two arguments coincide. The solution is therefore the I-projection of $Q_0$ onto the intersection of the constraint sets.

For the constraint $R_i(y_i)$, the standard update is

$Q_k(x) = Q_{k-1}(x) \frac{R_i(y_i)}{Q_{k-1}(y_i)}.$

The multiplier is constant over configurations sharing the same value of $y_i$. Consequently,

$R_i(y_i)$0

after the update. A subsequent update can disturb a previously fitted marginal, so IPF cycles through all constraints. If the constraints are jointly consistent, ordinary IPF converges to the I-projection onto their intersection.

The multiplicative form follows from KL geometry. Given a current distribution $R_i(y_i)$1, minimizing $R_i(y_i)$2 subject to one affine marginal constraint produces an optimizer of the form

$R_i(y_i)$3

where $R_i(y_i)$4 indicates membership in the constrained marginal cell and $R_i(y_i)$5 is the target-to-current marginal ratio. IPF is thus a sequence of Bregman projections, or equivalently KL projections, onto affine constraint sets. This connects it to iterative scaling, information geometry, maximum entropy, matrix balancing, and entropy-regularized optimal transport (Loukas et al., 2022).

For categorical variables, IPF can be initialized with the uniform distribution. When the supplied marginal constraints are consistent, the limiting distribution is the unique maximum-entropy distribution satisfying them. Maximum entropy is equivalent to minimizing KL divergence from the uniform distribution:

$R_i(y_i)$6

where $R_i(y_i)$7 is uniform over the admissible microstates. The resulting distribution has an exponential-family representation, although the associated Lagrange multipliers need not be unique when the constraint system is redundant.

2. Matrix scaling, support, and convergence

For a nonnegative matrix $R_i(y_i)$8 with prescribed positive row and column marginals $R_i(y_i)$9 and $Y_i$0, IPF alternates row and column normalization:

$Y_i$1

followed by

$Y_i$2

The odd iterates have row marginals $Y_i$3, and the even iterates have column marginals $Y_i$4. When a matrix with the target marginals and support contained in that of $Y_i$5 exists, the sequence converges to the unique I-projection. If a full-support biproportional fitting exists, the limit has the form

$Y_i$6

where $Y_i$7 and $Y_i$8 are positive diagonal matrices. In this case convergence is at least geometric, and $Y_i$9 uniquely minimizes I-divergence over feasible matrices supported by $R_i$0 (Brossard et al., 2016).

Support is fundamental. IPF preserves zeros at every finite iteration:

$R_i$1

It cannot create mass outside the initial support. If target margins are feasible only after some entries disappear, the sequence can converge to a boundary point whose support is strictly smaller than the initial support. If no feasible matrix exists with support contained in that of $R_i$2, the full sequence may fail to converge while its even and odd subsequences converge separately to distinct limits. These limits have different marginal properties and are characterized through block-triangular support structures and I-projection problems.

For sparse matrices, feasibility can be tested by a Hall-type condition. If $R_i$3 denotes the columns adjacent in the support graph to a set of rows $R_i$4, then feasibility requires

$R_i$5

for every row subset $R_i$6. Violations identify blocking sets. The stochastic-matrix analysis of alternating normalization establishes convergence properties through backward products of stochastic matrices with diagonals bounded away from zero and controlled ratios of off-diagonal entries (Brossard et al., 2016).

The convergence behavior of general IPF procedures can also be quantified through the geometry of their constraint spaces. If the constraints arise from linear spaces $R_i$7, cyclic KL projections can converge exponentially when log-densities are uniformly bounded and the sum of the constraint spaces is closed. The contraction depends on the strong-convexity and smoothness ratio of the dual objective, the number of constraint blocks, and the condition number of the sum operator. For two spaces, the latter is governed by their Friedrichs angle: nearly orthogonal spaces yield faster convergence, while nearly parallel spaces produce slower convergence (Eckstein et al., 27 Feb 2025).

3. Bayesian networks, graphical models, and generalized iterative fitting

Ordinary IPF operates on an unrestricted joint table. Applying it directly to a Bayesian network can satisfy supplied marginals while violating the conditional independences encoded by the network’s directed acyclic graph. Extracting conditional probability tables from such a joint distribution and reconstructing the network need not reproduce the IPF table.

Extended IPF (E-IPFP) addresses this issue by adding a structural constraint. If the original Bayesian network factors as

$R_i$8

the structural update reconstructs the joint distribution from conditional tables extracted according to the original graph:

$R_i$9

E-IPFP cycles through the supplied probability constraints and this structural constraint. Under consistency, it converges to the minimum-I-divergence distribution satisfying both the supplied probability information and the original Bayesian-network factorization (Peng et al., 2012).

Decomposed IPF (D-IPFP) reduces computational cost by replacing a global E-IPFP computation with local E-IPFP problems on subnets. A constraint is local when it concerns one variable and zero or more of its parents. It can then be enforced by modifying only the associated conditional probability table. A non-local constraint involves variables associated with multiple conditional probability tables and is handled on a subnet containing the constrained variables and their external parents.

D-IPFP preserves the network structure and satisfies the constraints, but its normalized local updates are not generally global I-projections. It therefore usually produces a solution with larger KL divergence from the original network than E-IPFP. Its computational advantage arises because, for binary variables, global updates have approximate cost

$\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$0

whereas local updates operate on a subnet of size $\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$1, with approximate cost

$\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$2

In the reported 15-variable experiments, D-IPFP reduced runtime from $\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$3 seconds to $\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$4 seconds for four constraints and from $\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$5 seconds to $\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$6 seconds for sixteen constraints, while producing moderately larger divergences (Peng et al., 2012).

IPF also generalizes beyond ordinary undirected graphical models and Bayesian networks. Discrete chain factor graphs represent joint distributions as products of locally normalized conditional factors, encompassing undirected graphical models, Bayesian networks, chain graphs, and sigmoid belief networks. A generalized IPF algorithm updates one local potential at a time using expected sufficient statistics and coordinate-wise auxiliary functions. For incomplete data, it combines an EM-style posterior completion step with IPF-like potential updates. It reduces to standard EM for Bayesian networks and to classical IPF for undirected models (Wiegerinck et al., 2012).

For relational models without an overall effect, ordinary IPF, generalized iterative scaling, and improved iterative scaling can fail because the maximum-likelihood sufficient statistics are reproduced only up to a positive adjustment factor $\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$7. The generalized IPF$\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$8 algorithm fits

$\mathcal C_i=\{Q:Q(y_i)=R_i(y_i)\}.$9

while preserving the relational multiplicative structure. For probability models without an overall effect, $I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$0 is selected by bisection so that the fitted vector sums to one. This yields a generalized IPF procedure for curved exponential families (Klimova et al., 2013).

4. Applications in synthesis, balancing, and inference

IPF is widely used when individual-level or high-dimensional joint data are unavailable but lower-dimensional margins are known. In spatial microsimulation, survey microdata provide individual records and small-area aggregates provide constraints. IPF assigns a weight to every survey record for every geographic zone and repeatedly adjusts the weights to fit zonal totals.

The resulting weights are generally fractional. Integerization is required when the output must consist of whole simulated individuals, as in dynamic microsimulation or agent-based modeling. The “truncate, replicate, sample” method separates each IPF weight into an integer replication component and a fractional remainder:

$I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$1

Individuals are first replicated according to the integer part, and the remaining population is sampled using the fractional remainders. In the reported case study, TRS had total absolute error $I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$2, compared with $I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$3 for proportional-probability sampling and $I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$4 for simple rounding (Lovelace et al., 2013).

In mortality simulation, IPF disaggregates population or death tables across age, gender, smoker status, and region. It reconciles an initial mortality or population table with demographic margins, hazard-ratio information, and regional totals. Monte Carlo simulation is then applied separately to generate Poisson death realizations and uncertainty intervals. IPF determines the fitted demographic allocation; Monte Carlo generates stochastic outcomes conditional on that allocation (Nalmpatian et al., 12 Feb 2025).

IPF is also used to infer dynamic networks from an aggregate adjacency matrix and time-varying row and column marginals. Given an aggregate baseline $I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$5, the inferred network is

$I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$6

A biproportional Poisson model gives this procedure a maximum-likelihood interpretation:

$I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$7

The row and column marginals are sufficient statistics for the parameters $I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$8 and $I(Q\|Q_0) = \sum_x Q(x)\log\frac{Q(x)}{Q_0(x)}.$9, and IPF recovers their maximum-likelihood estimates when a finite scaling exists. Exact recovery of the realized network occurs only when the true network is itself a biproportional scaling of $Q_0$0. Estimation accuracy depends on total expected traffic and the connectivity of the weighted support graph, measured in part by the Fiedler eigenvalue. For sparse infeasible networks, ConvIPF uses maximum flow to identify blocking sets and adds a minimal set of edges needed to restore feasibility (Chang et al., 2024).

In geographically explicit population synthesis, IPF reweights survey records separately within census tracts so that demographic characteristics such as age, race or ethnicity, gender, income, and education match American Community Survey controls. Health attributes are inherited from the selected seed records after integerization. Experiments using zero-shot LLM-generated survey records show that IPF can reduce some upstream errors and amplify others; it cannot recover a category absent from the seed data because multiplicative updates preserve zero cells (Anderson et al., 23 Apr 2026).

5. Entropy-regularized optimal transport and Schrödinger bridges

In entropy-regularized optimal transport, IPF is commonly called Sinkhorn’s algorithm. Given a positive kernel $Q_0$1 and prescribed marginals $Q_0$2, the coupling has the form

$Q_0$3

where the scaling vectors are alternately updated so that row sums and column sums equal the target marginals. In continuous spaces, the same procedure uses positive scaling functions and KL projections onto the sets of measures with prescribed first and second marginals.

The Schrödinger bridge problem seeks the path measure closest in relative entropy to a reference process while satisfying endpoint constraints:

$Q_0$4

IPF alternately projects onto the initial- and terminal-marginal constraint sets. In the exact setting, it converges to the Schrödinger bridge under standard existence and absolute-continuity assumptions. For Brownian reference processes, the static endpoint problem is equivalent to entropy-regularized optimal transport with quadratic cost (Deligiannidis et al., 2021).

For compact metric spaces and Lipschitz costs, the IPF sequence is quantitatively stable with respect to perturbations of both marginals. If $Q_0$5 and $Q_0$6 are IPF iterates generated from $Q_0$7 and $Q_0$8, respectively, then

$Q_0$9

uniformly in $R_i(y_i)$0, with

$R_i(y_i)$1

The same bound applies to the limiting Schrödinger bridges. The proof uses uniform potential bounds, Hilbert–Birkhoff contraction, perturbation estimates for integral operators, and a geometric-series argument (Deligiannidis et al., 2021).

Recent Schrödinger-bridge work distinguishes IPF from Iterative Markovian Fitting (IMF). IPF alternates endpoint-marginal projections, whereas IMF alternates reciprocal and Markovian projections. A practical bidirectional method combines these operations into Iterative Proportional Markovian Fitting (IPMF), which incorporates endpoint corrections into IMF. In exact Gaussian settings, IPMF has parameter-level exponential convergence results; general neural implementations remain approximate (Kholkin et al., 2024).

For tree-structured entropy-regularized optimal transport and Wasserstein barycentres, IPF-based methods such as TreeDSB are contrasted with IMF-based TreeDSBM. The latter preserves prescribed marginals during idealized iterations and uses edgewise bridge matching, while TreeDSB is based on alternating IPF-like marginal projections. The comparison concerns computational and neural implementation properties rather than a new IPF convergence theorem (Howard et al., 20 Jun 2025).

6. Statistical interpretation, limitations, and methodological controversies

IPF has several distinct statistical interpretations, and they depend on the problem specification.

In maximum-entropy reconstruction, IPF selects the least-informative joint distribution compatible with supplied marginals. In contingency-table completion, it is the KL projection of a sample table onto tables with known population margins. In the $R_i(y_i)$2 case, alternating row and column scaling preserves the seed table’s odds ratio:

$R_i(y_i)$3

This provides the maximum-likelihood solution when the seed is a random sample from the same population whose margins are being imposed (Naszodi, 2023).

That interpretation does not automatically extend to counterfactual comparisons between different populations. In such settings, IPF preserves an odds-ratio-based association while changing the margins, but the seed table is not necessarily a sample from the counterfactual population. The NM method instead preserves an ordinal Liu–Lu association measure and is designed for counterfactual tables involving ordered categories. The distinction can change inferred historical trends in educational homogamy; survey evidence in the cited application supported the NM method over IPF for that counterfactual purpose (Naszodi, 2023, Naszodi, 2023).

Several limitations recur across applications:

• Constraint consistency: if no distribution satisfies all supplied constraints, exact IPF convergence is impossible. Inconsistent systems may oscillate or diverge.
• Support compatibility: zeros in the seed remain zero. Sampling zeros can therefore be mistaken for structural zeros and prevent feasible fitting.
• Seed dependence: when only partial margins are known, the seed table supplies association information. An independence seed imposes different higher-order structure from an observed or model-based seed.
• Dimensionality: direct joint-table implementations scale with the number of microstates, which grows as $R_i(y_i)$4.
• Marginal versus joint accuracy: exact fitted margins do not imply that unobserved interactions, conditional independences, or latent associations are correct.
• Model misspecification: in network inference, Poisson independence and stationarity assumptions may fail; in mortality synthesis, hazard-ratio transferability and regional uniformity are modeling assumptions.
• Computational iteration: IPF variants for graphical models can require repeated clique or block updates. For high-dimensional Gaussian graphical models, the Constrained Cholesky Approach has been proposed as a non-iterative alternative to IPF and related methods (Khare et al., 2024).

The appropriate use of IPF therefore depends on what is being estimated and what association structure is intended to remain invariant. It is well founded for KL projection, maximum-entropy fitting, population-table completion under a sample-to-population relationship, matrix scaling, and entropy-regularized transport. Its use for counterfactual populations, causal decompositions, or synthetic data should be accompanied by explicit justification of the seed, support, constraints, and preserved notion of association.