The entries of the Sinkhorn limit of an $m \times n$ matrix
Eric Rowland Jason Wu


The entries of the Sinkhorn limit of an m × n matrix]The entries of the Sinkhorn 
 limit of an m × n matrix


	Department of Mathematics 

	Hofstra University 

	Hempstead, NY 

	USA



	Half Hollow Hills High School East 

	Dix Hills, NY 

	USA

Department of Mathematics 

	Cornell University 

	Ithaca, NY 

	USA



§ ABSTRACT 
We use Gröbner bases to compute the Sinkhorn limit of a positive 3 × 3 matrix A, showing that the entries are algebraic numbers with degree at most 6.
The polynomial equation satisfied by each entry is large, but we show that it has a natural representation in terms of linear combinations of products of minors of A.
We then use this representation to interpolate a polynomial equation satisfied by an entry of the Sinkhorn limit of a positive 4 × 4 matrix.
Finally, we use the PSLQ algorithm and 1.5 years of CPU time to formulate a conjecture, up to certain signs that we have not been able to identify, for a polynomial equation satisfied by an entry of the Sinkhorn limit of a positive m × n matrix.
In particular, we conjecture that the entries are algebraic numbers with degree at most m + n - 2m - 1.
This degree has a combinatorial interpretation as the number of minors of an (m - 1) × (n - 1) matrix, and the coefficients reflect new combinatorial structure on sets of minor specifications.


