how to balance a matrix by proportional adjustment of rows and columns

Question

Given a matrix:

ClearAll[mat, refR, refC];
mat = {
   {21.2, 29, 19.9},
   {9, 48.8, 22.4},
   {49.8, 62.2, 7.8}
    };

I want to adjust the elements of mat to construct a new matrix satisfying the following row refR and column refC sums:

refR = {75, 82, 125};
refC = {87, 150, 45};

At 1st stage I adjust the rows by the following code:

mat1 = Table[(mat[[i, j]]/Total[mat, {2}][[i]])*refR[[i]], {i,3}, {j, 3}];

Then adjust the columns by:

mat2 = Table[(mat1[[j, i]]/Total[mat1][[i]])*refC[[i]], {j, 3}, {i, 3}];

Thereafter, I repeatedly carry out row and column adjustments:

mat3 = Table[(mat2[[i, j]]/Total[mat2, {2}][[i]])*refR[[i]], {i, 3}, {j, 3}];
mat4 = Table[  (mat3[[j, i]]/Total[mat3][[i]])*refC[[i]], {j, 3}, {i, 3}];
mat5 = Table[  (mat4[[i, j]]/Total[mat4, {2}][[i]])*refR[[i]], {i, 3}, {j, 3}];
mat6 = Table[  (mat5[[j, i]]/Total[mat5][[i]])*refC[[i]], {j, 3}, {i, 3}];
mat7 = Table[  (mat6[[i, j]]/Total[mat6, {2}][[i]])*refR[[i]], {i, 3}, {j, 3}];
mat8 = Table[  (mat7[[j, i]]/Total[mat7][[i]])*refC[[i]], {j, 3}, {i, 3}];
mat9 = Table[  (mat8[[i, j]]/Total[mat8, {2}][[i]])*refR[[i]], {i, 3}, {j, 3}];
mat10 = Table[  (mat9[[j, i]]/Total[mat9][[i]])*refC[[i]], {j, 3}, {i, 3}];
mat11 = Table[  (mat10[[i, j]]/Total[mat10, {2}][[i]])*refR[[i]], {i, 3}, {j, 3}];
mat12 = Table[  (mat11[[j, i]]/Total[mat11][[i]])*refC[[i]], {j, 3}, {i, 3}];
mat13 = Table[  (mat12[[i, j]]/Total[mat12, {2}][[i]])*refR[[i]], {i, 3}, {j, 3}];
mat14 = Table[  (mat13[[j, i]]/Total[mat13][[i]])*refC[[i]], {j, 3}, {i, 3}];

until the following conditions hold TRUE:

Total[mat14, {2}] == refR
Total[mat14] == refC

There must be a clever way of doing these iterations to balance the matrix mat to the new reference row and column sums as above.

Can a single function be created? Such as

matrixBalance[m_Matrix, refR_, refC_, mp_]:= 

where m_Matrix is the original unbalanced matrix, refR is the reference vector of row totals, refC is the reference vector of column totals, and mp is the number of matrix iterations. The outcome of this function should be the final balanced matrix satisfying refR and refC. In economics, this repeated iteration method is called RAS method and is widely used in the balancing of Input-Output matrix and/or Social Accounting Matrix.

Any improvement on the above code is appreciated...

Answer 1

Inputs:

mat = {{21.2, 29, 19.9}, {9, 48.8, 22.4}, {49.8, 62.2, 7.8}};
refR = {75, 82, 125};
refC = {87, 150, 45};

Iteration steps:

