2

我正在尝试基于cvxopt包在 python 中实现半定嵌入算法（参见此处），以解决半定编程。我在将半定程序的定义映射到 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) < p 或 Rank([G; A]) < n)。任何提示这里出了什么问题？

4

2 回答 2

0

此错误来自求解器，它对A. A当太稀疏时，我在 cvxopt 的所有求解器中都遇到了这个问题。A*x=y在某些情况下，我通过采用等式约束并在其上附加相同等式的噪声版本成功地获得了结果，Ap*x=y其中 A 的零点充满了低级噪声。可以通过设置来做到这一点

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

并将它们传递给等式约束。您可能会发现此链接很有用： 使用 cvxopt 最小化核规范

于 2013-01-29T04:39:45.247 回答
0

这是一个旧线程，但这可能会帮助人们在未来寻找 MVU 的 python 实现：

https://github.com/buquicchiol/MaximumVarianceUnfolding

于 2018-11-14T14:03:46.280 回答