[
    Jason Wu
    September 4, 2024
=====================






§ INTRODUCTION


For a video introduction to this paper, see <https://youtu.be/-uIwboK4nwE>.

In a 1964 paper, Sinkhorn <cit.> considered the following iterative scaling process.
Let A be a square matrix with positive entries.
Scale the rows so that each row sum is 1.
Then scale the columns so that each column sum is 1; generically, this changes the row sums.
To restore the row sums 1, scale the rows again, then scale the columns, and so on.
Sinkhorn showed that the sequence of matrices obtained through this process converges to a matrix whose row and column sums are 1 (in other words, a doubly stochastic matrix).
We call this matrix the Sinkhorn limit of A and denote it (A).
Sinkhorn also showed that (A) is the unique doubly stochastic matrix S with the same size as A such that S = R A C for some diagonal matrices R and C with positive diagonal entries.
Here R can be taken to be the product of the row-scaling matrices and C the product of the column-scaling matrices.

The literature on questions related to the iterative scaling process is large;
Idel's extensive survey <cit.> covers results as of 2016.
Mathematical applications include preconditioning linear systems to improve numerical stability, approximating the permanent of a matrix, and determining whether a graph has a perfect matching.
In a number of other areas, iterative scaling was discovered independently multiple times <cit.>, and it is used in machine learning to efficiently compute optimal transport distances <cit.>.
As a result of its ubiquity and importance, many authors have been interested in fast algorithms for approximating Sinkhorn limits numerically <cit.>.

However, until recently, nothing was known about the exact values of the entries of Sinkhorn limits.
For a 2 × 2 matrix

    A =
    	[ a b; c d ]

with positive entries, Nathanson <cit.> showed that

    (A) =
    	1/√(a d) + √(b c)[ √(a d) √(b c); √(b c) √(a d) ].

In particular, the top left entry x of (A) satisfies

    (a d - b c) x^2 - 2 a d x + a d = 0.


For 3 × 3 matrices, analogous descriptions of the entries of (A) were only known in special cases.
In the case that A is a symmetric 3 × 3 matrix containing exactly 2 distinct entries, a formula for (A) was obtained by Nathanson <cit.>.
For a symmetric 3 × 3 matrix, Ekhad and Zeilberger <cit.> used Gröbner bases to compute, for each entry x of (A), a degree-4 polynomial of which x is a root.

For general 3 × 3 matrices, Chen and Varghese <cit.> used numeric experiments to conjecture that the entries of (A) generically have degree 6 over the field generated by the entries of A.


Applying the iterative scaling process to

    A =
    	[ 3 9 1; 3 2 9; 5 3 4 ]

gives the approximation

    (A) ≈[ 0.27667 0.64804 0.07527; 0.25194 0.13113 0.61692; 0.47138 0.22081 0.30780 ].

Let x be the top left entry of (A).
The PSLQ integer relation algorithm <cit.> can be used to recognize an algebraic number, given a sufficiently high-precision approximation.
For the approximation

    x ≈ 0.2766771162103280503525099931476512576251224460918253185145079454

with target degree 6, PSLQ produces the guess

    374752 x^6 - 220388 x^5 - 844359 x^4 - 125796 x^3 + 210897 x^2 + 14346 x - 12312 = 0.

This guess remains stable when we increase the precision or the target degree.


In Section <ref>, we prove the conjecture of Chen and Varghese by carrying out a Gröbner basis computation to obtain an explicit polynomial equation satisfied by an entry of (A) for a general 3 × 3 matrix A with positive entries.
In particular, we obtain formulas for the coefficients in Equation (<ref>).
We show that each coefficient in this polynomial can be written as a linear combination of products of minors of A (that is, determinants of square submatrices); this substantially reduces the amount of information required to specify the coefficients.

In Section <ref>, we use the form of the equation for 3 × 3 matrices to infer the form of the equation for general n × n matrices A.
In particular, we conjecture that the entries of (A) generically have degree 2 n - 2n - 1.
For a general 4 × 4 matrix, the Gröbner basis computation is infeasible, so instead we use PSLQ to recognize entries of Sinkhorn limits for enough integer matrices to solve for each coefficient in the equation.
We also obtain several coefficients in the equation for 5 × 5 matrices with this method.

In Section <ref>, we generalize to non-square matrices by defining the Sinkhorn limit of an m × n matrix A to be the matrix obtained by iteratively scaling so that each row sum is 1 and each column sum is m/n.
Again we use PSLQ and solve for coefficients in equations for matrices of various sizes.
We then interpolate formulas for these coefficients as functions of m and n.
By identifying combinatorial structure in these formulas, we build up to the following conjecture.
Let D(m, n) be the set of all pairs (R, C) where R ⊆{2, 3, …, m}, C ⊆{2, 3, …, n}, and R = C.
The set D(m, n) consists of the specifications of all minors of an m × n matrix A that do not involve the first row or first column.
For each S ⊆ D(m, n), let M(S) be the product of minors of A defined in Section <ref>.
Let _S, σ(S)(m, n) be the S×S matrix defined in Section <ref>;
this matrix resembles an adjacency matrix, and its entries are linear functions of m and n with signs determined by the sign alteration σ(S).


Let m ≥ 1 and n ≥ 1.
For each S ⊆ D(m, n), there exists a sign alteration σ(S) such that, for every m × n matrix A with positive entries, the top left entry x of (A) satisfies

    ∑_S ⊆ D(m, n)(_S, σ(S)(m, n)) M(S) x^S = 0.

In particular, x is algebraic over the field generated by the entries of A, with degree at most m + n - 2m - 1.


We conclude in Section <ref> with several open questions.
In particular, we mention that a more general process of iteratively scaling was introduced in 1937 by Kruithof in the context of predicting telephone traffic <cit.>.
We conjecture that the entries of this more general limit also have degree at most m + n - 2m - 1.

Our Mathematica package SinkhornPolynomials <cit.> uses the results of this paper to compute Sinkhorn limits rigorously for 3 × 3 matrices and conjecturally for 2 × n, 3 × 4, 3 × 5, and 4 × 4 matrices as well as their transposes.




§ THE SINKHORN LIMIT OF A 3 × 3 MATRIX


We refer to a matrix with positive entries as a positive matrix.
In this section, we determine (A) for a general positive 3 × 3 matrix

    A =
    	[ a_11 a_12 a_13; a_21 a_22 a_23; a_31 a_32 a_33 ],

first in Theorem <ref> as an explicit function of the entries of A and then in Theorem <ref> in a form that shows more structure.
We also introduce notation that will be used throughout the rest of the paper.

It suffices to describe the top left entry of (A).
This is because the iterative scaling process isn't sensitive to the order of the rows or the order of the columns.
Therefore, if R is the permutation matrix swapping rows 1 and i and C is the permutation matrix swapping columns 1 and j, then R (A) C = (R A C).
In particular, the (i, j) entry of (A) is equal to the (1, 1) entry of (R A C).
For example, the (2, 3) entry of (A) is equal to the (1, 1) entry of

    ([ a_23 a_22 a_21; a_13 a_12 a_11; a_33 a_32 a_31 ]).


To compute a polynomial equation satisfied by the top left entry of (A), we set up three matrices

    S =
    	[ s_11 s_12 s_13; s_21 s_22 s_23; s_31 s_32 s_33 ],
    	  
    	R =
    	[ r_1   0   0;   0 r_2   0;   0   0 r_3 ],
    	  
    	C =
    	[ c_1   0   0;   0 c_2   0;   0   0 c_3 ].

The matrix equation S = R A C gives the 9 equations

    s_11   = r_1 a_11 c_1		   	s_12   = r_1 a_12 c_2		   	s_13   = r_1 a_13 c_3 
    
    	s_21   = r_2 a_21 c_1		   	s_22   = r_2 a_22 c_2		   	s_23   = r_2 a_23 c_3 
    
    	s_31   = r_3 a_31 c_1		   	s_32   = r_3 a_32 c_2		   	s_33   = r_3 a_33 c_3,

and we obtain 6 equations from the requirement that S is doubly stochastic:

    s_11 + s_12 + s_13   = 1	   	s_11 + s_21 + s_31   = 1 
    
    	s_21 + s_22 + s_23   = 1	   	s_12 + s_22 + s_32   = 1 
    
    	s_31 + s_32 + s_33   = 1	   	s_13 + s_23 + s_33   = 1.

We would like to eliminate the 14 variables s_12, s_13, …, s_33, r_1, r_2, r_3, c_1, c_2, c_3 from this system of 15 polynomial equations, resulting in a single equation in the variables s_11, a_11, a_12, …, a_33.
In principle, this can be done by computing a suitable Gröbner basis.
In practice, the runtime is significantly affected by the algorithm used.
Mathematica's  function <cit.> with certain settings[Namely, .] computes a single polynomial in a couple seconds, whereas with other settings the computation does not finish after several hours.
The output gives the following result.


Let A be a positive 3 × 3 matrix.
The top left entry x of (A) satisfies b_6 x^6 + … + b_1 x + b_0 = 0, where the coefficients b_k appear in Table <ref> in factored form.




In particular, the degree of each entry of (A) is at most 6.
For the matrix A in Example <ref>, Theorem <ref> states that the top left entry x of (A) satisfies

    81 (374752 x^6 - 220388 x^5 - 844359 x^4 - 125796 x^3 + 210897 x^2 + 14346 x - 12312) = 0.

This agrees with Equation (<ref>).
This polynomial is irreducible, and this confirms the conjecture of Chen and Varghese <cit.> that the entries of (A) for positive 3 × 3 matrices A generically have degree 6.

Let f(x) = b_6 x^6 + b_5 x^5 + … + b_1 x + b_0 be the polynomial in Theorem <ref>.
Project A to a symmetric matrix by setting a_21 = a_12, a_31 = a_13, and a_32 = a_23.
Then the projection of f(x) factors as -((a_11 a_23 - a_12 a_13) x - a_11 a_23)^2 g(x), where g(x) is the degree-4 polynomial computed by Ekhad and Zeilberger <cit.> for the top left entry.

An obvious question is whether there is a better way to write the polynomial f(x) in Theorem <ref>.
In fact there is, using determinants.
We first observe that b_k contains the factor a_11^5 - k for each k ∈{0, 1, …, 5}.
This suggests that the scaled polynomial a_11 f(x) is more natural than f(x), since the coefficient of x^k in a_11 f(x) contains the factor a_11^6 - k not just for k ∈{0, 1, …, 5} but for all k ∈{0, 1, …, 6} (where in fact the coefficient of x^6 also contains a_11^1).

Now the leading coefficient a_11 b_6 of a_11 f(x) is the product of 6 minors of A.
More specifically, it is the product of all the minors of A that involve a_11.
Furthermore, the constant coefficient a_11 b_0 of a_11 f(x) is the product of a_11^6 and the 6 minors of A that do not involve the first row or first column (including the determinant 1 of the 0 × 0 matrix).

To rewrite the other coefficients, we would like to interpolate between these products for a_11 b_6 and a_11 b_0.
One possibility is that a_11 b_k is a linear combination of products of minors, where each product consists of


  * k minors that involve the first row and first column,

  * 6 - k minors that do not involve the first row or first column, and

  * a_11^6 - k (equivalently, one factor a_11 for each of the latter).

We have seen that this is the case for a_11 b_6 and a_11 b_0.
Remarkably, it turns out that a_11 b_5, a_11 b_4, …, a_11 b_1 can be written this way as well.
Some notation for these minors will be useful.


Define

    D(m, n) = {(R, C) : R ⊆{2, 3, …, m} and C ⊆{2, 3, …, n} and R = C}.

Define A_R, C to be the submatrix of A obtained by extracting the rows indexed by R and the columns indexed by C.
For each (R, C) ∈ D(m, n), define

    Δ([ R; C ])
    	   =  A_{1}∪ R, {1}∪ C
    Γ([ R; C ])
    	   = a_11 A_R, C.

The minor
Δ([ R; C ])
involves a_11, and
Γ([ R; C ])
is the product of a_11 and a minor that does not involve the first row or first column.
This notation does not reflect the dependence on A, but the matrix will be clear from context.
Each subset S ⊆ D(m, n) specifies a monomial in the expressions
Δ([ R; C ])
and
Γ([ R; C ]), namely

    M(S)
    	= ∏_(R, C) ∈ SΔ([ R; C ])
    	·∏_(R, C) ∈ D(m, n) ∖ SΓ([ R; C ]).

Each element of D(m, n) contributes to the monomial M(S) as the argument of either Δ or Γ; indeed we specify each monomial by the elements that appear as arguments of Δ.
To make a subset S = {(R_1, C_1), (R_2, C_2), …, (R_k, C_k)} easier to read, we format it as

    S =
    	
    		R_1     R_2    ⋯    R_k
    
    		C_1     C_2    ⋯    C_k.




For m = 3 and n = 3, we have
D(3, 3) =
		{}   {2}   {2}   {3}   {3}   {2, 3}
{}   {2}   {3}   {2}   {3}   {2, 3}.
Moreover, from the expressions for b_6 and b_0 in Table <ref>, we have

    a_11 b_6    =
    		Δ([ {}; {} ])
    		Δ([ {2}; {2} ])
    		Δ([ {2}; {3} ])
    		Δ([ {3}; {2} ])
    		Δ([ {3}; {3} ])
    		Δ([ {2, 3}; {2, 3} ])
    		=
    		M({}   {2}   {2}   {3}   {3}   {2, 3}
    {}   {2}   {3}   {2}   {3}   {2, 3}) 
    
    	a_11 b_0    =
    		Γ([ {}; {} ])
    		Γ([ {2}; {2} ])
    		Γ([ {2}; {3} ])
    		Γ([ {3}; {2} ])
    		Γ([ {3}; {3} ])
    		Γ([ {2, 3}; {2, 3} ])
    		=
    		M(
    ).

Each of these monomials involves Δ or Γ but not both, whereas the monomial

    M({}   {3}
    {}   {2})
    	=
    		Δ([ {}; {} ])
    		Γ([ {2}; {2} ])
    		Γ([ {2}; {3} ])
    		Δ([ {3}; {2} ])
    		Γ([ {3}; {3} ])
    		Γ([ {2, 3}; {2, 3} ])

involves both Δ and Γ, for example.


We continue to let A be the general 3 × 3 matrix in Equation (<ref>), and let x be the top left entry of (A).
To rewrite the coefficients a_11 b_k, we look for coefficients c_S ∈, indexed by subsets S ⊆ D(3, 3), such that

    a_11 b_k = ∑_S ⊆ D(3, 3) 
    S = k c_S M(S)

for each k ∈{0, 1, …, 6}.
This will give the polynomial equation

    ∑_k = 0^6(∑_S ⊆ D(3, 3) 
    S = k c_S M(S)) x^k = 0.

The expressions for a_11 b_6 and a_11 b_0 in Example <ref> allow us to choose c_D(3, 3) = 1 and c_{} = 1.
For each k ∈{1, 2, 4, 5}, one checks that a_11 b_k can be written uniquely as a linear combination of the monomials M(S) where S = k.
This determines c_S for such subsets S.
For example,

    a_11 b_1 =
    		-3 M({}
    {})
    		- M({2}
    {2})
    		- M({2}
    {3})
    		- M({3}
    {2})
    		- M({3}
    {3})
    		+ M({2, 3}
    {2, 3}).


However, the coefficient a_11 b_3 has multiple representations.
The family of such representations is 1-dimensional, due to the relation

    M({}   {2}   {3}
    {}   {2}   {3})
    	+ M({}   {2}   {2, 3}
    {}   {3}   {2, 3})
    	+ M({}   {3}   {2, 3}
    {}   {2}   {2, 3}) 
    
    	+ M({2}   {2}   {3}
    {2}   {3}   {2})
    	+ M({2}   {3}   {2, 3}
    {2}   {3}   {2, 3})
    	+ M({2}   {3}   {3}
    {3}   {2}   {3}) 
    
    	=
    	M({}   {2}   {2, 3}
    {}   {2}   {2, 3})
    	+ M({}   {2}   {3}
    {}   {3}   {2})
    	+ M({}   {3}   {2, 3}
    {}   {3}   {2, 3}) 
    
    	+ M({2}   {2}   {3}
    {2}   {3}   {3})
    	+ M({2}   {3}   {3}
    {2}   {2}   {3})
    	+ M({2}   {3}   {2, 3}
    {3}   {2}   {2, 3}).

Therefore the expression for b_3 in Table <ref> does not uniquely determine the coefficients c_S where S = 3.
However, there is additional information we can use to obtain uniqueness.
The polynomial f(x) possesses symmetries arising from the following invariance properties.

Since the iterative scaling process isn't sensitive to row order, the top left entry of (A) is invariant under row permutations of A that fix the first row.
Similarly for column permutations.
Additionally, the top left entry of (A) is invariant under transposition of A; this follows from Sinkhorn's result that (A) is the unique doubly stochastic matrix S such that S = R A C for some diagonal matrices R and C.
This suggests the following equivalence relation.


Let S and T be subsets of D(m, n).
We write T ≡ S if the set of minors specified by T is transformed into the set of minors specified by S by some composition of


  * row permutations that fix the first row,

  * column permutations that fix the first column, and

  * transposition if m = n.

For each S ⊆ D(m, n), define the class sum

    Σ(S) = ∑_T ⊆ D(m, n) 
     T ≡ S M(T)

to be the sum of the monomials corresponding to the elements in the equivalence class of S.



Let m = 3 and n = 3, and consider the subset
S =
	{2}
{2} of size 1.
The equivalence class of S is

    {{2}
    {2},  {2}
    {3},  {3}
    {2},  {3}
    {3}}

since these 4 specifications of 1 × 1 minors can be obtained from each other by row and column permutations.
This equivalence class is reflected in the linear combination (<ref>) for a_11 b_1.
Namely, the four monomials

    M({2}
    {2}),  
    	M({2}
    {3}),  
    	M({3}
    {2}),  
    	M({3}
    {3})

all have the same coefficient -1.
In particular, their contribution to a_11 b_1 is

    -Σ({2}
    {2})
    	=
    	-M({2}
    {2})
    	-
    	M({2}
    {3})
    	-
    	M({3}
    {2})
    	-
    	M({3}
    {3}).



Motivated by Example <ref>, we add the constraint that c_T = c_S when T ≡ S, so that the coefficients c_S reflect the symmetries of f(x).
With this constraint, the coefficient a_11 b_3 has a unique representation as a linear combination of monomials M(S) where S = 3.
Writing each coefficient in Table <ref> using class sums gives the following improvement of Theorem <ref>.
Here d_k = a_11 b_k.


Let A be a positive 3 × 3 matrix.
The top left entry x of (A) satisfies d_6 x^6 + d_5 x^5 + d_4 x^4 + d_3 x^3 + d_2 x^2 + d_1 x + d_0 = 0, where

    d_6    =
    		Σ({}   {2}   {2}   {3}   {3}   {2, 3}
    {}   {2}   {3}   {2}   {3}   {2, 3}) 
    
    	d_5    =
    		-3 Σ({}   {2}   {2}   {3}   {3}
    {}   {2}   {3}   {2}   {3})
    		- Σ({}   {2}   {2}   {3}   {2, 3}
    {}   {2}   {3}   {2}   {2, 3})
    		+ Σ({2}   {2}   {3}   {3}   {2, 3}
    {2}   {3}   {2}   {3}   {2, 3}) 
    
    	d_4    =
    		4 Σ({}   {2}   {2}   {3}
    {}   {2}   {3}   {2})
    		+ Σ({}   {2}   {3}   {2, 3}
    {}   {2}   {3}   {2, 3})
    		- 3 Σ({2}   {2}   {3}   {3}
    {2}   {3}   {2}   {3}) 
    
    	d_3    =
    		-4 Σ({}   {2}   {2}
    {}   {2}   {3})
    		- 5 Σ({}   {2}   {3}
    {}   {2}   {3})
    		+ Σ({}   {2}   {2, 3}
    {}   {2}   {2, 3})
    		+ Σ({2}   {2}   {3}
    {2}   {3}   {2})
    		- Σ({2}   {3}   {2, 3}
    {2}   {3}   {2, 3}) 
    
    	d_2    =
    		4 Σ({}   {2}
    {}   {2})
    		- 3 Σ({}   {2, 3}
    {}   {2, 3})
    		+ Σ({2}   {3}
    {2}   {3}) 
    
    	d_1    =
    		-3 Σ({}
    {})
    		- Σ({2}
    {2})
    		+ Σ({2, 3}
    {2, 3}) 
    
    	d_0    =
    		Σ(
    ).



Several class sums do not appear in Theorem <ref>, namely

    Σ({}   {2}   {2}   {2, 3}
    {}   {2}   {3}   {2, 3}),  Σ({2}   {2}   {3}   {2, 3}
    {2}   {3}   {2}   {2, 3}),  Σ({2}   {2}   {2, 3}
    {2}   {3}   {2, 3}),  Σ({2}   {2}
    {2}   {3}),  Σ({2}   {2, 3}
    {2}   {2, 3}).

These class sums get assigned the coefficient 0 when the coefficients d_4, d_3, d_2 are written as linear combinations of Σ(S).
In total, there are 24 equivalence classes of subsets S ⊆ D(3, 3), so we have compressed the information in Table <ref> down to a function from the set of these 24 equivalence classes to the set {-5, -4, -3, -1, 0, 1, 4}.

For the particular matrix A in Example <ref>, Theorem <ref> gives 91064736 x^6 - 53554284 x^5 - 205179237 x^4 - 30568428 x^3 + 51247971 x^2 + 3486078 x - 2991816 = 0 for the top left entry.
This equation can be computed quickly with SinkhornPolynomials <cit.>.
Dividing this equation by 243 produces Equation (<ref>).

Some special cases can be obtained from Theorem <ref>.
For example, the following corollary shows that the degree drops if a minor is 0.
(If multiple minors are 0, the degree can drop further.)


Let A be a positive 3 × 3 matrix, and let a_11 be the (1, 1) entry of A.
If one of the 2 × 2 minors involving a_11 is 0 and all minors not involving a_11 are not 0, then the top left entry of (A) has degree at most 5.



The coefficient of x^6 in Theorem <ref> is

    d_6
    	=
    	Σ({}   {2}   {2}   {3}   {3}   {2, 3}
    {}   {2}   {3}   {2}   {3}   {2, 3})
    	=
    	M({}   {2}   {2}   {3}   {3}   {2, 3}
    {}   {2}   {3}   {2}   {3}   {2, 3})
    	=
    	∏_(R, C) ∈ D(3, 3)Δ([ R; C ]).

Since one of the minors involving a_11 is 0, this product is 0, so d_6 = 0.
On the other hand, since all minors not involving a_11 are not 0, we have

    d_0
    	=
    	Σ(
    )
    	=
    	M(
    )
    	=
    	∏_(R, C) ∈ D(3, 3)Γ([ R; C ])
    	≠ 0.

Therefore the polynomial in Theorem <ref> is not the 0 polynomial, and its degree in x is at most 5.


If a matrix is sufficiently degenerate, then Theorem <ref> is vacuously true and does not immediately give any information about (A).
For example, let

    A =
    	[ 4 5 6; 1 2 3; 2 4 6 ].

Here rows 2 and 3 are scalar multiples of each other, so
Δ([ {2, 3}; {2, 3} ])
	=  A
	= 0
and
Γ([ {2, 3}; {2, 3} ])
	= a_11 A_{2, 3}, {2, 3}
	= 0.
Since each monomial M(S) contains one of these two determinants as a factor, we have M(S) = 0 for all S ⊆ D(3, 3).
Therefore d_k = 0 for each k.
However, we can still use Theorem <ref> to determine the entries of (A), as the following result shows.


Let A be a positive 3 × 3 matrix, and let a_ij be the (i, j) entry of A.
If rows 2 and 3 are scalar multiples of each other, then the top left entry x of (A) satisfies e_3 x^3 + e_2 x^2 + e_1 x + e_0 = 0, where

    e_3    = a_11 (a_11 a_22 - a_12 a_21) (a_11 a_23 - a_13 a_21) 
    
    	e_2    = a_11 (a_11 a_12 a_21 a_23 + a_11 a_13 a_21 a_22 - 3 a_11^2 a_22 a_23 + a_12 a_13 a_21^2) 
    
    	e_1    = 3 a_11^3 a_22 a_23
    
    	e_0    = -a_11^3 a_22 a_23.

In particular, x is independent of row 3.


The analogous result holds if columns 2 and 3 are scalar multiples of each other.
For the matrix in Equation (<ref>), Corollary <ref> gives 24 (3 x^3 - 25 x^2 + 48 x - 16) = 0.


The idea is to carefully factor out the minors that are 0.
First we describe how to obtain the polynomial e_3 x^3 + e_2 x^2 + e_1 x + e_0, and then we justify it.

Begin with a general 3 × 3 matrix A with symbolic entries;
in particular, there are no algebraic relations between the entries.
Let r, s be symbols; eventually we will set r to be the scalar factor a_31/a_21.
Apply Theorem <ref> to A to obtain a polynomial equation satisfied by the top left entry of (A).
Then replace each instance of A in this polynomial with s  A_{2, 3}, {2, 3}.
By the definition of M(S), each monomial now contains A_{2, 3}, {2, 3} as factor; divide by this factor.
Finally, replace a_3,j with r a_2,j for each j ∈{1, 2, 3}.
The resulting polynomial factors as a product two cubic polynomials in x.
One of these cubic factors is independent of r and s;
this is the polynomial e_3 x^3 + e_2 x^2 + e_1 x + e_0.

To justify the construction, let A now be the matrix in the statement of the corollary, let r = a_31/a_21 be the scalar factor, and fix a real number s ≠ 0.
We approximate A by matrices with no 0 minors.
Namely, let B(t) be a continuous (3 × 3 matrix)-valued function such that lim_t → 0^+ B(t) = A and, for all t > 0, no minor of B(t) is 0.
Further, we assume for all t > 0 that B(t)/ B(t)_{2, 3}, {2, 3} = s.
Apply Theorem <ref> to B(t) where t > 0.
Since B(t) = s  B(t)_{2, 3}, {2, 3} by assumption, the polynomial given by Theorem <ref> is of the form ( B(t)_{2, 3}, {2, 3}) g(t, x).
The entries of (B(t)) are continuous functions of the entries of B(t), and the roots of a polynomial are continuous functions of its coefficients, so g(x) lim_t → 0^+ g(t, x) is a polynomial for the top left entry of (A).
Moreover, since only one of the two cubic factors of g(x) is independent of s, the top left entry of (A) is a root of that cubic factor.





§ SQUARE MATRICES AND SOLVING FOR COEFFICIENTS


In this section, we use Theorem <ref> to infer the form of an equation satisfied by the top left entry of the Sinkhorn limit of a positive n × n matrix.
We then interpolate the coefficients in the equation for 4 × 4 matrices to obtain Conjecture <ref>.

For a positive 2 × 2 matrix A, the top left entry x of (A) satisfies Equation (<ref>), namely (a_11 a_22 - a_12 a_21) x^2 - 2 a_11 a_22  x + a_11 a_22 = 0.
Multiplying by a_11, this can be written using Δ and Γ as

    Δ([ {}; {} ])
    	Δ([ {2}; {2} ])
    	x^2
    	- 2
    	Δ([ {}; {} ])
    	Γ([ {2}; {2} ])
    	x
    	+
    	Γ([ {}; {} ])
    	Γ([ {2}; {2} ])
    	= 0

or, equivalently,

    M({}   {2}
    {}   {2})
    	x^2
    	- 2
    	M({}
    {})
    	x
    	+
    	M(
    )
    	= 0.

Along with Theorem <ref>, this suggests that, for an n × n matrix, the coefficient of x^k is a -linear combination of the monomials M(S) where S = k.
In particular, we expect the degree of the equation to be D(n, n) = ∑_k = 0^n - 1n - 1k^2 = 2 n - 2n - 1, so that generically the entries of (A) for a 4 × 4 matrix have degree 63 = 20 and for a 5 × 5 matrix have degree 84 = 70.
We generalize the notation for the coefficients c_S from the previous section to c_S(n) for an n × n matrix, where c_S(3)  c_S.


Let n ≥ 1.
There exist integers c_S(n), indexed by subsets S ⊆ D(n, n), such that, for every positive n × n matrix A, the top left entry x of (A) satisfies

    ∑_k = 0^2 n - 2n - 1(∑_S ⊆ D(n, n) 
    S = k c_S(n) M(S)) x^k = 0

or, equivalently,

    ∑_S ⊆ D(n, n) c_S(n) M(S) x^S = 0.



This conjecture predicts that, since both
Δ([ R; C ]) and Γ([ R; C ])
are homogeneous degree-(R + 1) polynomials in the entries of A, each monomial M(S) (and therefore also the coefficient of each x^k in Equation (<ref>)) is a homogeneous polynomial with degree

    ∑_(R, C) ∈ D(n, n) (R + 1)
    	= ∑_i = 0^n - 1n - 1i^2 (i + 1)
    	= n + 1/22 n - 2n - 1.

For n = 1, 2, 3, …, this degree is 1, 3, 12, 50, 210, 882, 3696, 15444, … A092443.

We mention two surprising properties of the polynomial in Theorem <ref> that we expect to generalize to the polynomial in Conjecture <ref>.
The first is a symmetry.
Each coefficient d_k in Theorem <ref> can be obtained from d_6 - k by replacing

    Δ([ R; C ])
    	↦Γ([ {2, 3}∖ R; {2, 3}∖ C ])
    	  and  Γ([ R; C ])
    	↦Δ([ {2, 3}∖ R; {2, 3}∖ C ]).

The second is that, for each k, the sum of the coefficients c_S(3) for S = k is a signed binomial coefficient, namely

    ∑_S ⊆ D(3, 3) c_S(3) x^S = (x - 1)^6.

Analogous statements also hold for the quadratic polynomial in Equation (<ref>);
each coefficient d_k is related to d_2 - k by a symmetry, and

    ∑_S ⊆ D(2, 2) c_S(2) x^S = (x - 1)^2.

We do not have explanations for either of these properties.

It remains to determine the coefficients c_S(n).
They are not uniquely determined by the conditions in Conjecture <ref> since we can scale the polynomial.
To remove this source of non-uniqueness, we define c_{}(n) = 1 for all n ≥ 1.
Based on the polynomials for n = 2 and n = 3, we conjecture that c_D(n, n)(n) = 1 for all n ≥ 2.
For n = 3 we determined the coefficients c_S(3) from the output of a Gröbner basis computation, but for n ≥ 4 this computation seems to be infeasible.
For a general 4 × 4 matrix, we aborted the computation after 1 week.

Instead, we generate many pseudorandom n × n matrices A, identify the top left entry of (A) as an algebraic number for each, and set up systems of linear equations in c_S(n).
We describe these three steps next.

For the first step, we generate matrices with entries from {1, 2, …, 20}.
Since the entries are integers, each M(S) is also an integer; this will be important in the third step.
For each matrix, we check that none of its minors are 0, since a 0 minor implies M(S) = 0 for several S, removing the dependence on the corresponding coefficients c_S(n).
If any minors are 0, we discard that matrix.

In the second step, for each matrix A generated in the first step, we determine a polynomial equation satisfied by the top left entry x of (A).
There are two possible methods.
One method is to use Gröbner bases;
this is faster than the Gröbner basis computation for a matrix with symbolic entries, but for n ≥ 5 it is still slow.
Therefore we use another method, which is to apply the iterative scaling process to obtain a numeric approximation to (A) and then guess a polynomial for its top left entry.
We begin by numericizing the integer entries of A to high precision.
For n = 4 we use precision 2^12, and for n = 5 we use precision 2^15.
Then we iteratively scale until we reach a fixed point.
The precision of the entries drops during the scaling process, but with the initial precisions 2^12 and 2^15 we get entries with sufficiently high precision that we can reliably recognize them using PSLQ.
The expected degree of the polynomial is d = 2 n - 2n - 1 according to Conjecture <ref>.
Building in redundancy, we use Mathematica's  <cit.> to approximate x by an algebraic number with target degree d + 2.
When the output has degree d, which is almost always the case, this is strong evidence that the approximation is in fact the exact algebraic number we seek.
Occasionally the output has degree less than d, in which case we discard the matrix; for example, one 4 × 4 matrix produced a polynomial with degree 8 rather than 20, presumably because the general degree-20 polynomial, when evaluated at the entries of the matrix, is reducible.
Finally, we perform a check on the output by computing the ratio of its leading coefficient to its constant coefficient.
This ratio should be M(D(n, n))/M({}), assuming Conjecture <ref> is correct and c_D(n, n)(n) = 1.
All outputs passed this test.
We record the matrix A along with the polynomial equation satisfied by x, and this will give us 1 equation in the third step.
(In fact we can get n^2 equations by applying PSLQ to each entry of the numeric approximation to (A) and recording each polynomial along with the matrix obtained by swapping the appropriate rows and columns of A.
For n = 5 this is worthwhile, since our implementation of the iterative scaling process takes roughly 3 minutes to reach a fixed point.)

The third step is to determine the coefficients c_S(n) in the coefficient of x^k in Conjecture <ref>.
For a given k, we do this by solving a system of linear equations involving the coefficients c_S(n) where S = k.
We assume c_T(n) = c_S(n) if T ≡ S, so it suffices to determine c_S(n) for one representative S of each equivalence class, analogous to Theorem <ref>.
This reduces the number of unknown coefficients, which reduces the number of equations we need, which reduces the number of matrices A we apply the iterative scaling process to in the second step above.
However, first we must partition {S ⊆ D(n, n) : S = k} into its equivalence classes under ≡.
Some care must be taken to do this efficiently;
we make use of the fact that row permutations commute with column permutations, so it suffices to apply row permutations first.
Once we have computed the equivalence classes, we take each polynomial computed in the second step above, scale it by an integer so that its leading coefficient is M(D(n, n)) (and therefore its constant coefficient is M({})), extract the coefficient of x^k, and set this coefficient equal to ∑_S c_S(n) Σ(S) where the sum is over one representative from each equivalence class of size-k subsets.
The number of equivalence classes tells us how many such equations we need in order to solve for the unknown coefficients c_S(n).
We include more equations than necessary in the system, building in redundancy, so that if a solution is found then we can be confident that the conjectured form is correct.
Then we solve the system.
In practice, even setting up the system can be computationally expensive.
For n = 5 and k = 4, our initial implementation took 6 days to set up the system (before solving!) because there are 1518 unknown coefficients c_S(n), so we need at least that many equations, and each equation involves 704 = 916895 monomials, each of which is a product of 70 determinants.
For small k, a more efficient way to construct each equation is to take advantage of the fact that, for each pair S, T of subsets of D(5, 5), the products M(S) and M(T) have most factors in common, and almost all are Γ factors.
Therefore, we form an equivalent but substantially simpler equation by dividing both sides by
M({})
	= ∏_(R, C) ∈ D(5, 5)Γ([ R; C ]).
On the left, we compute M({}) once and divide the extracted coefficient of x^k by M({}).
On the right, we divide each M(S) by M({}).
Instead of computing M(S)/M({}) from definitions, we precompute the ratio
Δ([ R; C ])
	/
	Γ([ R; C ])
for each (R, C) ∈ D(5, 5); then, for each S with S = k, we use the precomputed ratios to compute

    M(S)/M({})
    	=
    	∏_(R, C) ∈ SΔ([ R; C ])/Γ([ R; C ]).

Here we're using the fact that no minor is 0 to divide by
Γ([ R; C ]).
This method is also used by SinkhornPolynomials to compute polynomials more quickly when no minor is 0.

We now carry out these three steps for n = 4 to obtain a polynomial for the top left entry x of (A) for 4 × 4 matrices A.
This polynomial is analogous to the polynomial in Theorem <ref> for 3 × 3 matrices.
The number of equivalence classes of size-k subsets of D(4, 4) for k = 0, 1, …, 20 is

    1, 4, 12, 40, 123, 324, 724, 1352, 2108, 2760, 3024, 2760, …, 4, 1.

The number of unknown coefficients c_S(4) is the sum of these numbers, which is 17920.
We use the definition c_{}(4) = 1 and the conjecture c_D(4, 4)(4) = 1, and we solve the remaining 19 systems of linear equations, the largest of which requires 3024 equations.
Unfortunately, for each k ∈{4, 5, …, 16} the system has multiple solutions.
For example, when we solve the system for k = 4, only 104 of the 123 unknown coefficients c_S(4) with S = 4 are uniquely determined;
the other 19 are parameterized by 2 free variables.
For k = 10, the solution space has dimension 1141.
This is not a weakness of the interpolation strategy but rather implies that the coefficients of x^4, x^5, …, x^16 have multiple representations as linear combinations of class sums Σ(S) and therefore that relations exist among these class sums.
However, by setting the free variables to 0 (or any other values), we obtain variable-free coefficients c_S(4).


The 17920 coefficients

    c_D(4, 4)(4)    = 1 
         ⋮
    
    	c_{}
    {}(4)    = -4 
    
    	c_{2}
    {2}(4)    = -2 
    
    	c_{2, 3}
    {2, 3}(4)    = 0 
    
    	c_{2, 3, 4}
    {2, 3, 4}(4)    = 2 
    
    	c_{}(4)    = 1

determine a polynomial equation satisfied by the top left entry x of the Sinkhorn limit of a 4 × 4 matrix, namely

    ∑_k = 0^20(∑ c_S(4) Σ(S)) x^k = 0

where the inner sum is over one representative S from each equivalence class of size-k subsets of D(4, 4).


The full list of coefficients is included in SinkhornPolynomials <cit.>.


Consider the matrix

    A =
    	[ 3 1 2 2; 2 2 2 1; 1 2 3 2; 1 4 2 3 ],

and let x be the top left entry of (A).
Conjecture <ref> gives

    382625520076800 x^20 - 15753370260418560 x^19 + 224644720812019200 x^18
    
    	- 1949693785825830912 x^17 + 11625683820163305984 x^16 - 50547801347982259200 x^15
    
    	+ 165827284134596798976 x^14 - 419342005165888558080 x^13 + 828111699533723747328 x^12
    
    	- 1284220190788992755712 x^11 + 1558933050581256001536 x^10 - 1456458194243244008448 x^9 
    
    	+ 999710159534823121920 x^8 - 435645828109071673344 x^7 + 31060141423020794880 x^6 
    
    	+ 122853118332060905472 x^5 - 110123924197151416320 x^4 + 53612068706701295616 x^3 
    
    	- 16383341182381572096 x^2 + 2975198930601246720 x - 246790694704250880
    	= 0.



The values of the coefficients c_S(4) in Conjecture <ref> are not all canonical, since we don't have natural conditions under which they are uniquely determined.
However, we continue under the assumption that there is a unique natural function c_S(n) and seek to identify it.

In principle, we can use the same method to interpolate a polynomial equation for n × n matrices for any given n.
However, the computation is formidable.
For 5 × 5 matrices, we were able to compute the coefficients for k = 4 and k = 5 (where the families of representations are respectively 8-dimensional and 44-dimensional), but for k = 6 we could not get the computation of equivalence classes to finish (aborting after 2 weeks).
Completing the computation up to the halfway point k = 1/284 = 35, after which we could use the conjectural symmetry, seems to be infeasible.

For 6 × 6 matrices, PSLQ must recognize algebraic numbers with degree 105 = 252.
We couldn't get any of these computations to finish with sufficiently high precision.




§ RECTANGULAR MATRICES AND COMBINATORIAL STRUCTURE


The coefficients c_S(n) for n ∈{2, 3, 4, 5} that we have computed so far are not enough to guess the general formula for c_S(n).
In this section, we expand the scope to include matrices that are not necessarily square.
This allows us to compute enough coefficients to identify formulas in special cases.
These formulas reveal the relevant combinatorial structure on subsets S, and this allows us to piece together all but one detail of the general formula for the coefficients, resulting in Conjecture <ref>.

Doubly stochastic matrices are necessarily square.
This is because, in every matrix, the sum of the row sums is equal to the sum of the column sums.
Therefore we must generalize how we scale.


Let A be a positive m × n matrix.
The Sinkhorn limit of A is the matrix obtained by iteratively scaling so that each row sum is 1 and each column sum is m/n.
Its existence was established (in a more general form) by Sinkhorn in a 1967 paper <cit.>.


Conjecture <ref> generalizes to m × n matrices as follows.
We extend the coefficients c_S(n) in the previous section to coefficients c_S(m, n), where c_S(n, n) = c_S(n).
For each m and n, we scale the coefficients so that c_{}(m, n) = 1.
We have D(m, n) = ∑_k = 0^min(m, n) - 1m - 1kn - 1k = m + n - 2m - 1.


Let m ≥ 1 and n ≥ 1.
There exist rational numbers c_S(m, n), indexed by subsets S ⊆ D(m, n), such that, for every positive m × n matrix A, the top left entry x of (A) satisfies

    ∑_S ⊆ D(m, n) c_S(m, n) M(S) x^S = 0.

In particular, x has degree at most m + n - 2m - 1.


To interpolate values of c_S(m, n), we use the method described in Section <ref>.
Overall, we used 1.5 years of CPU time to iteratively scale matrices of various sizes and recognize 102000 algebraic numbers.
An additional month of CPU time was spent setting up and solving systems of linear equations in the coefficients c_S(m, n), which resulted in the identification of 63000 rational coefficients (and an additional 56000 coefficients parameterized by free variables).

Rather than fixing m and n and varying S as in previous sections, we change our perspective now by fixing S and working to identify c_S(m, n) as a function of m and n.
Accordingly, we define c_S(m, n) to be 0 if S ⊈ D(m, n) (that is, S contains row or column indices that are larger than an m × n matrix supports).
We begin with subsets S where S = 1.


Let S = {}
{}.
Equation (<ref>) implies c_S(2, 2) = -2, Theorem <ref> implies c_S(3, 3) = -3, and Conjecture <ref> implies c_S(4, 4) = -4.
The values of c_S(m, n) we computed for several additional matrix sizes appear in the following table.

    [       n = 1     2     3     4; m = 1    -1    -2    -3    -4;     2    -1    -2    -3    -4;     3    -1    -2    -3    -4;     4    -1    -2    -3    -4 ]

This suggests the formula c_S(m, n) = -n.


There is an asymmetry in our definition of (A) for non-square matrices, since (A) has row sums 1 and column sums m/n.
However, we expect some symmetry in the coefficients c_S(m, n) since m (A) = n (A).
We can obtain this symmetry by considering 1/m(A), which has row sums 1/m and column sums 1/n.
The form of the polynomial for the top left entry y of 1/m(A) can be obtained from Conjecture <ref> by substituting x = m y.
Therefore, the coefficients c_S(m, n) should satisfy m^S c_S(m, n) = n^S c_S(n, m), where S is defined by

    (
    		R_1     R_2    ⋯    R_k
    
    		C_1     C_2    ⋯    C_k)
    	=
    	
    		C_1     C_2    ⋯    C_k
    
    		R_1     R_2    ⋯    R_k.



As in the previous example, let S = {}
{}.
Since S = S, we have m c_S(m, n) = n c_S(n, m).
In other words, m c_S(m, n) is a symmetric function of m and n.
Indeed, the table of values suggests m c_S(m, n) = -m n.


The coefficients for other size-1 subsets S also seem to be simple polynomial functions of m and n.


For S = {2}
{2}, the data suggests m c_S(m, n) = m + n - m n (for all m ≥ 2 and n ≥ 2).
For S = {2}
{3}, the data suggests the same formula m c_S(m, n) = m + n - m n); indeed this is expected, since {2}
{2}≡{2}
{3}.
Since each subset
R 

		C with R = i is equivalent to S = {2, 3, …, i + 1}
{2, 3, …, i + 1}, it suffices to consider the latter.
For S = {2, 3}
{2, 3}, the data suggests m c_S(m, n) = 2 m + 2 n - m n.
For S = {2, 3, 4}
{2, 3, 4}, it suggests m c_S(m, n) = 3 m + 3 n - m n.


