Projection onto Birkhoff Polytope

Question

Suppose we would like to compute the Euclidean projection of an arbitrary matrix $A$ onto the Birkhoff polytope, the set of doubly-stochastic matrices.

1. Under some conditions on $A$, Sinkhorn's algorithm returns two diagonal matrices $D_1,D_2$ such that $D_1 A D_2$ is a doubly-stochastic matrix (which is unique). How does this matrix compare to the Euclidean projection, and if not, does this solution correspond to the projection under any particular metric?
2. Are there known methods for solving this efficiently besides just formulating the problem as a quadratic program and applying a generic QP algorithm?

Answer 1

It is known that symmetric Sinkhorn algorithm in fact minimizes KL divergence [2,3]. In 1, authors present a method to minimize the Euclidean distance. This is called BBS (Bregmanian Bi-Stochastication). It is reported that BBS algorithm using the Euclidian distance has noticeably better potential of producing good clustering results than the SK algorithm, while Sinkhorn performs better in terms of doubly stochasticity.

(1) Learning a Bi-Stochastic Data Similarity Matrix, Fei Lang, Ping Li, Arnd Christian Konig

(2) J. N. Darroch and D. Ratcliff. Generalized iterative scaling for log-linear models. The Annals of Mathematical Statistics, 43(5):1470–1480, 1972.

(3) G. W. Soules. The rate of convergence of Sinkhorn balancing. Linear Algebra and its Applications, 150:3 – 40, 1991.

Regarding the projection of an arbitrary (arguably non-positive matrix): First, it is quite hard to project a matrix with negative entries. I would first project it onto the orthogonal matrices and then from there round it to the Birkhoff polytope. It is not a great procedure though, but a reasonable approach is given in :

(4) Approximating Orthogonal Matrices by Permutation Matrices, Alexander Barvinok, 2005

(5) Reparameterizing the Birkhoff Polytope for Variational Permutation Inference - See Doubly Stochastic Matrices in Stan.

Answer 2

Defining the problem as:

$$ \begin{align*} \arg \min_{\boldsymbol{X}} \quad & \frac{1}{2} {\left\| \boldsymbol{X} - \boldsymbol{Y} \right\|}_{2}^{2} \\ \text{subject to} \quad & \begin{aligned} \boldsymbol{X} \boldsymbol{1} & = \boldsymbol{1} \\ {X}_{i, j} & \geq 0 \\ \boldsymbol{X} & = \boldsymbol{X}^{T} \end{aligned} \end{align*} $$

One could use Orthogonal Projection onto the Intersection of Convex Sets in order to calculate the projection.

Each projection on its own is pretty trivial:

• $ \mathcal{C} = \left\{ \boldsymbol{X} \in \mathbb{R}^{n \times n} \mid \boldsymbol{X} \boldsymbol{1} \leq \boldsymbol{1} \right\} \Rightarrow {P}_{\mathcal{C}} \left( \boldsymbol{Y} \right) = \boldsymbol{Y} - \frac{1}{n} \left( \boldsymbol{Y} \boldsymbol{1} - \boldsymbol{1} \right) \boldsymbol{1} $.
• $ \mathcal{C} = \left\{ \boldsymbol{X} \in \mathbb{R}^{n \times n} \mid \boldsymbol{1}^{T} \boldsymbol{X} = \boldsymbol{1} \right\} \Rightarrow {P}_{\mathcal{C}} \left( \boldsymbol{Y} \right) = \boldsymbol{Y} - \frac{1}{n} \boldsymbol{1} \left( \boldsymbol{1}^{T} \boldsymbol{Y} - \boldsymbol{1}^{T} \right) $.
• $ \mathcal{C} = \left\{ \boldsymbol{X} \in \mathbb{R}^{m \times n} \mid \boldsymbol{X}_{i, j} \geq 0 \right\} \Rightarrow {P}_{\mathcal{C}} \left( \boldsymbol{Y} \right) = \max \left\{ \boldsymbol{Y}, \boldsymbol{0} \right\} $ where $\boldsymbol{0}$ is a matrix of zeros.
• $ \mathcal{C} = \left\{ \boldsymbol{X} \in \mathbb{R}^{n \times n} \mid \boldsymbol{X} = \boldsymbol{X}^{T} \right\} \Rightarrow {P}_{\mathcal{C}} \left( \boldsymbol{Y} \right) = 0.5 \left( \boldsymbol{Y} + \boldsymbol{Y}^{T} \right) $.

The full code is available on my StackExchange Mathematics GitHub Repository (Look at the Mathematics\Q838813 folder).

Remarks:

• The definition of Doubly Stochastic Matrix does not require symmetry.
• I used symmetry to follow the Learning a Bi-Stochastic Data Similarity Matrix paper. In the paper context they require similarity as the matrix origins is similarity matrix.
• The paper applies alternating projection (Although claims to use Dykstra's Projection Algorithm) though the sets are not sub spaces. It skipped proving the convergence to optimality. My code use Dykstra's Projection Algorithm.
• I apply projection both onto rows and columns as it made convergence faster.

Answer 3

Besides the Sinkhorn-Knopp algorithm, I believe it is possible to apply the alternating projection algorithm to the matrix $A$. Specifically, I use the projection onto the probability simplex as described by Wang et al. (2013), which is easy to implement. The process involves alternately projecting each row and then each column onto the probability simplex until convergence. Empirically, it converges very quickly and yields results that are close to those of Dykstra's algorithm, which guarantees the exact projection onto the intersection of convex sets.