Iterative proportional fitting

{{Short description|Estimates values in an N-dimensional matrix}} The '''iterative proportional fitting procedure''' ('''IPF''' or '''IPFP''', also known as '''biproportional fitting''' or '''biproportion''' in statistics or economics (input-output analysis, etc.), '''RAS algorithm'''<ref>{{cite journal |last=Bacharach|first=M. |year=1965 |title=Estimating Nonnegative Matrices from Marginal Data |journal=International Economic Review |volume=6|pages=294–310 |doi=10.2307/2525582 |jstor=2525582 |issue=3 |publisher=Blackwell Publishing }}</ref> in economics, '''raking''' in survey statistics, and '''matrix scaling''' in computer science) is the operation of finding the fitted matrix <math>X</math> which is the closest to an initial matrix <math>Z</math> but with the row and column totals of a target matrix <math>Y</math> (which provides the constraints of the problem; the interior of <math>Y</math> is unknown). The fitted matrix being of the form <math>X=PZQ</math>, where <math>P</math> and <math>Q</math> are diagonal matrices such that <math>X</math> has the margins (row and column sums) of <math>Y</math>. Some algorithms can be chosen to perform biproportion. We have also the entropy maximization,<ref>Jaynes E.T. (1957) Information theory and statistical mechanics, Physical Review, 106: 620-30.</ref><ref>Wilson A.G. (1970) Entropy in urban and regional modelling. London: Pion LTD, Monograph in spatial and environmental systems analysis.</ref> information loss minimization (or cross-entropy)<ref>Kullback S. & Leibler R.A. (1951) On information and sufficiency, Annals of Mathematics and Statistics, 22 (1951) 79-86.</ref> or RAS which consists of factoring the matrix rows to match the specified row totals, then factoring its columns to match the specified column totals; each step usually disturbs the previous step's match, so these steps are repeated in cycles, re-adjusting the rows and columns in turn, until all specified marginal totals are satisfactorily approximated. However, all algorithms give the same solution.<ref>{{cite journal |last=de Mesnard|first=L.|year=1994 |title=Unicity of Biproportion |journal=SIAM Journal on Matrix Analysis and Applications |volume=15 |issue=2 |pages=490–495 |doi=10.1137/S0895479891222507 }}https://www.researchgate.net/publication/243095013_Unicity_of_Biproportion</ref> In three- or more-dimensional cases, adjustment steps are applied for the marginals of each dimension in turn, the steps likewise repeated in cycles.

== History ==