These formulas lead to the following.


Let m ≥ 1 and n ≥ 1, and let
S =
	
		R 

		C ⊆ D(m, n).
Then m c_S(m, n) = R (m + n) - m n.


Next we consider subsets S where S = 2.


Let S = {}   {2}
{}   {2}.
Several values of m^2 c_S(m, n) appear in the following table.

    [       n = 2     3     4     5; m = 2     4    12    24    40;     3    12    36    72   120;     4    24    72   144   240;     5    40   120   240   400 ]

This suggests that m^2 c_S(m, n) = (m - m n) (n - m n).



Let S = {2}   {2}
{2}   {3}.
The two minor specifications in S involve the same row but different columns.
Therefore there is no subset T in the equivalence class of S such that T = T, so we do not expect m^2 c_S(m, n) to be a symmetric function of m and n.
Here are several of its values (where the bottom right entry is not included because we have no data for 6 × 6 matrices):

    [       n = 3     4     5     6; m = 2    -3     0     5    12;     3     0    16    40    72;     4     9    48   105   180;     5    24    96   200   336;     6    45   160   325 ]

This suggests that m^2 c_S(m, n) = (n - m n) (2 m + n - m n).


By interpolating polynomial formulas for additional subsets S, one conjectures that the scaled coefficient m^S c_S(m, n) is a polynomial function of m and n with degree S in each variable.
Moreover, there seem to be two basic relationships between minor specifications in S that play a central role in the form of the polynomial function.
These are illustrated by the previous two examples.


