title: Semidefinite embedding in python based on cvxopt
body: <p>I am trying to implement the <strong>semidefinite embedding</strong> algorithm (see <a href="https://en.wikipedia.org/wiki/Semidefinite_embedding" rel="nofollow">here</a>) in python based on the package <strong>cvxopt</strong> for solving semidefinite programming. I am having some problems mapping the definition of the semidefinite program to cvxopt's interface (see <a href="http://abel.ee.ucla.edu/cvxopt/userguide/coneprog.html#semidefinite-programming" rel="nofollow">this</a>).</p>
  
  <p>This is my current implementation:</p>
  
  <pre><code>from numpy import array, eye
  from numpy.linalg import norm
  from numpy.random import multivariate_normal
  from sklearn.neighbors import NearestNeighbors
  from cvxopt import matrix, spmatrix, solvers
  
  n = 10 # The number of data points
  
  #### Generate data and determine nearest neighbor structure ####
  # Sample data
  X = multivariate_normal([0, 0], [[10, -10], [2, -1]] ,n)
  
  # Determine 5-nearest neighbors graph
  N = NearestNeighbors(n_neighbors=5).fit(X).kneighbors_graph(X).todense()
  N = array(N)
  
  # Generate set of index-pairs, where each pair corresponds to a eighbor pair
  neighs = set(zip(*N.nonzero()))
  neighs.union(set(zip(*(N.T.dot(N)).nonzero())))
  
  # Plot data points and nearest neighbor graph
  #import pylab
  #for i, j in neighs:
      #pylab.plot([X[i, 0], X[j, 0]], [X[i, 1], X[j, 1]], 'k')
  #pylab.plot(X[:, 0], X[:, 1], 'ro')
  #pylab.show()
  
  #### Create kernel using semidefinite programming ####
  # We want to maximize the trace, i.e. minimize the negative of the trace
  c = matrix(-1 * eye(n).reshape(n**2, 1))
  
  ## x must be positive semidefinite
  # Note -Gs*x should be the n*n matrix representation of x (which is n**2, 1)
  Gs = [spmatrix([-1.0]*n**2, range(n**2), range(n**2), (n**2, n**2))]
  hs = [spmatrix([], [], [], (n, n))]  #  Zero matrix
  
  ## Equality constraints
  # Kernel must be centered, i.e. sum of all kernel elements needs to be 0
  A = [[1] * (n**2)]
  b = [[0.0]]
  
  # Add one row to A and b for each neighbor distance-preserving constraint
  for i, j in neighs:
      if i &gt; j and (j, i) in neighs: 
          continue # skip since this neighbor-relatonship is symmetric
      i = int(i)
      j = int(j)
      # TODO: Use sparse matrix for A_
      #spmatrix([1.0, -2.0, 1.0], [0, 0, 0], [i*(n + 1), i*n + j, j*(n + 1)], 
      #        (1, n**2))
      A_ = [0.0 for k in range(n**2)]
      A_[i*n + i] = +1 # K_ii
      A_[i*n + j] = -2 # K_ij
      A_[j*n + j] = +1 # K_jj
      b_ = [norm(X[i] - X[j])**2] # squared distance between points with index i and j
      A.append(A_)
      b.append(b_)
  
  # Solve positive semidefinite program
  A = matrix(array(A, dtype=float))
  b = matrix(array(b, dtype=float))
  
  sol = solvers.sdp(c=c, Gs=Gs, hs=hs, A=A, b=b)
  </code></pre>
  
  <p>Executing the program throws a <code>ValueError</code> (Rank(A) &lt; p or Rank([G; A]) &lt; n). Any hints what is going wrong here?</p>