Solving the multi-dimensional scaling problem in CVXPY

convex optimizationcvxpyoptimizationpythonsemidefinite-programming

I want to use the CVXPY Python toolbox to solve the simple classical multi-dimensional scaling problem.

Problem: Given $n$ points in $d$-dimensional space, i.e. $\mathbf{x}_i \in \mathbb{R}^d$. The position coordinates of the points are stacked in $X \in \mathbb{R}^{d \times n}$. The Euclidean Distance Matrix, i.e. the matrix containing the square of point-wise distance between the $n$ points is given as $D \in \mathbb{R}^{n \times n}$. The problem is to recover $X$ given $D$, i.e.

$$\hat{X} = \arg \min_{X} \| D – X^T X \|_F^2$$

MDS solution: Without going into too much details, a solution can be obtained using eigenvalue decomposition, i.e.

$$\hat{X} = \Lambda^{1/2} V^T$$

where $\Lambda$ is a diagonal matrix containing largest $d$ eigenvalues of centered matrix $\dfrac{1}{2} J D J$, with $J = I – \dfrac{\mathbf{1} \mathbf{1}^T}{n}$ and $V$ is the matrix with the corresponding eigenvectors.

My Problem:I want to formulate it as an SDP and this is how I do it:

\begin{align}
\hat{Z} &= \arg \min_{Z} \quad \textbf{tr} (D – Y)\\
\text{s.t.} & \begin{pmatrix}
Y & Z^T\\
Z & I
\end{pmatrix} \geq 0 \\
& \quad Y \geq 0
\end{align}

I am a bit unsure if the second constraint is needed though.

Question 1: Is this SDP formulation correct? I am getting $\infty$ as the solution.

Note: The part below has already been solved. The code I used is also laid put for anyone who wants to try.

Question 2: This is how I wrote my python code:

import numpy as np
import cvxpy as cvx

nDim = 2
n = 4

def edm(X):
    d, n = X.shape

    one = np.ones((n, 1))
    G = X.transpose().dot(X)
    g = G.diagonal().reshape((n, 1))
    D = g.dot(one.transpose()) + one.dot(g.transpose()) - 2.0 * G

    return D

X = np.array([[0.0, 0.0], [10.0, 5.0], [10.0, 20.0], [0.0, 10.0]]).transpose()
D = edm(X)

# Setting up the cvx problem
Y = cvx.Variable((n, n), symmetric=True)
Z = cvx.Variable((nDim, n))

# Schur complement
I = np.identity(nDim)
M = np.block([[Y, Z.T], [Z, I]])

constraints = [M >> 0]
constraints += [Y >> 0]

prob = cvx.Problem(cvx.Minimize(cvx.trace(D - Y)), constraints)
prob.solve()

I get the following error:

ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 1 and the array at index 1 has size 2

I have checked the dimensions and there should not be any problem with the block matrix inequality in terms of dimensions of $Y$, $Z$ and $I$.

I will appreciate if someone can help me with the syntax or the right formulation. Let me know if more information is needed!

Edit #1: The second question is syntax related and I found the solution here: Github Issue. Though the solver is not able to solve it and gives Inf.

Edit #2: Added more context to make it accessible to more people.

Edit #3: I am accepting the answer from @LinAlg as this user helped to formulate the problem correctly which was one of the intents of the question. The other intent was to directly get the coordinates from the CVX problem. In my short search, eigenvalue decomposition is the only way to get some coordinates, as is also indirectly pointed out by @LinAlg. But I will always encourage other people to provide/share a solution/nuance that helps everyone.

Best Answer

  1. You need $\text{tr}(Y-D)$ as the objective (currently you can let $Y \to \infty I$). Note that $D$ is just a constant you can omit so currently you do not use your input data. You should add the constraints from slide 10 (or 12). About the constraint on $Y$: if a matrix is positive semidefinite, principal submatrices are also positive semidefinite (which is easy to show using the $x^TAx \geq 0$ definition of psd).

  2. Use cvx.bmat instead of np.block.

Related Question