We say that two minor specifications (R_1, C_1), (R_2, C_2) ∈ D(m, n) are linked if either of the following conditions holds.


  * Their sizes differ by 1, and the smaller is a subset of the larger.
That is,


  * R_1 = R_2 ∖{i} and C_1 = C_2 ∖{j} for some i ∈ R_2 and j ∈ C_2, or

  * R_2 = R_1 ∖{i} and C_2 = C_1 ∖{j} for some i ∈ R_1 and j ∈ C_1.

In this case, we say that they form a type-1 link.

  * Their sizes are the same, and they differ in exactly 1 row index or 1 column index.
That is,


  * R_1 = R_2 and C_1 ∖{i} = C_2 ∖{j} for some i ∈ C_1 and j ∈ C_2 such that i ≠ j, or

  * C_1 = C_2 and R_1 ∖{i} = R_2 ∖{j} for some i ∈ R_1 and j ∈ R_2 such that i ≠ j.

In this case, we say that they form a type-2 link.



The minor specifications ({}, {}) and ({2}, {2}) in Example <ref> form a type-1 link, and the minor specifications ({2}, {2}) and ({2}, {3}) in Example <ref> form a type-2 link.


Let S = {2}   {2, 3}
{2}   {2, 3}.
The two minor specifications ({2}, {2}) and ({2, 3}, {2, 3}) form a type-1 link.
The values of m^2 c_S(m, n) suggest m^2 c_S(m, n) = (m + 2 n - m n) (2 m + n - m n).



