SVM by OpenOpt.QP

example9_2.py

#SVM by OpenOpt
#use 'classification.txt'
from scipy import *
from scipy.linalg import norm
import numpy as np
from openopt import QP
from pylab import *

C = 0.5 # 
SIGMA = 0.45 # the parameter of gaussian_kernel

def gaussian_kernel(x, y):
    return exp(-norm(x-y)**2 / (2*(SIGMA**2)))

kernel = gaussian_kernel

def func(x, a, t, b, X, S):
    return b + sum((a[i] * t[i] * kernel(x, X[i]) for i in arange(len(S))))

def example9():
    data = np.loadtxt('classification.txt')
    X = data[:, 0:2] # training set
    t = data[:, 2] * 2 - 1.0 # label
    
    # make Gram matrix: K
    N = len(t)
    K = array([[t[i]*t[j]*kernel(X[i],X[j]) for j in range(N)] for i in range(N)])
    
    qp = QP(H = K,
             f = -ones(N),
             lb = zeros(N),
             ub = ones(N)*C,
             Aeq=t,
             beq=0)
    sol = qp.solve('nlp:ralg')
    a = sol.xf # Lagrange multiplier
    
    S = [i for i in range(N) if 0 < a[i]]
    M = [i for i in range(N) if 0 < a[i] < C]
    
    b = 0.0
    for i in M:
        b += t[i]
        for j in S:
            b -= a[j] * t[j] * kernel(X[i], X[j])
    b = b / len(M)
    
    # plot
    for i in range(len(X)):
        c = 'r' if t[i] > 0 else 'b'
        scatter([X[i][0]], [X[i][1]], color=c)
    
    X1, X2 = meshgrid(np.linspace(-2.5, 2.5, 50), np.linspace(-2.5, 2.5, 50))
    w, h = X1.shape
    X1.resize(X1.size)
    X2.resize(X2.size)
    Z = array([func(array([x1, x2]), a, t, b, X, S) for (x1, x2) in zip(X1, X2)])
    X1.resize((w,h))
    X2.resize((w,h))
    Z.resize((w,h))
    
    CS = contour(X1, X2, Z, [0.0],
                     colors = 'k',
                     linewidths = 3,
                     origin='lower')
        
    xlim(-2.5, 2.5)
    ylim(-2.5, 2.5)
    show()
    
if __name__ == '__main__':
    example9()