python - 基于cvxopt的python中的半定嵌入

我试图基于cvxopt软件包在python中实现半定嵌入算法（请参见here），以解决半定编程问题。我在将半定程序的定义映射到cvxopt的接口时遇到一些问题（请参见this）。

这是我当前的实现：

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 > 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)

执行程序将引发ValueError（Rank（A）

最佳答案

此错误来自求解器，求解器对A进行排名要求。当A太稀疏时，所有cvxopt的求解器都遇到了这个问题。在某些情况下，我通过采用等式约束A*x=y并在其后附加一个等式的嘈杂版本Ap*x=y来成功获得结果，其中A的零点被低电平噪声填充。通过设置可以做到这一点

A1 = matrix([A,Ap])
y1 = matrix([y,y])

并通过平等约束。您可能会发现此链接很有用：minimizing nuclear norm with cvxopt

关于python - 基于cvxopt的python中的半定嵌入，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/13526606/