Let S = {2, 3}   {2, 4}
{2, 3}   {2, 3}, which consists of a type-2 link.
The data suggests m^2 c_S(m, n) = (2 m + n - m n) (2 m + 3 n - m n).



Let S = {2}   {2, 3, 4}
{2}   {2, 3, 4}, whose two elements are not linked.
The data suggests m^2 c_S(m, n) = (m + n - m n) (3 m + 3 n - m n).


These examples, along with others, suggest general formulas for size-2 subsets.


Let m ≥ 1 and n ≥ 1, and let
S =
	
		R_1     R_2

		C_1     C_2⊆ D(m, n).


  * If (R_1, C_1) and (R_2, C_2) form a type-1 link with R_1 + 1 = R_2, then

    m^2 c_S(m, n)
    	= (R_1 (m + n) + n - m n) (R_1 (m + n) + m - m n).


  * If (R_1, C_1) and (R_2, C_2) form a type-2 link with R_1 = R_2, then

    m^2 c_S(m, n)
    	= (R_1 (m + n) - m - m n) (R_1 (m + n) + m - m n).


  * If (R_1, C_1) and (R_2, C_2) form a type-2 link with C_1 = C_2, then

    m^2 c_S(m, n)
    	= (R_1 (m + n) - n - m n) (R_1 (m + n) + n - m n).


  * If (R_1, C_1) and (R_2, C_2) are not linked, then

    m^2 c_S(m, n)
    	= (R_1 (m + n) - m n) (R_2 (m + n) - m n).




