SinkhornPolynomials 1.00

documentation notebook

Eric Rowland and Jason Wu
https://ericrowland.github.io/packages.html

Introduction

SinkhornPolynomials is a package for computing the polynomial equation satisfied by an entry of the Sinkhorn limit of a matrix.
This introduction gets you started with a few features of the package; the next section provides a complete list of package symbols along with their usage messages and further examples.

To use SinkhornPolynomials, first you will need to load the package by evaluating the following cell. (If you need help, see loading a package.)

<<SinkhornPolynomials`

Square matrices

The function SinkhornFixedPoint computes an approximation to the Sinkhorn limit of a matrix:

SinkhornFixedPoint[{{5,5,4},{3,8,5},{9,3,2}}]//MatrixForm

The function SinkhornPolynomial gives a polynomial whose roots include the top left entry of the Sinkhorn limit:

SinkhornPolynomial[{{5,5,4},{3,8,5},{9,3,2}},x]
Factor[%]

The top left entry of the Sinkhorn limit is one of these two numbers:

Select[NSolveValues[%==0,x],#∈Reals&&0<#<1&]

For 4×4 matrices, unpacking the coefficient data takes several seconds the first time:

EchoTiming@SinkhornPolynomial[{{3,1,2,2},{2,2,2,1},{1,2,3,2},{1,4,2,3}},x]
Factor[%]

Subsequent uses are faster:

EchoTiming@SinkhornPolynomial[{{3,1,2,2},{2,2,2,1},{1,2,3,2},{1,4,2,3}},x]

The coefficients can be computed one at a time with SinkhornPolynomialCoefficient:

EchoTiming@SinkhornPolynomialCoefficient[{{3,1,2,2},{2,2,2,1},{1,2,3,2},{1,4,2,3}},3]

The Sinkhorn polynomial and its coefficients can be computed for matrices with symbolic entries:

SinkhornPolynomialCoefficient[{{a11,a12,a13},{a21,a22,a23},{a31,a32,a33}},2]

Row and column scaling

By default, rows are scaled so the row sums are 1.
A different row sum can be specified as an additional argument:

SinkhornFixedPoint[{{3,9,1},{3,2,9},{5,3,4}},(1/3)]//MatrixForm
Total/@%

SinkhornPolynomial[{{3,9,1},{3,2,9},{5,3,4}},x,(1/3)]
Select[NSolveValues[%==0,x],#∈Reals&&0<#<(1/3)&]

Rectangular matrices

A 3×4 matrix:

SinkhornPolynomial[{{9,5,4,8},{9,6,1,5},{5,1,2,6}},x]
NSolveValues[%==0,x]

SinkhornFixedPoint[{{9,5,4,8},{9,6,1,5},{5,1,2,6}}]

The polynomials computed by SinkhornPolynomial have been proven to be correct for 2×n matrices where n<=12 (and their transposes) and for 3×3 matrices.
For 2×n (where n>=13), 3×4, 3×5, and 4×4 matrices (and their transposes), the polynomials computed by SinkhornPolynomial are conjecturally correct.
Larger matrices are not supported because the coefficients for the corresponding general matrix have never been computed.

SinkhornPolynomial[{{6,2,9,8,2},{7,2,4,5,9},{4,7,5,6,9},{6,9,6,7,7}},x]

With sufficiently high precision, an entry of a Sinkhorn limit can be recognized by RootApproximant:

EchoTiming@Module[{matrix={{6,2,9,8,2},{7,2,4,5,9},{4,7,5,6,9},{6,9,6,7,7}},precision=(2^13),rowcount,columncount},
{rowcount,columncount}=Dimensions[matrix];
MinimalPolynomial[RootApproximant[SinkhornFixedPoint[N[matrix,precision]][[1,1]],Binomial[rowcount+columncount-2,rowcount-1]],x]
]

Package symbols

SinkhornFixedPoint

?SinkhornFixedPoint

SinkhornFixedPoint[{{1,4,9},{16,25,36},{49,64,81}}]//MatrixForm

To get a fixed point with higher precision, numericize the entries with higher precision:

SinkhornFixedPoint[N[{{1,4,9},{16,25,36},{49,64,81}},50]]//MatrixForm

By default, the row sums are 1, so for non-square m×n matrices the column sums are m/n:

SeedRandom[1234];
MatrixForm[matrix=RandomInteger[{1,9},{5,8}]]

MatrixForm[fixedpoint=SinkhornFixedPoint[matrix]]

Total/@fixedpoint
Total/@Transpose[fixedpoint]

Scale so that the row sums are 1/m and the column sums are 1/n:

MatrixForm[scaledfixedpoint=SinkhornFixedPoint[matrix,(1/5)]]

Total/@scaledfixedpoint
Total/@Transpose[scaledfixedpoint]

The two fixed points differ by a scalar multiple:

MatrixForm[fixedpoint/scaledfixedpoint]

Clear[matrix,fixedpoint,scaledfixedpoint]

The matrix entries are allowed to be negative:

SinkhornFixedPoint[{{-7,9},{6,-5}}]//MatrixForm

However, negative entries can result in division by 0:

SinkhornFixedPoint[{{-6,8},{3,1}}]//MatrixForm

SinkhornPolynomial

?SinkhornPolynomial

Compute the Sinkhorn polynomial for a general 2×2 matrix:

SinkhornPolynomial[{{a,b},{c,d}},x]

Find the roots of the Sinkhorn polynomial for a specific 3×3 matrix:

SinkhornPolynomial[{{1,4,9},{16,25,36},{49,64,81}},x]
NSolveValues[%==0,x]

Change the row scaling:

SinkhornPolynomial[{{1,4,9},{16,25,36},{49,64,81}},x,(1/3)]
NSolveValues[%==0,x]

This is equivalent to replacing x->3x:

SinkhornPolynomial[{{1,4,9},{16,25,36},{49,64,81}},3x]

The Sinkhorn polynomial is 0 if certain rows or columns are scalar multiples of each other:

SinkhornPolynomial[{{4,5,6},{1,2,3},{2,4,6}},x]

A nonzero polynomial for the top left entry of the Sinkhorn limit can be obtained by carefully factoring out the determinants that are 0, as follows.

With[{matrix={{a11,a12,a13},{a21,a22,a23},{a31,a32,a33}}},With[{determinant=Det[matrix],minor=Det[matrix[[{2,3},{2,3}]]]},Factor[Cancel[(1/minor)(SinkhornPolynomial[matrix,x]/.determinant->detratio minor)]/.Thread[matrix[[3]]->rowscale matrix[[2]]]]]]
Collect[FactorList[%],x,Factor]
%[[4,1]]

Evaluate the polynomial for a particular matrix:

-(a11^2) a22 a23+3 (a11^2) a22 a23 x+(a12 a13 (a21^2)+a11 a13 a21 a22+a11 a12 a21 a23-3 (a11^2) a22 a23) (x^2)+(a12 a21-a11 a22) (a13 a21-a11 a23) (x^3)/.Thread[Join@@Take[{{a11,a12,a13},{a21,a22,a23},{a31,a32,a33}},2]->Join@@Take[{{4,5,6},{1,2,3},{2,4,6}},2]]
NSolveValues[%==0,x]

One of the roots is the top left entry of the Sinkhorn limit:

SinkhornFixedPoint[{{4,5,6},{1,2,3},{2,4,6}}]//MatrixForm

By default, free variables that arose when solving for the coefficients are set to 0.
They can be shown by specifying a function for GeneratedParameters:

SinkhornPolynomial[{{9,5,4,8},{9,6,1,5},{5,1,2,6}},x,GeneratedParameters->C]

The free variables cancel each other out when expanded:

Expand[%]

SinkhornPolynomialCoefficient

?SinkhornPolynomialCoefficient

SinkhornPolynomialCoefficient[{{1,4,9},{16,25,36},{49,64,81}},4]

SinkhornPolynomialCoefficient[{{1,4,9},{16,25,36},{49,64,81}},4,(1/3)]