1

Given a matrix $\boldsymbol{Y} \in \mathbb{R}^{m \times n}$ where $n \leq m$.
I want to find its projection onto the Convex Hull of Permutation Matrices:

$$ \mathcal{P} = \left\{ \boldsymbol{P} \mid \boldsymbol{P} \boldsymbol{1} \leq \boldsymbol{1}, \boldsymbol{1}^{T} \boldsymbol{P} = \boldsymbol{1}, \boldsymbol{P} \geq 0 \right\} $$

Namely, each matrix is a set of columns of a Permutation Matrix (More accurately, from the convex hull of the permutation matrices).

The actual projection is given by:

$$\begin{align*} \arg \min_{\boldsymbol{X}} \quad & \frac{1}{2} {\left\| \boldsymbol{X} - \boldsymbol{Y} \right\|}_{F}^{2} \\ \text{subject to} \quad & \begin{aligned} \boldsymbol{X} & \in \mathcal{P} \end{aligned} \end{align*}$$

Is there an efficient algorithm to find such $\boldsymbol{X}$?

I thought about 2 approaches:

  1. Least Squares with Linear Constraints (See MATLAB's lsqlin()).
  2. Projection onto an intersection of convex sets (See Orthogonal Projection onto the Intersection of Convex Sets).

Yet, it seems to me there should be something more efficient to take advantage of the structure of the set.

There is the Sinkhorn Algorithm which solve the problem in the case of a square matrix and equality $\boldsymbol{P} \boldsymbol{1} = \boldsymbol{1}$.

Royi
  • 10,050
  • Might be related to https://math.stackexchange.com/questions/4315197. – Royi Jun 23 '24 at 19:27
  • The set $\mathcal P$ contains more than just projection matrices. Somehow the condition $P_{ij} \in {0,1}$ is missing. – daw Jun 26 '24 at 06:50

2 Answers2

1

For now, I implemented the solution using DCP Solver (Convex.jl) and Dykstra's Projection Algorithm.

My naive Julia implementation of Dykstra's Projection Algorithm was order of magnitude faster than the DCP bases solver (Using the ECOS solver).

The projections are given by:

  • $ \mathcal{C} = \left\{ \boldsymbol{X} \in \mathbb{R}^{m \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}_{\mathbb{I}}^{T} $ where $\boldsymbol{1}_{\mathbb{I}}^{T}$ is a vector with value $1$ where the sum of the row of $ \boldsymbol{Y}$ is larger than 1 and $0$ elsewhere.
  • $ \mathcal{C} = \left\{ \boldsymbol{X} \in \mathbb{R}^{m \times n} \mid \boldsymbol{1}^{T} \boldsymbol{X} = \boldsymbol{1} \right\} \Rightarrow {P}_{\mathcal{C}} \left( \boldsymbol{Y} \right) = \boldsymbol{Y} - \frac{1}{m} \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.

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

Remarks:

  • I will try to add another solution based on RipQP.jl.
  • The implementation can be farther optimized. Mainly to reduce allocations.
  • I still think some clever variation of the Sinkhorn algorithm could be utilized here to make the process much faster.
Royi
  • 10,050
0

Are you looking for a custom solver? Because if you're not, any QP solver will do. Even the high-level interface of CVXPY will.

inmport cvxpy as cp

Y = ...
X = cp.Variable(Y.shape)
objective = cp.sum_squares(Y - X)
constraints = [
    cp.sum(X, axis=1) <= 1,
    cp.sum(X, axis=0) == 1,
    X >= 0]
prob = cp.Problem(cp.Minimize(objective), constraints)
prob.solve()
print(X.value)
Royi
  • 10,050
Alex Shtoff
  • 168
  • I am looking for something which is specialized on this case. For instance lsqlin() will be faster than CVXPY on this. But I'm pretty sure one can be much faster. – Royi Jun 24 '24 at 07:08