Conjectures <ref> and <ref> imply that, if S≤ 2, then the polynomial m^S c_S(m, n) is a product of factors that are linear in m and linear in n.
When we consider subsets S with S≥ 3, we find additional polynomials with this property.
In particular, when S has no linked pairs, there seems to be a simple description of m^S c_S(m, n) as follows.


Let m ≥ 1 and n ≥ 1, and let
S =
	
		R_1     R_2    ⋯    R_k

		C_1     C_2    ⋯    C_k⊆ D(m, n).
If S contains no linked pairs, then

    m^k c_S(m, n)
    	= ∏_i = 1^k (R_i (m + n) - m n).



This conjecture suggests more generally that the structure of the polynomial m^k c_S(m, n) is determined by the connected components of linked pairs in S and, moreover, that the contributions of the components are independent of each other.
This is the content of Conjecture <ref> below.
We make this precise by defining the following graph.


For each S ⊆ D(m, n), let G_S be the graph whose vertex set is S and whose edges connect pairs of linked vertices.



Let S = {}   {2}   {2, 3}   {2, 4}
{}   {2}   {3, 4}   {3, 4}.
The first two elements of S form a type-1 link, and the last two form a type-2 link.
These are the only links, so G_S is the graph on 4 vertices with two non-adjacent edges.
The two connected components are equivalent to the subsets in Examples <ref> and <ref>, respectively.
The values we computed of m^4 c_S(m, n) are as follows.

    [       n = 4     5     6     7; m = 4 -2304 -5040 -7200 -6552;     5 -2880     0;     6     0;     7 10080 ]

This is consistent with m^4 c_S(m, n) = (m - m n) (n - m n) (2 m + n - m n) (2 m + 3 n - m n), which is the product of the formulas in Examples <ref> and <ref>.



Let m ≥ 1 and n ≥ 1.
For every S ⊆ D(m, n),

    m^S c_S(m, n) = ∏_T m^T c_T(m, n)

where the product is over the connected components T of S.
Moreover, since S = ∑_T T, this implies c_S(m, n) = ∏_T c_T(m, n).


Conjecture <ref> is supported by all the rational coefficients we computed.
Assuming it is true, it suffices to determine m^S c_S(m, n) for subsets S consisting of a single connected component.

When a connected component consists of more than one linked pair, the corresponding polynomial is not necessarily a product of linear factors.


Let S = {2}   {3}   {2, 3}
{2}   {3}   {2, 3}.
The minor specifications ({2}, {2}) and ({3}, {3}) are not linked, but each of the other two pairs is.
Therefore G_S consists of a single connected component.
Here are several values of m^3 c_S(m, n):

    [       n = 3     4     5     6; m = 3   -27   -70  -161  -324;     4   -70  -256  -682 -1456;     5  -161  -682 -1875 -4028;     6  -324 -1456 -4028 ]

Without the value for 6 × 6 matrices, we cannot interpolate a cubic polynomial in m and n.
However, by interpolating cubic polynomials in n for the first three rows, we see that m + n - m n is likely a factor.
Dividing each value in the table by this factor, we then interpolate a quadratic polynomial to obtain m^3 c_S(m, n) = (m + n - m n) (2 m^2 + 2 n^2 + 6 m n - 3 m^2 n - 3 m n^2 + m^2 n^2).
This quadratic factor is irreducible.


An obvious question is whether there is a better way to write such polynomials, so that we can see the general structure.
In fact there is, using determinants.
Since the determinant of a block diagonal matrix is the product of the determinants of the blocks, determinant formulas are good candidates for functions that decompose as products over connected components.
In this direction, we next rewrite Conjecture <ref> using determinant formulas; these are likely more natural than the factorizations in Conjecture <ref>.


Let m ≥ 1 and n ≥ 1, and let
S =
	
		R_1     R_2

		C_1     C_2⊆ D(m, n).


  * If (R_1, C_1) and (R_2, C_2) form a type-1 link with R_1 + 1 = R_2, then

    m^2 c_S(m, n)
    	= [ R_1 (m + n) - m n                 m;                -n R_2 (m + n) - m n ].


  * If (R_1, C_1) and (R_2, C_2) form a type-2 link with R_1 = R_2, then

    m^2 c_S(m, n)
    	= [ R_1 (m + n) - m n                -m;                -m R_2 (m + n) - m n ].


  * If (R_1, C_1) and (R_2, C_2) form a type-2 link with C_1 = C_2, then

    m^2 c_S(m, n)
    	= [ R_1 (m + n) - m n                -n;                -n R_2 (m + n) - m n ].


  * If (R_1, C_1) and (R_2, C_2) are not linked, then

    m^2 c_S(m, n)
    	= [ R_1 (m + n) - m n                 0;                 0 R_2 (m + n) - m n ].




For size-1 subsets
S =
	
		R 

		C ⊆ D(m, n), we can rewrite Conjecture <ref> as the 1 × 1 determinant m c_S(m, n) = [ R (m + n) - m n ].
For the size-0 subset S = {}, the definition c_{}(m, n) = 1 is consistent with c_{}(m, n) being the determinant of the 0 × 0 matrix.

Along with Conjecture <ref>, this suggests a determinant formula for m^S c_S(m, n) for an arbitrary subset S.
Since the matrix in this determinant formula resembles an adjacency matrix, we introduce the following notation.


For each
S =
	
		R_1     R_2    ⋯    R_k

		C_1     C_2    ⋯    C_k, define _S(m, n) to be the k × k matrix with the property that, for all i, j satisfying 1 ≤ i < j ≤ k, the 2 × 2 submatrix (_S(m, n))_{i, j}, {i, j} is the matrix in Conjecture <ref> for (R_i, C_i) and (R_j, C_j).
