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

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}];


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...

• What if no matrix has the given row sums and column sums? The example here is a case in point. What result is expected? Some form of optimization? Commented Jul 27, 2023 at 17:15
• @Daniel Lichtblau: As I indicated below, there is a lower bound for the resulting matrix elements. They should be all non-negative and Total[refR]=Total[refC]. Commented Jul 27, 2023 at 22:59
• with your inputs refR and refC the condition Total[refR]=Total[refC] is not satisfied (and your iterative scheme oscillates indefinitely). Is it acceptable to pre-process refR and refC to ensure Total[refR]=Total[refC]?
– kglr
Commented Jul 27, 2023 at 23:31
• @kglr: Yes, you are right that the example in my question above has a typo, refC[[1]]=87 is the correct input. It is necessary to pre-process the condition Total[refR] = Total[refC]. The type of matrices I have has non-negative elements (i.e., input coefficients of a production function), meaning (1) if the original matrix element is positive, the resulting number cannot be 0 and must be positive but maybe very small, (2) if the original matrix element is equal to zero, it has to remain as zero. Zero elements cannot be positive after iterations. Commented Jul 27, 2023 at 23:52

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] &;

• Thank you very much. It works fine with the small matrix but when I try it with SeedRandom[123]; mat = RandomInteger[{0, 100}, {9, 9}]; refR = {150, 120, 140, 160, 140, 160, 150, 170, 150}; refC = {140, 170, 150, 140, 110, 150, 160, 150, 170}; it is very very slow. Is there anyway to speed it up? Commented Jul 28, 2023 at 10:55
• I did not understand what "M" refers to in gridDisplay[mat, refR, refC, msoln, "M", iter]. It seems to refer to imbalance, but why "M"? Commented Jul 28, 2023 at 11:18
• I got it, sorry...it is the notation for the final matrix with the iteration number as superscript. Commented Jul 28, 2023 at 11:20
• @TugrulTemel, re "anyway to speed it up", warp mat with N (ie use mat = N@ mat). Re M^8 it is meant to be the label for the for the form of the initital matrix at 8th iteration.
– kglr
Commented Jul 28, 2023 at 11:24
• Thanks for the tip. Now it is very very quick... Commented Jul 28, 2023 at 11:30

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

• Thanks a lot for this Code. It works fine. It is quite slow when matrix dimension goes over 5. Even with my very powerful computer, ComposeList is extremely slow and in fact I cannot get a solution because it runs for ever. Is there a way to speed up the operations? Commented Jul 27, 2023 at 14:33
• Thanks for the edit. It works just fine with the given example. For large matrices, it runs for ever with no result. Commented Jul 28, 2023 at 11:16
• When I wrap mat with N (suggested by kglr), your code is also very quickly converging. Thanks a lot for your efforts. Commented Jul 28, 2023 at 11:33
• The code you developed above is very very slow. Is there a way to speed it up? I am running a system with 75 by 75, and the code continuously runs with no output. It should not be that way, I. suppose. Commented Dec 3, 2023 at 19:08

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 ---

• Thanks for this efficient code. It works like charm. In fact, when I introduce a necessary condition, which is Total[refR]==Total[refC], then I get {1.65192*10^-19, {{23.8444, 32.6444, 18.5111}, {10.6111, 51.4111, 19.9778}, {52.5444, 65.9444, 6.51111}}} a very tiny error. Commented Jul 27, 2023 at 22:17
• When I tried it for a matrix (22,22), some of the resulting matrix entries become negative, which are not acceptable in the type of problems I address. It means that a condition needs to be imposed on the final solution matrix that all entries are non-negative. Is it possible to add this condition to your code? Commented Jul 27, 2023 at 22:39
• See edit for imposing constraints. Commented Jul 27, 2023 at 23:05