4 แบบฝึกหัด

class primal_SVM:
    def __init__(self):
        self.bopt = None
        self.wopt = None
        self.phi = None
            
    def train(self, x, y, wb0, 
                    options={'ftol': 1e-9, 'disp': True}):
        assert (x.shape[1] == y.shape[1]) and (y.shape[0] == 1) 

        N = x.shape[1]
        # projected x to z
        z = self.phi(x)
        assert z.shape[1] == x.shape[1]         
        D = z.shape[0] # number of projected dimensions
        
        if wb0 is None: # w: (D,1), b: scalar
            wb0 = np.random.normal(0, 1, D+1) 

        def primal_ineq(w, b):
            w = w.reshape((-1,1))
            b = np.asscalar(b)            
            cineq = y.reshape((1,-1)) * \
               (np.dot(w.T,z.reshape((D,-1))) + b) -1  
            return cineq.reshape((-1,)) # (M,) 
            
        def primal_dc(w, b):
            gradw = y.reshape((1,-1)) * z.reshape((D,-1))
            gradb = y.reshape((1,-1))                    
            dc = np.hstack( (gradw.T, gradb.T) )            
            return dc.reshape((-1,D+1)) # (M,D+1) 
        
        # inequality constraints: y[i] * ( w.T x[i] + b)-1 >= 0 
        ineq_cons = {'type': 'ineq', 
               'fun' : lambda wb: primal_ineq(wb[:-1], wb[-1]),     
               'jac' : lambda wb: primal_dc(wb[:-1], wb[-1])} 
        
        # Objective        
        def primal_f(wb): # primal loss: minimization form
            w = wb[:-1].reshape((-1,1)) # D x 1
            L = 0.5 * np.dot(w.T, w) # scalar
            return np.asscalar(L)

        def primal_grad(wb): # gradient of primal loss
            w = wb[:-1].reshape((-1,1)) # D x 1
            pgrad = np.vstack((w, 0))            
            return pgrad.reshape((-1,))   # (D,)
        
        minf = lambda wb: primal_f(wb)
        gradf =  lambda wb: primal_grad(wb)
        res = optimize.minimize(minf, wb0, method='SLSQP', 
          jac=gradf, constraints=ineq_cons, options=options)
        
        if not res.success:
            return res

        # Train succeeds.
        self.wopt = res.x[:-1].reshape((-1,1))
        self.bopt = res.x[-1]
        return res
    
    def decision_score(self, x):        
        # Project to feature space
        z = self.phi(x)
        assert z.shape[0] == self.wopt.shape[0]

        # Compute the score
        yhat = np.dot(self.wopt.T, z) + self.bopt
        return yhat        # 1 x N

class cSVM:
    def __init__(self):
        self.epsilon = 0.001
        self.bopt = None
        self.sv = None
        self.sy = None
        self.salpha = None
        self.kernel = None
            
    def train(self, x, y, C=1, a0=None,
              options={'ftol': 1e-9, 'disp': True}):
        assert (x.shape[1] == y.shape[1]) and (y.shape[0] == 1) 

        N = x.shape[1]
        if a0 is None:
            a0 = np.random.normal(0, 1, N)

        # parameter bounds: 0 <= alpha_i <= C for all i's
        bounds = optimize.Bounds([0 for i in range(N)], 
                                 [C for i in range(N)])
        
        # inequality constraints: dummy
        ineq_cons = {'type': 'ineq', 
                     'fun' : lambda a: np.array([0]),  # (1,)
                     'jac' : lambda a: np.zeros((1,N))}   # (M,N)
        
        # equality constraint: sum_i y_i alpha_i = 0
        eq_cons = {'type': 'eq', 'fun' : lambda a: np.dot(y, 
                         a.reshape((-1, 1))).reshape((-1,)),
           'jac' : lambda a: y.reshape((1,N))}                      
        
        # Objective
        K = self.kernel(x, x)     # N x N
        H = np.dot(y.T, y) * K    # N x N
        
        def dual_minf(a, H): # dual loss in minimization form
            a = a.reshape((-1,1))
            Q = 0.5 * np.dot(a.T, np.dot(H, a)) - np.sum(a) 
            return np.asscalar(Q)

        def dual_grad(a, H, N): # gradient of dual loss
            dQ = np.dot(H, a.reshape((-1,1))) - np.ones((N,1)) 
            return dQ.reshape((-1,))   # (N,)
        
        minf = lambda a: dual_minf(a, H)
        gradf =  lambda a: dual_grad(a, H, N)
        res = optimize.minimize(minf, a0, method='SLSQP', 
              jac=gradf, constraints=[ineq_cons, eq_cons], 
              bounds=bounds, options=options)
        
        if not res.success:
            return res
        
        alpha = res.x.reshape((-1,1))
        sv_ids = np.where(alpha > self.epsilon)[0]
        svp_ids = np.where( np.logical_and(alpha > self.epsilon, 
                                  alpha < C - self.epsilon) )[0]
        Ns = len(sv_ids)
        Np = len(svp_ids)
        
        if Np == 0: # No support vector on the edge
            print('# No support vector on the edge.')
            svp_ids = sv_ids
            Np = Ns
        
        # support vectors
        self.sv = x[:, sv_ids] # D x Ns
        self.sy = y[0, sv_ids].reshape((1,-1)) # 1 x Ns        
        spy = y[0, svp_ids].reshape((1,-1)) # 1 x Np
        
        # support alphas
        self.salpha = alpha[sv_ids]
        
        # optimal b
        sK = K[np.repeat(sv_ids, Np), 
               np.tile(svp_ids, Ns)].reshape((Ns, Np)) 
        self.bopt = np.mean(spy - \
                       np.dot(self.salpha.T * self.sy, sK))
        return res

    def decision_score(self, x):
        assert x.shape[0] == self.sv.shape[0] # x in D x N
        # Compute support kernels
        sK = self.kernel(x, self.sv) # N x Ns
        assert (sK.shape[0] == x.shape[1]) 
                      and (sK.shape[1] == self.sv.shape[1])

        yhat = np.dot(sK, self.salpha*self.sy.reshape((-1,1))) \
               + self.bopt
        return yhat.T        # 1 x N

def gaussian(xq, xp, sigma=1):
    assert xq.shape[0] == xp.shape[0]   # xq: D x Nx, xp: D x Ns 
    Nx = xq.shape[1]
    Ns = xp.shape[1]
    K = np.zeros((Nx, Ns))
    
    c = -1/(2*sigma**2)    
    for i in range(Nx):
        vi = xq[:,[i]] - xp # (D,1)-(D,Ns): broadcast to (D,Ns)
        K[i,:] = np.exp(c*np.sum(vi**2, axis=0))  # (Ns,)
    
    return K # Nx x Ns