In particular, the ith diagonal entry of _S(m, n) is R_i (m + n) - m n, and the off-diagonal entries are elements of {-m, -n, 0, m}.



As in Example <ref>, let S = {2}   {3}   {2, 3}
{2}   {3}   {2, 3}.
We have

    _S(m, n)
    	= [     m + n - m n               0               m;               0     m + n - m n               m;              -n              -n 2 m + 2 n - m n ].

Indeed, _S(m, n) is equivalent to the formula in Example <ref> for m^3 c_S(m, n).


Unfortunately, this construction doesn't always work.


Let S = {2, 3}   {2, 3}   {2, 3}
{2, 3}   {2, 4}   {2, 5}.
Each pair of elements in S forms a type-2 link, and the common column is the same for all three links.
We have

    _S(m, n)
    	= [ 2 m + 2 n - m n              -m              -m;              -m 2 m + 2 n - m n              -m;              -m              -m 2 m + 2 n - m n ].

The formula _S(m, n) does not produce the following values we computed for m^3 c_S(m, n).

    [       n = 5     6     7     8     9    10    11; m = 3    28    54    80   100   108    98    64;     4   216   256   200     0;     5   500   338;     6   784 ]

However, if we alter the signs of the off-diagonal terms, we can get a determinant formula that produces these values, namely

    m^3 c_S(m, n)
    	= [ 2 m + 2 n - m n               m               m;               m 2 m + 2 n - m n               m;               m               m 2 m + 2 n - m n ].



We conjecture below that the signs of the off-diagonal terms in _S(m, n) can always be altered in such a way, independent of m and n, that its determinant gives the value of m^S c_S(m, n).
Exactly how to alter them is the detail we have not been able to determine.
This alteration is not unique, since if we pick a set V ⊆{1, 2, …, k} of indices and negate the rows and columns indexed by V then the determinant does not change.
An equivalence class of correct sign alterations is generated in this way, all of which lead to the value of m^S c_S(m, n).
This equivalence class seems to depend only on the link structure of S and not on the sizes of its elements.


Let S = {2, 3, 4}   {2, 3, 4}   {2, 3, 4}
{2, 3, 4}   {2, 3, 5}   {2, 3, 6}.
The elements of S are related to each other in the same way as in the previous example;
each pair forms a type-2 link, and the common columns are the same for all three links.
Using the same signs as in the previous example, the determinant formula

    m^3 c_S(m, n)
    	= [ 3 m + 3 n - m n               m               m;               m 3 m + 3 n - m n               m;               m               m 3 m + 3 n - m n ]

agrees with the 4 values we computed.
These same signs also work for S = {2}   {2}   {2}
{2}   {3}   {4}, which also has the same link structure.


Our most general result is therefore Conjecture <ref>, which appears in Section <ref> and which we restate below with more information.
We introduce one last bit of notation to enable sign alterations.
We also take the opportunity to scale the entries of _S(m, n) by 1/m; this allows us to dispense with the factor m^S that would otherwise appear in the summand in Conjecture <ref>.


Let S ⊆ D(m, n), and let k = S.
Let σ{1, 2, …, k}^2 →{-1, 1}.
Define _S, σ(m, n) to be the k × k matrix whose (i, j) entry is 1/mσ((i, j)) times the (i, j) entry of _S(m, n).



Let m ≥ 1 and n ≥ 1.
For each S ⊆ D(m, n), there exists a function σ(S) {1, 2, …, S}^2 →{-1, 1} with σ(S)((i, i)) = 1 for all i such that, for every positive m × n matrix A, the top left entry x of (A) satisfies

    ∑_S ⊆ D(m, n)(_S, σ(S)(m, n)) M(S) x^S = 0.

Moreover, the sign alterations σ(S) can be chosen so that they


  * satisfy σ(T) = σ(S) for all T ≡ S,

  * are independent of m and n, and

  * depend only on the link structure of S and not on the sizes of its elements.



We conclude this section by identifying the correct sign alterations for 2 × n matrices for general n, leading to an explicit formula for the equation satisfied by x.
For each exponent k, there are at most 2 equivalence classes of subsets — an equivalence class containing subsets S such that ({}, {}) ∈ S and another containing the rest.
We consider them separately.


Let m = 2 and S = {}   {2}   {2}   ⋯   {2}
{}   {2}   {3}   ⋯   {k}.
Then

    _S(2, n)
    	= [  -2 n     2     2     ⋯     2;    -n 2 - n    -2     ⋯    -2;    -n    -2 2 - n     ⋯    -2;     ⋮     ⋮     ⋮     ⋱     ⋮;    -n    -2    -2     ⋯ 2 - n ].

