Kernel Multivariate Analysis¶

KK = np.matrix([[np.dot(X[i,:],X[j,:]) for i in range(X.shape[0])]for j in range(X.shape[0])])
print(np.dot(X[0,:],X[0,:]), '\n', squareform(pdist(X, lambda x,y:np.dot(x,y)))) # must be equal but...

if not all(np.equal(squareform(pdist(X, lambda x,y:np.dot(x,y))), KK)):
  print('pdist provides zero elements for the diagnal elements of the output.')
if all(np.equal(cdist(X, X, lambda x,y:np.dot(x,y)), KK)):
  print('cdist is acceptable, but slow.')

# ラグランジュ乗数を二次計画法（Quadratic Programming）で求める
N = X.shape[0]
K = np.zeros((N,N))
K_x = squareform(pdist(X, lambda x,y:np.dot(x,y)))
for i in range(N):
    for j in range(N):
        K[i, j] = y[i] * y[j] * np.dot(X[i,:],X[j,:])
# L_dual = gamma * alpha - alpha.T * (y.T * K * y) * alpha -> max
# - L_dual -> min in solvers.qp
Q = cvxopt.matrix(K)
p = cvxopt.matrix(-np.ones(N))             # -1がN個の列ベクトル
G = cvxopt.matrix(np.diag([-1.0]*N))       # 対角成分が-1のNxN行列
h = cvxopt.matrix(np.zeros(N))             # 0がN個の列ベクトル
A = cvxopt.matrix(array(y.tolist()).astype('double'), (1,N))# N個の教師信号が要素の行ベクトル（1xN）
# Cast A to double type, thanks to a post on stack overflow https://stackoverflow.com/questions/36510859/cvxopt-qp-solver-typeerror-a-must-be-a-d-matrix-with-1000-columns
b = cvxopt.matrix(0.0)                     # 定数0.0
sol = cvxopt.solvers.qp(Q, p, G, h, A, b)  # 二次計画法でラグランジュ乗数aを求める
a = array(sol['x']).reshape(N)             # 'x'がaに対応する
print(a)

threshold = .00001

SV = np.where(a>=threshold)[0]
SV

import numpy as np
from scipy.linalg import norm
import cvxopt
import cvxopt.solvers
from pylab import *

N = 100  # データ数

def f(x1, w, b):
    return - (w[0] / w[1]) * x1 - (b / w[1])

def kernel(x, y):
    return np.dot(x, y)  # 線形カーネル


# 訓練データを作成
cls1 = []
cls2 = []

mean1 = [-1, 2]
mean2 = [1, -1]
cov = [[1.0,0.8], [0.8, 1.0]]
half_n = int(N/2)
cls1.extend(np.random.multivariate_normal(mean1, cov, half_n))
cls2.extend(np.random.multivariate_normal(mean2, cov, half_n))

# データ行列Xを作成
X = vstack((cls1, cls2))

# ラベルtを作成
t = []
for i in range(half_n):
    t.append(1.0)   # クラス1
for i in range(half_n):
    t.append(-1.0)  # クラス2
t = array(t)

# ラグランジュ乗数を二次計画法（Quadratic Programming）で求める
K = np.zeros((N, N))
for i in range(N):
    for j in range(N):
        K[i, j] = t[i] * t[j] * kernel(X[i], X[j])
# L_dual = gamma * alpha - alpha.T * (y.T * K * y) * alpha -> max
# - L_dual -> min in solvers.qp
Q = cvxopt.matrix(K)
p = cvxopt.matrix(-np.ones(N))             # -1がN個の列ベクトル
G = cvxopt.matrix(np.diag([-1.0]*N))       # 対角成分が-1のNxN行列
h = cvxopt.matrix(np.zeros(N))             # 0がN個の列ベクトル
A = cvxopt.matrix(t, (1,N))                # N個の教師信号が要素の行ベクトル（1xN）
b = cvxopt.matrix(0.0)                     # 定数0.0
sol = cvxopt.solvers.qp(Q, p, G, h, A, b)  # 二次計画法でラグランジュ乗数aを求める
a = array(sol['x']).reshape(N)             # 'x'がaに対応する
# print( a)



# サポートベクトルのインデックスを抽出
S = []
for i in range(len(a)):
    if a[i] < 0.00001: continue
    S.append(i)

# wを計算
w = np.zeros(2)
for n in S:
    w += a[n] * t[n] * X[n]

# bを計算
sum = 0
for n in S:
    temp = 0
    for m in S:
        temp += a[m] * t[m] * kernel(X[n], X[m])
    sum += (t[n] - temp)
b = sum / len(S)

print( S, b)

# 訓練データを描画
x1, x2 = np.array(cls1).transpose()
plot(x1, x2, 'rx')

x1, x2 = np.array(cls2).transpose()
plot(x1, x2, 'bx')

# サポートベクトルを描画
for n in S:
    scatter(X[n,0], X[n,1], s=80, c='c', marker='o')

# 識別境界を描画
x1 = np.linspace(-6, 6, 1000)
x2 = [f(x, w, b) for x in x1]
plot(x1, x2, 'g-')

xlim(-6, 6)
ylim(-6, 6)
show()

S