Iterative proportional fitting
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 in economics, raking in survey statistics, and matrix scaling in computer science) is the operation of finding the fitted matrix formula_1 which is the closest to an initial matrix formula_2 but with the row and column totals of a target matrix formula_3 (which provides the constraints of the problem; the interior of formula_3 is unknown). The fitted matrix being of the form formula_5, where formula_6 and formula_7 are diagonal matrices such that  formula_1 has the margins (row and column sums) of formula_3. Some algorithms can be chosen to perform biproportion. We have also the entropy maximization, information loss minimization (or cross-entropy) 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. 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 in relation to telephone traffic ("Kruithof’s double factor method"), Deming and Stephan in 1940 for adjusting census crosstabulations, and G.V. Sheleikhovskii for traffic as reported by Bregman. (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". Early proofs of uniqueness and convergence came from Sinkhorn (1964), Bacharach (1965), Bishop (1967), and  Fienberg (1970). 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). found necessary and sufficient conditions for general tables having zero entries. Pukelsheim and Simeone (2009) An exhaustive treatment of the algorithm and its mathematical foundations can be found in the book of Bishop et al. (1975). Idel (2016) 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. Biproportion. Biproportion, whatever the algorithm used to solve it, is the following concept: formula_2 , matrix formula_3 and matrix formula_1 are known real nonnegative matrices of dimension formula_13; the interior of formula_3 is unknown and formula_1 is searched such that formula_1 has the same margins than formula_3, i.e. formula_18 and formula_19 (formula_20 being the sum vector, and such that formula_1 is closed to formula_2 following a given criterion, the fitted matrix being of the form formula_23, where formula_6 and formula_7 are diagonal matrices. formula_26 s.t.formula_27, ∀formula_28 and formula_29, ∀formula_30. The Lagrangian is formula_31. Thus formula_32, for ∀formula_33, which, after posing formula_34 and formula_35, yields formula_36, ∀formula_33, i.e., formula_5, with formula_39, ∀formula_28 and formula_41, ∀formula_30. formula_43  and formula_44 form a system that can be solve iteratively: formula_45, ∀formula_28 and formula_47, ∀formula_30. The solution formula_1 is independent of the initialization chosen (i.e., we can begin by formula_50, ∀formula_30 or byformula_52, ∀formula_28. If the matrix formula_2 is “indecomposable”, then this process has a unique fixed-point because it is deduced from program 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. 335-336 (freely available)). Some properties (see de Mesnard (1994)): Lack of information: if formula_2 brings no information, i.e., formula_56, ∀formula_33 then formula_58. Idempotency: formula_59 if formula_3 has the same margins than formula_2. Composition of biproportions: formula_62; formula_63. Zeros: a zero in formula_2 is projected as a zero in formula_1. 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 formula_2 is premutiplied by a diagonal matrix and/or postmultiplied by a diagonal matrix, then the  solution is unchanged. Theorem of "unicity": If formula_67 is any non-specified algorithm, with formula_68, formula_69 and formula_70 being unknown, then formula_69 and formula_70 can always be changed into the standard form of formula_6 and formula_7. 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" × "J")-table formula_75, we wish to estimate a new table formula_76 for all "i" and "j" such that the marginals satisfy formula_77 and formula_78. Choose initial values formula_79, and for formula_80 set Repeat these steps until row and  column totals are sufficiently close to u and v. Notes: 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 formula_88, and for formula_80 set 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 formula_92 Notes: Discussion. The vaguely demanded 'similarity' between M and X can be explained as follows: IPFP (and thus RAS) maintains the crossproduct ratios, i.e. since formula_97 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 elementary operations in each iteration step (including a row and a column fitting step), factor estimation needs only 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 formula_100, and formula_101 for included cells and  formula_102 for excluded cells. For fully independent (complete) contingency tables, estimation with IPFP concludes exactly in one cycle. Existence and uniqueness of MLEs. Necessary and sufficient conditions for the existence and uniqueness of MLEs are complicated in the general case (see), but sufficient conditions for 2-dimensional tables are simple: 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 formula_104) 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. For executing the classical IPFP, we first adjust the rows: The first step exactly matched row sums, but not the column sums. Next we adjust the columns: 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: Implementation. The R package mipfp (currently in version 3.1) provides a multi-dimensional implementation of the traditional iterative proportional fitting procedure. 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 that can be installed via pip. The package supports numpy and pandas input objects.