mnt = Map[Normalize[#, Total] &];

fr = refR  mnt @ # &;

fc = Transpose[refC  mnt @ Transpose @ #] &;

frc = fc @* fr;

Stopping criteria:

imbalance = Norm[Join[refC, refR] - Join[Total@#, Total[#, {2}]]] &;

tolerance = 10^-5;

unbalancedQ = imbalance[#] > tolerance &;

maxIterations = 10^3;

Initialize:

iter = 1;

msoln = N @ mat;

Iterate while current imbalance exceeds tolerance and iter is below maxiterations:

While[iter++ <= maxIterations && unbalancedQ[msoln = frc @ msoln]]

Display results:

gridDisplay[mat, refR, refC, msoln, "M", iter]

Helper functions for display:

augM = Grid[Join[Join[#, List /@  Total[#, {2}], 2], {Total@#}], 
    Dividers -> {{-2 -> Red}, {-2 -> Red}}] &;

gridDisplay = Grid[{{"\nmat", "refR\nrefC", 
      "\n" <> ToString[Superscript[#5, #6], StandardForm],  
      "\nimbalance", "\niterations"}, 
   {augM @ #, Column[{#2, #3}], augM @ #4, imbalance @ #4, iter}}, 
   Dividers -> All] &;

Answer 2

A starting point.

Clear[refR,refC,m,f,g];
refR = {x1, x2, x3};
refC = {y1, y2, y3};
m = Array[Subscript[a, ##] &, {3, 3}];
f = Times[#, refR] &@*Map[Normalize[#, Total] &];
g = Transpose[#] &@Times[#, refC] &@*Map[Normalize[#, Total] &]@
    Transpose[#] &;

f@m == Table[(m[[i, j]]/Total[m, {2}][[i]])*refR[[i]], {i, 3}, {j, 3}]
g@m == Table[(m[[j, i]]/Total[m][[i]])*refC[[i]], {j, 3}, {i, 3}]

True.

True.

mat = {{21.2, 29, 19.9}, {9, 48.8, 22.4}, {49.8, 62.2, 7.8}};
refR = {75, 82, 125};
refC = {85, 150, 45};
ComposeList[{f, g, f, g, f, g, f, g, f, g, f, g, f, g}, mat]
%== {mat, mat1, mat2, mat3, mat4, mat5, mat6, mat7, mat8, mat9, mat10, 
  mat11, mat12, mat13, mat14}

True.

Edit

Since refC = {87, 150, 45}; now the above code work.

mat = {{21.2, 29, 19.9}, {9, 48.8, 22.4}, {49.8, 62.2, 7.8}};
refR = {75, 82, 125};
refC = {87, 150, 45};
result=FixedPoint[g@*f, mat]
Total[result, {2}] == refR
Total[result] == refC

{{24.1287, 32.4684, 18.4029}, {9.81314, 52.342, 19.8448}, {53.0581, 65.1896, 6.75229}}

True.

True

Answer 3

This is actually an optimization problem. One thing that can be minimized is the sum-of-squares-of-discrepancies error.

matrixBalance[m_, refR_, refC_] /; 
  Dimensions[m] === {Length[refR], Length[refC]} := Module[
  {dims = Dimensions[m], x, mplus, eqC, eqR, vals, min, soln},
  mplus = Array[x, dims];
  eqC = Total[mplus + m] - refC;
  eqR = Total[mplus + m, {2}] - refR;
  vals = Join[eqR, eqC];
  {min, soln} = FindMinimum[vals . vals, Flatten[mplus]];
  {min, m + mplus /. soln}
  ]

This example:

mat = {{21.2, 29, 19.9}, {9, 48.8, 22.4}, {49.8, 62.2, 7.8}};
refR = {75, 82, 125};
refC = {85, 150, 45};
mb = matrixBalance[mat, refR, refC]

(* Out[187]= {0.666667, {{23.2889, 32.7556, 18.6222}, {10.0556, 51.5222, 
   20.0889}, {51.9889, 66.0556, 6.62222}}} *)

This is in fact close to mat14 from the question. But that has sos error of 1.40707 instead of 2/3.

--- edit ---

If there are constraints on the result e.g. nonnegativity of matrix elements (as noted in a comment), one can use this optimization instead.

{min, soln} = 
  FindMinimum[{vals . vals, Thread[Flatten[mplus + m] >= 0]}, 
   Flatten[mplus]];

In general, just place the constraints in a list after the objective function. Since these particular constraints are linear, FindMinimum will likely use quadratic programming under the hood so it should still be fairly efficient. And one can minimize the 1-norm instead of Euclidean, and then cast it as a linear program. This is probably overkill but could matter when dimensions are large.

--- end edit ---