import numpy as np from numpy import ndarray from scipy.optimize import Bounds, LinearConstraint, minimize def norm_squared(vector: ndarray) -> float: """ Return the squared second norm of vector norm_squared(v) = sum(x * x for x in v) Args: vector (ndarray): input vector Returns: float: squared second norm of vector >>> norm_squared([1, 2]) 5 >>> norm_squared(np.asarray([1, 2])) 5 >>> norm_squared([0, 0]) 0 """ return np.dot(vector, vector) class SVC: """ Support Vector Classifier Args: kernel (str): kernel to use. Default: linear Possible choices: - linear regularization: constraint for soft margin (data not linearly separable) Default: unbound >>> SVC(kernel="asdf") Traceback (most recent call last): ... ValueError: Unknown kernel: asdf >>> SVC(kernel="rbf") Traceback (most recent call last): ... ValueError: rbf kernel requires gamma >>> SVC(kernel="rbf", gamma=-1) Traceback (most recent call last): ... ValueError: gamma must be > 0 """ def __init__( self, *, regularization: float = np.inf, kernel: str = "linear", gamma: float = 0.0, ) -> None: self.regularization = regularization self.gamma = gamma if kernel == "linear": self.kernel = self.__linear elif kernel == "rbf": if self.gamma == 0: raise ValueError("rbf kernel requires gamma") if not isinstance(self.gamma, (float, int)): raise ValueError("gamma must be float or int") if not self.gamma > 0: raise ValueError("gamma must be > 0") self.kernel = self.__rbf # in the future, there could be a default value like in sklearn # sklear: def_gamma = 1/(n_features * X.var()) (wiki) # previously it was 1/(n_features) else: msg = f"Unknown kernel: {kernel}" raise ValueError(msg) # kernels def __linear(self, vector1: ndarray, vector2: ndarray) -> float: """Linear kernel (as if no kernel used at all)""" return np.dot(vector1, vector2) def __rbf(self, vector1: ndarray, vector2: ndarray) -> float: """ RBF: Radial Basis Function Kernel Note: for more information see: https://en.wikipedia.org/wiki/Radial_basis_function_kernel Args: vector1 (ndarray): first vector vector2 (ndarray): second vector) Returns: float: exp(-(gamma * norm_squared(vector1 - vector2))) """ return np.exp(-(self.gamma * norm_squared(vector1 - vector2))) def fit(self, observations: list[ndarray], classes: ndarray) -> None: """ Fits the SVC with a set of observations. Args: observations (list[ndarray]): list of observations classes (ndarray): classification of each observation (in {1, -1}) """ self.observations = observations self.classes = classes # using Wolfe's Dual to calculate w. # Primal problem: minimize 1/2*norm_squared(w) # constraint: yn(w . xn + b) >= 1 # # With l a vector # Dual problem: maximize sum_n(ln) - # 1/2 * sum_n(sum_m(ln*lm*yn*ym*xn . xm)) # constraint: self.C >= ln >= 0 # and sum_n(ln*yn) = 0 # Then we get w using w = sum_n(ln*yn*xn) # At the end we can get b ~= mean(yn - w . xn) # # Since we use kernels, we only need l_star to calculate b # and to classify observations (n,) = np.shape(classes) def to_minimize(candidate: ndarray) -> float: """ Opposite of the function to maximize Args: candidate (ndarray): candidate array to test Return: float: Wolfe's Dual result to minimize """ s = 0 (n,) = np.shape(candidate) for i in range(n): for j in range(n): s += ( candidate[i] * candidate[j] * classes[i] * classes[j] * self.kernel(observations[i], observations[j]) ) return 1 / 2 * s - sum(candidate) ly_contraint = LinearConstraint(classes, 0, 0) l_bounds = Bounds(0, self.regularization) l_star = minimize( to_minimize, np.ones(n), bounds=l_bounds, constraints=[ly_contraint] ).x self.optimum = l_star # calculating mean offset of separation plane to points s = 0 for i in range(n): for j in range(n): s += classes[i] - classes[i] * self.optimum[i] * self.kernel( observations[i], observations[j] ) self.offset = s / n def predict(self, observation: ndarray) -> int: """ Get the expected class of an observation Args: observation (Vector): observation Returns: int {1, -1}: expected class >>> xs = [ ... np.asarray([0, 1]), np.asarray([0, 2]), ... np.asarray([1, 1]), np.asarray([1, 2]) ... ] >>> y = np.asarray([1, 1, -1, -1]) >>> s = SVC() >>> s.fit(xs, y) >>> s.predict(np.asarray([0, 1])) 1 >>> s.predict(np.asarray([1, 1])) -1 >>> s.predict(np.asarray([2, 2])) -1 """ s = sum( self.optimum[n] * self.classes[n] * self.kernel(self.observations[n], observation) for n in range(len(self.classes)) ) return 1 if s + self.offset >= 0 else -1 if __name__ == "__main__": import doctest doctest.testmod()