2

半正定プログラミングを解決するためのパッケージcvxoptに基づいて、Python で半正定埋め込みアルゴリズム (こちらを参照)を実装しようとしています。半正定値プログラムの定義を cvxopt のインターフェースにマッピングする際にいくつか問題があります (これを参照してください)。

これは私の現在の実装です:

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 or Rank([G; A]) < n)がスローされます。ここで何がうまくいかないのかヒントはありますか？

4

2 に答える 2

0

このエラーは、 にランク要件を課すソルバーから発生しますA。がまばらすぎると、cvxopt のすべてのソルバーでこの問題が発生しましたA。場合によっては、等式制約を取得し、A のゼロが低レベルのノイズで満たされてA*x=yいる同じ等式のノイズ バージョンを追加することで、結果を取得することに成功しました。Ap*x=y設定することでこれを行うことができます

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 に答える