# Project an orthogonal matrix onto the Birkhoff Polytope

It is known that the permutation matrices lie at the intersection of the orthogonal group $$\mathbb{O}^N$$ with the Birkhoff polytope $$\mathbb{DS}^N$$. It is also known that any non-negative matrix $$X\in\mathbb{R}^{N \times N}$$ can be projected onto $$\mathbb{DS}^N$$ by the $$\textit{Sinkhorn algorithm}$$. Finally, any matrix in $$\mathbb{O}^N$$ has to be a permutation matrix, if all the entries are positive (matrix is non-negative). This also means that, any orthogonal matrix $$O$$ will contain negative entries if it is not a permutation.

My qusetion is, how can I project an orthogonal matrix, which is allowed to contain negative entries onto $$\mathbb{DS}$$?

I see that this is a relatively uncharted territory. My answer to this will be Randomized Rounding [1]. I am not sure whether the cited resource completely describes the algorithm I came up with, so I will explain it here:

Let $$\mathbf{O}$$ denote the input orthogonal matrix, which we wish to project onto a point on the Birkhoff Polytope $$\mathbf{X} \in \mathcal{DS}$$. Let $$(\mathbf{x}, \mathbf{y})$$ be vectors in $$\mathbb{R}^d$$ with $$\mathbf{y}=\mathbf{O}\mathbf{x}$$ containing distinct coordinates. Let $$\mathbf{x}$$ be sampled at random from the standard Gaussian measure on $$\mathbb{R}^d$$. Barvinok observed that the action of a given orthogonal matrix on a random vector $$\mathbf{x}$$ with high probability is very close to a permutation of the coordinates:
The rounding of $$\mathbf{O}$$ at $$\mathbf{x}$$ is defined as the permutation $$\mathbf{P}_{(\mathbf{O},\mathbf{x})}=\sigma(\mathbf{O}, \mathbf{x})$$ s.t. $$\mathbf{P}_{(\mathbf{O},\mathbf{x})}$$ matches the $$k^{th}$$ smallest coordinate of $$\mathbf{x}$$ with the $$k^{th}$$ smallest coordinate of $$\mathbf{y}=\mathbf{O}\mathbf{x}$$. Yet this permutation varies with $$\mathbf{x}$$. To arrive at the $$\mathbf{x}$$-invariant projection onto $$\mathcal{DS}_d$$, one can easily construct small approximate non-commutative convex combinations using the Birkhoff von Neumann theorem [2]: $$$$(\mathbf{X} \in \mathcal{DS}_d) =\frac{1}{T}\sum\limits_i^{T}\mathbf{P}_{(\mathbf{O},\mathbf{x}_i)}$$$$ where $$T$$ is the number of iterates - the higher the more accurate. Finally, the permutation matrix corresponding to $$\mathbf{X}\in\mathcal{DP}_d$$ can be found via the famous Hungarian algorithm.

I am providing a sample MATLAB code of the procedure summarized above:

function DS = ON_2DS(O, iterations)

n = size(O,1);
if(~exist('iterations','var'))
iterations = 1000;
end

DS = zeros(size(O));
for i=1:iterations
x = randn(n,1);
%x=x./norm(x);
y = O*x;

[~,indx] = sort(x,'ascend');
[~,indy] = sort(y,'ascend');

ind = sub2ind(size(O), indy, indx);
P = zeros(size(O));
P(ind) = 1;
DS = DS + P;
end

DS = DS./iterations;

end


References
[1] Approximating orthogonal matrices by permutation matrices, Alexander Barvinok, 2005
[2] A SHORT PROOF OF THE BIRKHOFF-VON NEUMANN THEOREM, GLENN HURLBERT, 2008 http://www.people.vcu.edu/~ghurlbert/papers/SPBVNT.pdf