# Find a float matrix with given row and column sums and predetermined arbitrary values

Sorry for the naivety and poor exposition of the question but I've been banging my head against this one for way too long and, as you'll be able to tell, math is not my forte.

My problem is similar to this one, with some additional information.

The simple version of it can be stated like so: given an empty matrix with known row and column equal sums, find a matrix that satisfies the aforementioned sum constraints. Since there is no constraints for the numbers to be integers, a possible solution can be found using 𝑟𝑖𝑐𝑗/𝑇 in position (𝑖,𝑗), where 𝑟𝑖 is the sum for row 𝑖, 𝑐𝑗 the sum for column 𝑗 and 𝑇 the total. For instance:

$$\begin{array}{cc|c} a & b & 2 \\ c & d & 4 \\ \hline 5 & 1 & 6 \\ \end{array}$$

gives us:

$$\begin{array}{cc|c} 1.667 & 0.333 & 2 \\ 3.333 & 0.667 & 4 \\ \hline 5 & 1 & 6 \\ \end{array}$$

So far, so good. Now my actual problem is more complex as arbitrary values are set to 0 in the starting matrix, for instance:

$$\begin{array}{ccc|c} a & 0 & b & 4 \\ 0 & c & d & 8 \\ \hline 4 & 2 & 6 & 12\\ \end{array}$$

This is where I'm stuck. I can get the columns to add up but then the rows stop adding up, or vice-versa, doing e.g. something like: $$r_ic_j/(T-x)$$ where $$x$$ is the column sum corresponding to the 0 value in the row. Note that I would ideally keep with that kind of proportional approach as it fits the real-world scenario behind it best. I've looked into solving this using linear programming but to no avail.

Let me know if anything needs clarifying.

• You just have a system of linear equations. The non-zero matrix elements are the unknowns and you have one equation for each row, and one for each column. – saulspatz Aug 7 '20 at 13:18
• Take a look at this. Just solve over the reals, rather than $\Bbb F_2$. – Rodrigo de Azevedo Aug 7 '20 at 14:09
• I assume you know about the iterative proportional fitting algorithm, so I assume the issue is preserving those zero values. The column constraints immediately yield $\,(a=4,\,c=2)\,$ after which the first row constraint yields $\,(b=4-a=0)\,$ and the second yields $\,(d=8-c=6)\,$ – greg Aug 7 '20 at 15:20

The standard iterative proportional fitting (IPF) algorithm in Julia/Matlab looks like

function ipf(x,y)
u,v = x-x.+1, y-y.+1
A = u*v'
for k = 1:2000
B = A .* (x*v') ./ ((A*v)*v')
A = B .* (u*y') ./ (u*(u'*B))
if norm(A-B) < 1e-6; break; end
end
return A
end


where the $$(x,y)$$ vectors contain the marginal totals,
$$(u,v)$$ are all-ones vectors the same size as $$(x,y)$$, and
$$A$$ is an initial (unfit) table; let's use $$A=uv^T\,$$ (i.e. the all-ones matrix).

In matrix notation, the algorithm is \eqalign{ A_0 &= uv^T \\ B_{k+1} &= \left(\frac{A_k\odot xv^T}{A_k\cdot vv^T}\right) \\ A_{k+1} &= \left(\frac{uy^T\odot B_{k+1}}{uu^T\cdot B_{k+1}}\right) \\ }

By including a $$\{0,\!{\tt1}\}$$ pattern matrix to specify the location of the zero and non-zero elements, e.g. $$P = \pmatrix{1&0&1\\0&1&1}$$ the algorithm can be modified to preserve those zero elements

function ipfm(P,x,y)
u,v = x-x.+1, y-y.+1
A = u*v'
for k = 1:2000
B = A .* (x*v') ./ ((A*v)*v')
A = B .* (u*y') ./ (u*(u'*B)) .* P
if norm(A-B) < 1e-6; break; end
end
return A
end


or in matrix notation \eqalign{ A_0 &= uv^T \\ B_{k+1} &= \left(\frac{A_k\odot xv^T}{A_k\cdot vv^T}\right) \\ A_{k+1} &= \left(\frac{uy^T\odot B_{k+1}}{uu^T\cdot B_{k+1}}\right) \odot P \\ }

For the current problem, the algorithms yield

x,y,P = [4;8],  [4;2;6],  [1 0 1; 0 1 1];

# standard algorithm
S = ipf(x,y)
1.33333  0.666667  2.0
2.66667  1.33333   4.0

# modified algorithm
M = ipfm(P,x,y)
4.0  0.0  0.00133422
0.0  2.0  5.99867


### Update

Due to the presence of Hadamard products in the iteration, any zeros in the initial $$A$$ matrix will be propagated into the solution.

So there's no need for separate algorithms, just pass the initial matrix as an input parameter:

function ipf(A,x,y)
u,v = x-x.+1, y-y.+1
for k = 1:2000
B = A .* (x*v') ./ ((A*v)*v')
A = B .* (u*y') ./ (u*(u'*B))
if norm(A-B) < 1e-6; break; end
end
return A
end

x,y = [4;8],  [4;2;6];

ipf( [1 1 1; 1 1 1], x,y)
1.33333  0.666667  2.0
2.66667  1.33333   4.0

ipf( [1 0 1; 0 1 1], x,y)
4.0  0.0  0.00133422
0.0  2.0  5.99867