IPF has been "re-invented" many times, the earliest by Kruithof in 1937 <ref>{{cite journal |last1=Kruithof |first1=J |title=Telefoonverkeersrekening (Calculation of telephone traffic) |journal=De Ingenieur |date=February 1937 |volume=52 |issue=8 |pages=E15–E25 |url=https://resolver.kb.nl/resolve?urn=dts:2980078:mpeg21}}</ref> in relation to telephone traffic ("Kruithof’s double factor method"), Deming and Stephan in 1940<ref>{{cite journal |last1=Deming |first1=W. E.|author-link1=W. Edwards Deming |last2=Stephan |first2=F. F. |year=1940 |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 |issue=4 |pages=427–444 |mr=3527 |doi=10.1214/aoms/1177731829 |doi-access=free}}</ref> for adjusting census crosstabulations, and G.V. Sheleikhovskii for traffic as reported by Bregman.<ref>Lamond, B. and Stewart, N.F. (1981) Bregman's balancing method. Transportation Research 15B, 239-248.</ref> (Deming and Stephan proposed IPFP as an algorithm leading to a minimizer of the Pearson X-squared statistic, which Stephan later reported it ''does not'').<ref>{{cite journal |last=Stephan |first=F. F.|year=1942 |title=Iterative method of adjusting frequency tables when expected margins are known |journal=Annals of Mathematical Statistics |volume=13 |issue=2 |pages=166–178 |mr=6674 | zbl = 0060.31505 |doi=10.1214/aoms/1177731604 |doi-access=free}}</ref> Early proofs of uniqueness and convergence came from Sinkhorn (1964),<ref>Sinkhorn, Richard (1964). “A Relationship Between Arbitrary Positive Matrices and Doubly Stochastic Matrices”. In: Annals of Mathematical Statistics 35.2, pp. 876–879.</ref> Bacharach (1965),<ref>Bacharach, Michael (1965). “Estimating Nonnegative Matrices from Marginal Data”. In: International Economic Review 6.3, pp. 294–310.</ref> Bishop (1967),<ref>Bishop, Y. M. M. (1967). “Multidimensional contingency tables: cell estimates”. PhD thesis. Harvard University.</ref> and Fienberg (1970).<ref>{{cite journal |last=Fienberg |first=S. E.|author-link=Stephen Fienberg |year=1970 |title=An Iterative Procedure for Estimation in Contingency Tables |journal=Annals of Mathematical Statistics |volume=41 |issue=3 |pages=907–917 |mr=266394 | zbl = 0198.23401 | jstor = 2239244 |doi=10.1214/aoms/1177696968 |doi-access=free}}</ref> Bishop's proof that IPFP finds the maximum likelihood estimator for any number of dimensions extended a 1959 proof by Brown for 2x2x2... cases. Fienberg's proof by differential geometry exploits the method's constant crossproduct ratios, for strictly positive tables. Csiszár (1975).<ref>{{cite journal |last=Csiszár |first=I.|author-link=Imre Csiszár |year=1975 |title=''I''-Divergence of Probability Distributions and Minimization Problems |journal=Annals of Probability |volume=3 |issue=1 |pages=146–158 |mr=365798 | zbl = 0318.60013 | jstor = 2959270 |doi=10.1214/aop/1176996454 |doi-access=free}}</ref> found necessary and sufficient conditions for general tables having zero entries. Pukelsheim and Simeone (2009) <ref>{{cite web |title=On the Iterative Proportional Fitting Procedure: Structure of Accumulation Points and L1-Error Analysis |url=http://opus.bibliothek.uni-augsburg.de/volltexte/2009/1368/ |publisher= Pukelsheim, F. and Simeone, B. |access-date=2009-06-28}}</ref> give further results on convergence and error behavior.

An exhaustive treatment of the algorithm and its mathematical foundations can be found in the book of Bishop et al. (1975).<ref>{{cite book |title=Discrete Multivariate Analysis: Theory and Practice |last1=Bishop |first1=Y. M. M. |author1-link=Yvonne Bishop |first2=S. E. |last2=Fienberg |author-link2=Stephen Fienberg |first3=P. W. |last3=Holland |year=1975 |publisher=MIT Press |isbn=978-0-262-02113-5 |mr=381130 |url=https://archive.org/details/discretemultivar00bish }}</ref> Idel (2016)<ref>Martin Idel (2016) A review of matrix scaling and Sinkhorn’s normal form for matrices and positive maps arXiv preprint https://arxiv.org/pdf/1609.06349.pdf</ref> gives a more recent survey.

Other general algorithms can be modified to yield the same limit as the IPFP, for instance the Newton–Raphson method and the EM algorithm. In most cases, IPFP is preferred due to its computational speed, low storage requirements, numerical stability and algebraic simplicity.

Applications of IPFP have grown to include trip distribution models, Fratar or Furness and other applications in transportation planning (Lamond and Stewart), survey weighting, synthesis of cross-classified demographic data, adjusting input–output models in economics, estimating expected quasi-independent contingency tables, biproportional apportionment systems of political representation, and for a preconditioner in linear algebra.<ref>Bradley, A.M. (2010) Algorithms for the equilibration of matrices and their application to limited-memory quasi-newton methods. Ph.D. thesis, Institute for Computational and Mathematical Engineering, Stanford University, 2010</ref>

== Biproportion ==

Biproportion, whatever the algorithm used to solve it, is the following concept: <math>Z</math>, matrix <math>Y</math> and matrix <math>X</math> are known real nonnegative matrices of dimension <math>n,m</math>; the interior of <math>Y</math> is unknown and <math>X</math> is searched such that <math>X</math> has the same margins than <math>Y</math>, i.e. <math>Xs=Ys</math> and <math>s'X=s'Y</math> (<math>s</math> being the sum vector), and such that <math>X</math> is close to <math>Z</math> following a given criterion, the fitted matrix being of the form <math>X=K(Z,Y)=PZQ</math>, where <math>P</math> and <math>Q</math> are diagonal matrices.

<math>min \sum_i\sum_j x_{ij}\log(x_{ij}/z_{ij})</math> s.t.<math>\sum_j x_{ij}=y_{i.}</math>, ∀<math>i</math> and <math>\sum_i x_{ij}=y_{.j}</math>, ∀<math>j</math>. The Lagrangian is <math>L=\sum_i\sum_j x_{ij}\log(x_{ij}/z_{ij})-\sum_i p_i(y_{i.}-\sum_j x_{ij})-\sum_jq_j(y_{.j}-\sum_i x_{ij})</math>.

Thus <math>x_{ij}=z_{ij} \exp-(1+p_i+q_j)</math>, for ∀<math>i,j</math>,

which, after posing <math>P_i=\exp-(1+p_i)</math> and <math>Q_j=\exp-q_j</math>, yields

<math>x_{ij}=P_i z_{ij} Q_j</math>, ∀<math>i,j</math>, i.e., <math>X=PZQ</math>,

with <math>P_i=y_{i.}(\sum_j z_{ij}Q_j)^{-1}</math>, ∀<math>i</math> and <math>Q_j=y_{.j}(\sum_i z_{ij}P_i)^{-1}</math>, ∀<math>j</math>. <math>P_i</math> and <math>Q_j</math> form a system that can be solve iteratively:

<math>P_i^{(t+1)}=y_{i.}(\sum_j z_{ij}Q_j^{(t)})^{-1}</math>, ∀<math>i</math> and <math>Q_j^{(t+1)}=y_{.j}(\sum_i z_{ij}P_i^{(t+1)})^{-1}</math>, ∀<math>j</math>.

The solution <math>X</math> is independent of the initialization chosen (i.e., we can begin by <math>q_j^{(0)}=1</math>, ∀<math>j</math> or by <math>p_i^{(0)}=1</math>, ∀<math>i</math>. If the matrix <math>Z</math> is “indecomposable”, then this process has a unique fixed-point because it is deduced from a program where the function is a convex and continuously derivable function defined on a compact set. In some cases the solution may not exist: see de Mesnard's example cited by Miller and Blair (Miller R.E. & Blair P.D. (2009) Input-output analysis: Foundations and Extensions, Second edition, Cambridge (UK): Cambridge University Press, p.&nbsp;335-336 (freely available)).

Some properties (see de Mesnard (1994)):

Lack of information: if <math>Z</math> brings no information, i.e., <math>z_{ij}=z</math>, ∀<math>i,j</math> then <math>X=PQ</math>.

Idempotency: <math>X=K(Z,Y)=Z</math> if <math>Y</math> has the same margins than <math>Z</math>.

Composition of biproportions: <math>K(K(Z,Y_1),Y_2=K(Z,Y_2)</math>; <math>K(...K(Z,Y_1),Y_2)...Z_N)=K(Z,Y_N)</math>.

Zeros: a zero in <math>Z</math> is projected as a zero in <math>X</math>. Thus, a bloc-diagonal matrix is projected as a bloc-diagonal matrix and a triangular matrix is projected as a triangular matrix.

Theorem of separable modifications: if <math>Z</math> is premutiplied by a diagonal matrix and/or postmultiplied by a diagonal matrix, then the solution is unchanged.

Theorem of "unicity": If <math>K^{q}</math> is any non-specified algorithm, with <math>\hat{X}=K^{q}(Z,Y)=UZV</math>, <math>U</math> and <math>V</math> being unknown, then <math>U</math> and <math>V</math> can always be changed into the standard form of <math>P</math> and <math>Q</math>. The demonstrations calls some above properties, particularly the Theorem of separable modifications and the composition of biproportions.

== Algorithm 1 (classical IPF) ==

Given a two-way (''I'' &times; ''J'')-table <math>x_{ij}</math>, we wish to estimate a new table <math>\hat{m}_{ij} = a_i b_j x_{ij}</math> for all ''i'' and ''j'' such that the marginals satisfy <math>\sum_j \hat{m}_{ij}\ = u_{i},</math> and <math>\sum_i \hat{m}_{ij}\ = v_{j}</math>.

Choose initial values <math>\hat{m}_{ij}^{(0)} := x_{ij}</math>, and for <math>\eta \geq 1</math> set

: <math>\hat{m}_{ij}^{(2\eta - 1)} = \frac{\hat{m}_{ij}^{(2\eta-2)}u_{i}}{\sum_{k=1}^J \hat{m}_{ik}^{(2\eta-2)}}</math>

: <math>\hat{m}_{ij}^{(2\eta)} = \frac{\hat{m}_{ij}^{(2\eta-1)}v_{j}}{\sum_{k=1}^I \hat{m}_{kj}^{(2\eta-1)}}.</math>

Repeat these steps until row and column totals are sufficiently close to u and v.

Notes: * For the RAS form of the algorithm, define the diagonalization operator <math>diag: \mathbb{R}^k \longrightarrow \mathbb{R}^{k\times k}</math>, which produces a (diagonal) matrix with its input vector on the main diagonal and zero elsewhere. Then, for each row adjustment, let <math>R^{\eta}=diag(\frac{u_i}{\sum_j m_{ij}^{(2\eta-2)}})</math>, from which <math>M^{2\eta - 1}=R^{\eta}M^{2\eta-2}</math>. Similarly each column adjustment's <math>S^{\eta}=diag(\frac{v_i}{\sum_i m_{ij}^{(2\eta-1)}})</math>, from which <math>M^{2\eta}=M^{2\eta-1}S^{\eta}</math>. Reducing the operations to the necessary ones, it can easily be seen that RAS does the same as classical IPF. In practice, one would not implement actual matrix multiplication with the whole R and S matrices; the RAS form is more a notational than computational convenience.

== Algorithm 2 (factor estimation) ==

Assume the same setting as in the classical IPFP. Alternatively, we can estimate the row and column factors separately: Choose initial values <math>\hat{b}_j^{(0)} := 1</math>, and for <math>\eta \geq 1</math> set

: <math>\hat{a}_i^{(\eta)} = \frac{u_{i}}{\sum_j\ x_{ij}\hat{b}_j^{(\eta-1)}},</math>

: <math>\hat{b}_j^{(\eta)} = \frac{v_{j}}{\sum_i\ x_{ij}\hat{a}_i^{(\eta)}}</math>

Repeat these steps until successive changes of a and b are sufficiently negligible (indicating the resulting row- and column-sums are close to u and v).

Finally, the result matrix is <math>\hat{m}_{ij} = \hat{a}_i^{(\eta)}\hat{b}_j^{(\eta)}x_{ij}</math>

Notes:

* The two variants of the algorithm are mathematically equivalent, as can be seen by formal induction. With factor estimation, it is not necessary to actually compute each cycle's <math>\hat{m}_{ij}^{(\eta)}</math>. * The factorization is not unique, since it is <math>m_{ij} = a_i b_j x_{ij} = (\gamma a_i)(\frac{1}{\gamma}b_j)x_{ij}</math> for all <math>\gamma > 0</math>.

== Discussion ==

The vaguely demanded 'similarity' between M and X can be explained as follows: IPFP (and thus RAS) maintains the crossproduct ratios, i.e.

: <math>\frac{m^{(\eta)}_{ij}m^{(\eta)}_{hk}}{m^{(\eta)}_{ik}m^{(\eta)}_{hj}} = \frac{x_{ij}x_{hk}}{x_{ik}x_{hj}}\ \forall\ \eta \geq 0\text{ and }i\neq h,\quad j\neq k</math>

since <math>m^{(\eta)}_{ij} = a_i^{(\eta)}b_j^{(\eta)}x_{ij}.</math>

This property is sometimes called '''structure conservation''' and directly leads to the geometrical interpretation of contingency tables and the proof of convergence in the seminal paper of Fienberg (1970).

Direct factor estimation (algorithm 2) is generally the more efficient way to solve IPF: Whereas a form of the classical IPFP needs

: <math>IJ(2+J) + IJ(2+I) = I^2J + IJ^2 + 4IJ \, </math>

elementary operations in each iteration step (including a row and a column fitting step), factor estimation needs only

: <math>I(1+J) + J(1+I) = 2IJ + I + J \, </math>

operations being at least one order in magnitude faster than classical IPFP.

IPFP can be used to estimate expected quasi-independent (incomplete) contingency tables, with <math>u_i = x_{i+}, v_j = x_{+j}</math>, and <math>m^{0}_{ij}=1</math> for included cells and <math>m^{0}_{ij}=0</math> for excluded cells. For fully independent (complete) contingency tables, estimation with IPFP concludes exactly in one cycle.

=== Comparison with the NM-method === Similar to the IPF, the NM-method is also an operation of finding a matrix <math>X</math> which is the “closest” to matrix <math>Z</math> (<math> Z \in \mathbb{N}^{n \times m } </math>) while its row totals and column totals are identical to those of a target matrix <math>Y</math> <math>( Y \in \mathbb{N}^{n \times m }) </math>.

However, there are '''differences between the NM-method and the IPF'''. For instance, the '''NM-method defines closeness''' of matrices of the same size differently from the IPF.<ref name=NM2021>{{Cite journal |last1=Naszodi |first1=A. |last2=Mendonca |first2=F. |year=2021 |title=A new method for identifying the role of marital preferences at shaping marriage patterns |journal=Journal of Demographic Economics |volume=1 |issue=1 |pages=1–27 |doi=10.1017/dem.2021.1 |doi-access=free}}</ref> Also, the NM-method was developed to solve for matrix <math>X</math> in problems, where '''matrix''' <math>\boldsymbol{Z}</math> '''is not a sample''' from the population characterized by the row totals and column totals of matrix <math>Y</math>, but represents '''another population'''.<ref name=NM2021/> In contrast, '''matrix''' <math>\boldsymbol{Z}</math> '''is a sample''' from this population in problems where the IPF is applied as the '''maximum likelihood estimator'''.

Macdonald (2023)<ref name=KM2023>{{cite journal |last1=Macdonald |first1=K. |year=2023 |title=The marginal adjustment of mobility tables, revisited |journal=OSF |pages=1–19 |url=https://osf.io/z4u2d/download}}</ref> is at ease with the conclusion by Naszodi (2023)<ref name=N_IPF_NM_2023 >{{cite arXiv |last=Naszodi |first=A. |year=2023 |title=The iterative proportional fitting algorithm and the NM-method: solutions for two different sets of problems |class=econ.GN |eprint=2303.05515}}</ref> that the IPF is suitable for sampling correction tasks, but not for generation of counterfactuals. Similarly to Naszodi, Macdonald also questions whether the row and column proportional transformations of the IPF preserve the structure of association within a contingency table that allows us to study social mobility.

== Existence and uniqueness of MLEs ==

Necessary and sufficient conditions for the existence and uniqueness of MLEs are complicated in the general case (see<ref>{{cite book |title=The Analysis of Frequency Data |last=Haberman |first=S. J. |year=1974 |publisher=Univ. Chicago Press |isbn=978-0-226-31184-5 |url=https://archive.org/details/analysisoffreque0000habe }}</ref>), but sufficient conditions for 2-dimensional tables are simple:

* the marginals of the observed table do not vanish (that is, <math>x_{i+} > 0,\ x_{+j} > 0</math>) and * the observed table is inseparable (i.e. the table does not permute to a block-diagonal shape).

If unique MLEs exist, IPFP exhibits linear convergence in the worst case (Fienberg 1970), but exponential convergence has also been observed (Pukelsheim and Simeone 2009). If a direct estimator (i.e. a closed form of <math>(\hat{m}_{ij})</math>) exists, IPFP converges after 2 iterations. If unique MLEs do not exist, IPFP converges toward the so-called ''extended MLEs'' by design (Haberman 1974), but convergence may be arbitrarily slow and often computationally infeasible.

If all observed values are strictly positive, existence and uniqueness of MLEs and therefore convergence is ensured.

== Example ==

Consider the following table, given with the row- and column-sums and targets.

{| class="wikitable" style="margin:1em auto;" ! || 1 || 2 || 3 || 4 || TOTAL || TARGET |- | 1 || 40 || 30 || 20 || 10 || 100 || 150 |- | 2 || 35 || 50 || 100 || 75 || 260 || 300 |- | 3 || 30 || 80 || 70 || 120 || 300 || 400 |- | 4 || 20 || 30 || 40 || 50 || 140 || 150 |- | TOTAL || 125 || 190 || 230 || 255 || 800 || |- | TARGET|| 200 || 300 || 400 || 100 || || 1000 |- |}

For executing the classical IPFP, we first adjust the rows:

{| class="wikitable" style="margin:1em auto;" ! || 1 || 2 || 3 || 4 || TOTAL || TARGET |- | 1 || 60.00 || 45.00 || 30.00 || 15.00 || 150.00 || 150 |- | 2 || 40.38 || 57.69 || 115.38 || 86.54 || 300.00 || 300 |- | 3 || 40.00 || 106.67 || 93.33 || 160.00 || 400.00 || 400 |- | 4 || 21.43 || 32.14 || 42.86 || 53.57 || 150.00 || 150 |- | TOTAL || 161.81 || 241.50 || 281.58 || 315.11 || 1000.00 || |- | TARGET|| 200 || 300 || 400 || 100 || || 1000 |- |}

The first step exactly matched row sums, but not the column sums. Next we adjust the columns:

{| class="wikitable" style="margin:1em auto;" ! || 1 || 2 || 3 || 4 || TOTAL || TARGET |- | 1 || 74.16 || 55.90 || 42.62 || 4.76 || 177.44 || 150 |- | 2 || 49.92 || 71.67 || 163.91 || 27.46 || 312.96 || 300 |- | 3 || 49.44 || 132.50 || 132.59 || 50.78 || 365.31 || 400 |- | 4 || 26.49 || 39.93 || 60.88 || 17.00 || 144.30 || 150 |- | TOTAL || 200.00 || 300.00 || 400.00 || 100.00 || 1000.00 || |- | TARGET|| 200 || 300 || 400 || 100 || || 1000 |- |}

Now the column sums exactly match their targets, but the row sums no longer match theirs. After completing three cycles, each with a row adjustment and a column adjustment, we get a closer approximation:

{| class="wikitable" style="margin:1em auto;" ! || 1 || 2 || 3 || 4 || TOTAL || TARGET |- | 1 || 64.61 || 46.28 || 35.42 || 3.83 || 150.13 || 150 |- | 2 || 49.95 || 68.15 || 156.49 || 25.37 || 299.96 || 300 |- | 3 || 56.70 || 144.40 || 145.06 || 53.76 || 399.92 || 400 |- | 4 || 28.74 || 41.18 || 63.03 || 17.03 || 149.99 || 150 |- | TOTAL || 200.00 || 300.00 || 400.00 || 100.00 || 1000.00 || |- | TARGET|| 200 || 300 || 400 || 100 || || 1000 |- |}

== Implementation ==

The R package '''mipfp''' (currently in version 3.2) provides a multi-dimensional implementation of the traditional iterative proportional fitting procedure.<ref>{{cite web|last1=Barthélemy|first1=Johan|last2=Suesse|first2=Thomas|title=mipfp: Multidimensional Iterative Proportional Fitting|url=https://cran.r-project.org/web/packages/mipfp/index.html|publisher=CRAN|access-date=23 February 2015}}</ref> The package allows the updating of a ''N''-dimensional array with respect to given target marginal distributions (which, in turn can be multi-dimensional).

Python has an equivalent package, '''ipfn'''<ref>{{cite web|title=ipfn: pip|url=https://pypi.org/project/ipfn/}}</ref><ref>{{cite web|title=ipfn: github|website=GitHub |url=https://github.com/Dirguis/ipfn}}</ref> that can be installed via pip. The package supports numpy and pandas input objects.

== See also ==

* Data cleansing * Data editing * NM-method * Triangulation (social science) for quantitative and qualitative study data enhancement.

== References == {{Reflist}}

{{DEFAULTSORT:Iterative Proportional Fitting}} Category:Contingency table Category:Statistical algorithms