If k ∈{1, 2}, then 2^S c_S(2, n) = _S(2, n), so these signs are correct.
If k ≥ 3, they are not, but negating each -2 entry (that is, the off-diagonal entries that aren't in the first row or first column) gives a matrix whose determinant is 2^S c_S(2, n).



Let m = 2 and T = {2}   {2}   {2}   ⋯   {2}
{2}   {3}   {4}   ⋯   {k + 1}.
Then

    _T(2, n)
    	= [ 2 - n    -2     ⋯    -2;    -2 2 - n     ⋯    -2;     ⋮     ⋮     ⋱     ⋮;    -2    -2     ⋯ 2 - n ].

If k ∈{0, 1, 2}, then these signs are correct.
If k ≥ 3, then (again) negating the off-diagonal entries that aren't in the first row or first column gives a matrix whose determinant is 2^S c_T(2, n).


The determinants of the two previous sign-corrected matrices evaluate as follows.


Let n ≥ 1.
Let S_k = {}   {2}   {2}   ⋯   {2}
{}   {2}   {3}   ⋯   {k} for each k ∈{1, …, n - 1, n}, let T_k = {2}   {2}   {2}   ⋯   {2}
{2}   {3}   {4}   ⋯   {k + 1} for each k ∈{0, 1, …, n - 1}, and define the associated coefficients by

    2^k c_S_k(2, n)    =
    	
    		(-n)^k - 1 (2 k - 2 n - 2)	   if 1 ≤ k ≤ n
    
    		0					   if k = 0
    
    	2^k c_T_k(2, n)    =
    	
    		(-n)^k - 1 (2 k - n)	   if 0 ≤ k ≤ n - 1
    
    		0				   if k = n.

(We have defined c_S_0(2, n) = 0 and c_T_n(2, n) = 0 despite S_0 and T_n being undefined; this allows us to write the following sum simply.)
For every positive 2 × n matrix A, the top left entry x of (A) satisfies

    ∑_k = 0^n (
    		c_S_k(2, n) Σ(S_k)
    		+
    		c_T_k(2, n) Σ(T_k)
    	) x^k
    	= 0.



Gröbner basis computations are feasible for 2 × n matrices and establish that the previous conjecture is true for n ≤ 12.




§ OPEN QUESTIONS


The main open question is to identify the equivalence class of correct sign alterations σ(S) in Conjecture <ref> for each S, ideally in the form of a canonical representative.
Since _S(m, n) resembles an adjacency matrix, it seems likely that _S, σ(S)(m, n) has a combinatorial interpretation, and this may suggest the correct signs.
A natural candidate is a signed graph, which is a graph along with a function assigning -1 or 1 to each edge.
The switching class of a signed graph is the set of all signed graphs that can be obtained by choosing a subset V of its vertices and negating the signs of all edges incident to a vertex in V (with multiplicity if an edge is incident to multiple vertices in V) <cit.>.
The equivalence class of correct sign alterations for S corresponds to a switching class on the graph G_S defined in Section <ref>.
We haven't been able to use this to identify the correct signs, however.

Second, how can we prove Conjecture <ref>?
In the absence of a combinatorial proof, one could interpolate polynomial equations representing the entries of diagonal matrices R and C, as we did for the entries of (A), and hope to prove that all three matrices are correct by checking that they satisfy (A) = R A C.
This is how Nathanson established the Sinkhorn limit of a 2 × 2 matrix <cit.>.
However, this works for 2 × 2 matrices because we have explicit expressions for the entries of (A), rather than specifications as roots.
For 3 × 3 matrices, the diagonal entries of R and C are roots of degree-6 polynomials, and generically the product of two such roots has degree 6^2;
this is too big, since the entries of (A) have degree 6.
Therefore we would need to find a single degree-6 field extension that contains all entries of (A), R, and C so that we can perform arithmetic on them symbolically.
What is this extension?

Third, as we saw in Section <ref>, the lack of a unique representation of the coefficient of x^3 for 3 × 3 matrices is due to a relation among 12 monomials M(S) with S = 3.
What is the combinatorial structure of such relations?
Similarly, is there structure in relations among Σ(S)?
For example, the coefficient of x^5 for 3 × 4 matrices has a 1-dimensional family of representations due to the relation

    Σ({}   {2}   {2}   {3}   {2, 3}
    {}   {2}   {3}   {4}   {2, 4})
    	+ Σ({}   {2}   {2}   {2, 3}   {2, 3}
    {}   {2}   {3}   {2, 3}   {2, 4})
    	+ Σ({2}   {2}   {3}   {3}   {2, 3}
    {2}   {3}   {2}   {4}   {2, 3}) 
    
    	+ 2 Σ({}   {2}   {3}   {2, 3}   {2, 3}
    {}   {2}   {3}   {2, 4}   {3, 4})
    	+ 2 Σ({2}   {2}   {2}   {3}   {2, 3}
    {2}   {3}   {4}   {2}   {3, 4})
    	+ 2 Σ({2}   {2}   {3}   {2, 3}   {2, 3}
    {2}   {3}   {4}   {2, 4}   {3, 4}) 
    
    	=
    	Σ({}   {2}   {3}   {2, 3}   {2, 3}
    {}   {2}   {3}   {2, 3}   {2, 4})
    	+ Σ({2}   {2}   {2}   {3}   {2, 3}
    {2}   {3}   {4}   {2}   {2, 3})
    	+ Σ({2}   {2}   {3}   {2, 3}   {2, 3}
    {2}   {3}   {4}   {2, 3}   {2, 4}) 
    
    	+ 2 Σ({}   {2}   {2}   {3}   {2, 3}
    {}   {2}   {3}   {4}   {2, 3})
    	+ 2 Σ({}   {2}   {2}   {2, 3}   {2, 3}
    {}   {2}   {3}   {2, 4}   {3, 4})
    	+ 2 Σ({2}   {2}   {3}   {3}   {2, 3}
    {2}   {3}   {2}   {4}   {3, 4}).

It would be interesting to understand these better.

Fourth, how does Corollary <ref> generalize to m × n matrices with linear dependencies among their rows or columns?
The coefficients in Corollary <ref> have unique representations as linear combinations of class sums Σ(S) where S ⊆ D(2, 3), namely

    e_3    =
    		Σ({}   {2}   {2}
    {}   {2}   {3}) 
    
    	e_2    =
    		-2 Σ({}   {2}
    {}   {2})
    		+ Σ({2}   {2}
    {2}   {3}) 
    
    	e_1    =
    		3 Σ({}
    {}) 
    
    	e_0    =
    		-Σ(
    ).

What is the general formula?

Finally, there is a further generalization of Conjecture <ref> whose form is not known.
Decades before Sinkhorn's paper, the iterative scaling process was introduced by Kruithof <cit.> in the context of predicting telephone traffic.
In this application, rather than scaling to obtain row and column sums of 1, each row and column has a potentially different target sum.
Sinkhorn <cit.> showed that the limit exists.
We call this limit the Kruithof limit.


Kruithof <cit.> considered the matrix

    A =
    	[ 2000 1030  650  320; 1080 1110  555  255;  720  580  500  200;  350  280  210  160 ]

with target row sums V = [ 6000 4000 2500 1000 ] and target column sums W = [ 6225 4000 2340  935 ].
Let x be the top left entry of the Kruithof limit.
Numerically, x ≈ 3246.38700234.
A Gröbner basis computation gives an equation

    62   11170485642866385308015185014605806684592592997303612 x^20
       - 1911288675240357642608985257264441863326549355081446688219995 x^19
       + ⋯
       - 980316295756763597938629190043558577216660563425441394040234375 · 10^72 x 
       + 60077293526471262201893650291744622440239260152893558984375 · 10^79
    	= 0

satisfied by x.
In particular, x has degree 20.


The previous example suggests that the degrees of entries of Kruithof limits are the same as those of Sinkhorn limits of the same size.


Let m ≥ 1 and n ≥ 1.
Let A be a positive m × n matrix, let V be a positive m × 1 matrix, and let W be a positive 1 × n matrix such that the sum of the entries of V equals the sum of the entries of W.
The top left entry x of the Kruithof limit of A with target row sums V and target column sums W is algebraic over the field generated by the entries of A, V, and W, with degree at most m + n - 2m - 1.


Since the Kruithof limit specializes to the Sinkhorn limit when V = [ 1 1 ⋯ 1 ] and W = [ m/n m/n   ⋯ m/n ], we suspect that the entries in the determinant formulas in Conjecture <ref> generalize somehow to involve entries of V and W, and this should give a generalization of Conjecture <ref> to Kruithof limits.
Moreover, Gröbner basis computations suggest that the surprising property we mentioned in Section <ref> regarding the coefficients for square matrices satisfying

    ∑_S ⊆ D(n, n) c_S(n, n) x^S = (x - 1)^2 n - 2n - 1

generalizes as follows.


With the notation of Conjecture <ref>, let V = [ r_1 r_2   ⋯ r_m ] and W = [ c_1 c_2   ⋯ c_n ].
If the general equation satisfied by the top left entry x of the Kruithof limit of A with target row sums V and target column sums W is

    ∑_S ⊆ D(m, n) c_S(V, W) M(S) x^S = 0,

then

    ∑_S ⊆ D(m, n) c_S(V, W) x^S
    	= (x - r_1)^m + n - 3n - 1(x - c_1)^m + n - 3m - 1.

In particular, this sum is independent of r_2, …, r_m and c_2, …, c_n.



99

Allen-Zhu–Li–Oliveira–Wigderson
Zeyuan Allen-Zhu, Yuanzhi Li, Rafael Oliveira, and Avi Wigderson,
Much faster algorithms for matrix scaling,
58th Annual IEEE Symposium on Foundations of Computer Science (2017) 890–901.

Brown
David T. Brown,
A note on approximations to discrete probability distributions,
Information and Control 2 (1959) 386–392.

Chen–Varghese
Kevin Chen and Abel Varghese, personal communication, Summer 2019.

Cohen–Madry–Tsipras–Vladu
Michael B. Cohen, Aleksander Madry, Dimitris Tsipras, and Adrian Vladu,
Matrix scaling and balancing via box constrained Newton's method and interior point methods,
58th Annual IEEE Symposium on Foundations of Computer Science (2017) 902–913.

Cuturi
Marco Cuturi,
Sinkhorn distances: lightspeed computation of optimal transportation distances,
Advances in Neural Information Processing Systems 26 (2013) 2292–2300.

Deming–Stephan
W. Edwards Deming and Frederick F. Stephan,
On a least squares adjustment of a sampled frequency table when the expected marginal totals are known,
The Annals of Mathematical Statistics 11 (1940) 427–444.

Ekhad–Zeilberger
Shalosh B. Ekhad and Doron Zeilberger,
Answers to some questions about explicit Sinkhorn limits posed by Mel Nathanson,
<https://arxiv.org/abs/1902.10783> (6 pages).

Ferguson–Bailey
Helaman R. P. Ferguson and David H. Bailey,
A polynomial time, numerically stable integer relation algorithm,
RNR Technical Report RNR-91-032 (1992) (14 pages).

Franklin–Lorenz
Joel Franklin and Jens Lorenz,
On the scaling of multidimensional matrices,
Linear Algebra and its Applications 114–115 (1989) 717–735.

Idel
Martin Idel,
A review of matrix scaling and Sinkhorn's normal form for matrices and positive maps,
<https://arxiv.org/abs/1609.06349> (101 pages).

Kalantari–Khachiyan 1
Bahman Kalantari and Leonid Khachiyan,
On the rate of convergence of deterministic and randomized RAS matrix scaling algorithms,
Operations Research Letters 14 (1993) 237–244.

Kalantari–Khachiyan 2
Bahman Kalantari and Leonid Khachiyan,
On the complexity of nonnegative-matrix scaling,
Linear Algebra and its Applications 240 (1996) 87–103.

Kalantari–Lari–Ricca–Simeone
B. Kalantari, I. Lari, F. Ricca, and B. Simeone,
On the complexity of general matrix scaling and entropy minimization via the RAS algorithm,
Mathematical Programming 112 (2008) 371–401.

Kruithof
J. Kruithof,
Telefoonverkeersrekening,
De Ingenieur 52 (1937) E15–E25.
English translation by Pieter-Tjerk de Boer:
<https://wwwhome.ewi.utwente.nl/ ptdeboer/misc/kruithof-1937-translation.html>

Linial–Samorodnitsky–Wigderson
Nathan Linial, Alex Samorodnitsky, and Avi Wigderson,
A deterministic strongly polynomial algorithm for matrix scaling and approximate permanents,
Combinatorica 20 (2000) 545–568.

Nathanson 2x2
Melvyn B. Nathanson,
Alternate minimization and doubly stochastic matrices,
Integers 20A (2020) Article #A10 (17 pages).

Nathanson 3x3
Melvyn B. Nathanson,
Matrix scaling and explicit doubly stochastic limits,
Linear Algebra and its Applications 578 (2019) 111–132.

SinkhornPolynomials
Eric Rowland and Jason Wu,
SinkhornPolynomials,
a Mathematica package available as an ancillary file associated with this arXiv submission.

Sinkhorn
Richard Sinkhorn,
A relationship between arbitrary positive matrices and doubly stochastic matrices,
The Annals of Mathematical Statistics 35 (1964) 876–879.

Sinkhorn 1967
Richard Sinkhorn,
Diagonal equivalence to matrices with prescribed row and column sums,
The American Mathematical Monthly 74 (1967) 402–405.

OEIS
Neil Sloane et al.,
The On-Line Encyclopedia of Integer Sequences,
<http://oeis.org>.

GroebnerBasis
Wolfram Research,
GroebnerBasis,
Wolfram Language function,
<https://reference.wolfram.com/language/ref/GroebnerBasis.html> (1991, updated 2007).

RootApproximant
Wolfram Research,
RootApproximant,
Wolfram Language function,
<http://reference.wolfram.com/language/ref/RootApproximant.html> (2007, updated 2008).

Zaslavsky
Thomas Zaslavsky,
Characterizations of signed graphs,
Journal of Graph Theory 5 (1981) 401–406.




