Machine Learning Primer -- Basics

Claudius Gros, WS 2024/25

Institut für theoretische Physik
Goethe-University Frankfurt a.M.

SVM - Soft Margins & Kernel Trick

non-separable training data

Lagrange parameters $ \ \ a_k\ $ diverge
when data is not linearly seperable
classification condition
$\qquad\quad \mathbf{w}\cdot\mathbf{x}_i-b- l_i \ge 0, \qquad\quad l_i=\pm1 $
cannot be enforced

soft contraint

relax classification contraint
slack variables $\ \ \xi_i\ge0$

$$ \begin{array}{rclcrcl} \mathbf{w}\cdot\mathbf{x}_i&\ge& b+1-\xi_i &\quad\mbox{for}\quad & y_i&=&+1 \\ \mathbf{w}\cdot\mathbf{x}_i&\le& b-1+\xi_i&\quad\mbox{for}\quad & y_i&=&-1 \\ \end{array} $$

individual reduction of margin $\ \ (1-\xi_i)/|\mathbf{w}|$
may become negative

optimization condition

$$ \mathrm{min}\left[ \frac{1}{2}\mathbf{w}^2 +C\,\sum_i \xi_i \right]\,,\qquad\quad l_i\big(\mathbf{w}\cdot\mathbf{x}_i-b\big)- 1+\xi_i\ge 0 $$

weighting of miss-classification/margin: $\ \ C$

dual form

$$ L_P = \frac{\mathbf{w}^2}{2}+C\,\sum_i\xi_i -\sum_i a_i \Big[l_i\big(\mathbf{w}\cdot\mathbf{x}_i-b\big)-1+\xi_i\Big] -\sum_ib_i\xi_i $$

$b_i>0\ \ $ Lagrange parameter for constraint $\ \ \xi_i>0$
$a_i>0\ \ $ in analogy
stationarity: weight vector

$$ \frac{\partial L_P}{\partial\mathbf{w}} = 0 \qquad\quad \Rightarrow \qquad\quad \fbox{$\phantom{\big|} \mathbf{w} = \sum_i a_i l_i\mathbf{x}_i \phantom{\big|}$} $$

stationarity: hyperplane offset

$$ \frac{\partial L_P}{\partial b} = 0 \qquad\quad \Rightarrow \qquad\quad \fbox{$\phantom{\big|} \sum_i a_i l_i=0 \phantom{\big|}$} $$

stationarity: slack variables

$$ \frac{\partial L_P}{\partial \xi_i} = 0 \qquad \Rightarrow \qquad C=a_i+b_i \qquad \Rightarrow \qquad \fbox{$\phantom{\big|} 0\le a_i\le C \phantom{\big|}$} $$

$C \ \ $ limits runaway growth of Lagrange parameters $\ \ a_i$

soft margins - numerics

$$ \underset{\{a_l\}}{\mathrm{max}}\left[ \sum_i a_i-\frac{1}{2}\sum_{ij} a_iH_{ij}a_j\right]\,, \qquad 0\le a_i \le C, \qquad \sum_i a_il_i=0 $$

solve dual form
weight vector: $\ \ \mathbf{w} = \sum_i a_i l_i\mathbf{x}_i$

generalized support vectors

$$ a_i \quad\to\quad \left\{\begin{array}{rcl} 0 & & \mathrm{correctly\ classified} \\ ]0,C[ & & \mathrm{support\ vector} \\ C & & \mathrm{miss-classified} \\ \end{array} \right. $$

number of support vectors: $\ \ N_s$
hyperplane offset

$$ b=\frac{1}{N_S}\sum_{s\in S}\big(\mathbf{w}\cdot\mathbf{x}_s-l_s\big) $$

averaging over suppoort vectors if convergence is poor

SVM with soft margins

need to select free parameter $\ \ C$
: typically for machine learning
$C\ $ large: extended convergence times
$C\ $ small: support vectors / margins don't match
LDA classifier removed
:: stand-alone code

#!/usr/bin/env python3

import math                        # math
import matplotlib.pyplot as plt    # plotting
import numpy as np
from numpy.linalg import inv       # inverse matrix
from numpy import linalg as LA     # linear algebra


# ***
# *** covariance matrix from data
# ***

def covarianceMatrix(data):
  '''normalized covariance matrix of input data'''
  nRow = len(data)                     # number of data points
  nCol = len(data[0])                  # dimension
  mean = [0.0 for iCol in range(nCol)]
  for thisRow in data:                 # loop over data points
    mean += thisRow*1.0/nRow           # vector-mean
#
  coVar = np.zeros((nCol,nCol))        # empty matrix
  for thisRow in data:
    dataPoint = thisRow - mean         # shifted data point
    coVar += np.outer(dataPoint,dataPoint)/nRow    # outer product
#
  if (1==2):
    print("\ninput data\n", data)
    print("nRow, nCol\n", nRow, nCol)
    print("mean\n", mean)
    print("\n covariance matrix: measured \n",coVar)
#
  return mean, coVar                   # returning both

# ***
# *** generate data for plotting an ellipse
# ***

def plotEllipseData(data):
  "generates ellipse from symmetric 2x2 slice of the input matrix"
#
  slice22 = [inner[:2] for inner in data[:2]]   # slice
  slice22[0][1] = slice22[1][0]                 # symmetrize
                                                # ^ should not be necessary
  eigenValues, eigenVectors = LA.eig(slice22)   # eigensystem
#
  if (eigenValues[0]<0.0) or (eigenValues[1]<0.0):
    print("# plotEllipseData: only positive eigenvalues (variance)")
    return
#
  if (1==2):
    print("\n# coVar     \n", coVar  )
    print("\n# slice22     ", slice22)
    print("\n# eigenValues ", eigenValues)
    print("\n# eigenVectors\n",eigenVectors)
#
  a = math.sqrt(eigenValues[0])
  b = math.sqrt(eigenValues[1])
  cTheta = eigenVectors[0][0]
  sTheta = eigenVectors[1][0]
  x = []
  y = []
  for i in range(nPoints:=101):          # walrus assignment
    tt = i*2.0*math.pi/(nPoints-1)       # full loop
    cc = math.cos(tt)
    ss = math.sin(tt)
    xx = a*cTheta*cc - b*sTheta*ss
    yy = a*sTheta*cc + b*cTheta*ss
    x.append(xx)
    y.append(yy)
#   print(xx,yy)
  return x, y

# ***
# *** generate 2D test data
# ***

def testData(angle, var1, var2, nData, startMean=[0.0,0.0]):
  "2D Gaussian for a given angle and main variances.\
   A = \sum_i \lambda_i |lambda_i><lambda_i|"
#
  eigen1 = [math.cos(angle),-math.sin(angle)]
  eigen2 = [math.sin(angle), math.cos(angle)]
  startCoVar  = var1*np.outer(eigen1,eigen1)
  startCoVar += var2*np.outer(eigen2,eigen2)
#  print("\n covariance matrix: data generation \n",startCoVar)
  return np.random.multivariate_normal(startMean, startCoVar, nData)

# *******
# *** SVM 
# *******

def SVM_ascent(svmData, svmL, svmSoft):
  "Simple SVM code"
  nIter        =  10000         # fixed number of update iterations
  epsilon      =  0.001         # update rate
  svmThreshold =  0.001         # for support vectors
  nData = len(svmL)             # total number of labeled data points
  svmA  = np.random.rand(nData) # random [0,1] Lagrange parameters
  svmW  = np.zeros(2)           # 2D w-vector

  for iIter in range(nIter):                 # loop
    svmW = [0.0, 0.0]
    for ii in range(nData):   
      svmW += svmL[ii]*svmA[ii]*svmData[ii]
#   print(f'# {iIter:5d} {svmW[0]:6.2f} {svmW[1]:6.2f}')

    for ii in range(nData):                  # updating L-parameters
      svmA[ii] += epsilon*(1.0-svmL[ii]*np.dot(svmData[ii],svmW))

    factor = np.dot(svmA,svmL)*1.0/nData
    for ii in range(nData):                  # orthogonalization
      svmA[ii] = svmA[ii] - factor*svmL[ii]

    for ii in range(nData): 
      svmA[ii] = max(0.0,svmA[ii])           # positiveness
      svmA[ii] = min(svmSoft,svmA[ii])       # soft bound
#
# --- iteration loop finished
# --- threshold per support vector
#
  svmB  = np.zeros(nData) 
  BB_mean   = 0.0
  BB_number = 0
  for ii in range(nData):                    # positiveness
    if (svmA[ii]>svmThreshold):
      svmB[ii] = np.dot(svmData[ii],svmW) - svmL[ii]
      BB_mean   += svmB[ii]
      BB_number += 1
  BB_mean = BB_mean/BB_number

  print("# SVM data ")
# print(f'# {svmW[0]:6.2f} {svmW[1]:6.2f}')
  for ii in range(nData):                  
#   if (svmA[ii]>svmThreshold):
      print(f'{ii:5d} {svmA[ii]:10.6f} {svmL[ii]:5.1f} {svmB[ii]:10.4f} ',
            end="")
      print(f'{svmData[ii][0]:6.2f} {svmData[ii][1]:6.2f} ')
#
  return svmA, BB_mean, svmW     # Lagrange / threshold / w-vector
 

# ********
# *** main
# ********

dataMatrix_A  = testData( 0.2*math.pi, 1.0, 9.0, 20, [ 4.0,0.0])
dataMatrix_B  = testData(-0.3*math.pi, 1.0,12.0, 20, [-4.0,0.0])
mean_A, coVar_A = covarianceMatrix(dataMatrix_A)
mean_B, coVar_B = covarianceMatrix(dataMatrix_B)

#
# --- SVM preparation
#
dataLabel_A = [ 1.0 for _ in range(len(dataMatrix_A))]
dataLabel_B = [-1.0 for _ in range(len(dataMatrix_B))]
data_SVM   = np.vstack((dataMatrix_A,dataMatrix_B))    # vertical stacking
data_label = np.hstack((dataLabel_A,dataLabel_B))      # horizontal stacking
                                                       # do SVM
svmSoft = 15.0                                         # soft margin threshold
data_lagrange, data_BB, data_ww = SVM_ascent(data_SVM, data_label, svmSoft)

if (1==2):
  print()
  print(dataMatrix_A)
  print(dataLabel_A)
  print(dataMatrix_B)
  print(dataLabel_B)
  print(data_SVM)
  print(data_label)

#
# --- data points (normal, support)
#

xA_normal  = []
yA_normal  = []
xB_normal  = []
yB_normal  = []

xA_support = []
yA_support = []
xB_support = []
yB_support = []

xA_miss    = []
yA_miss    = []
xB_miss    = []
yB_miss    = []

svmThreshold = 0.001
for ii in range(len(data_label)):
  if (data_label[ii]>0) and (data_lagrange[ii]<svmThreshold):
   xA_normal.append(data_SVM[ii][0])
   yA_normal.append(data_SVM[ii][1])

  if (data_label[ii]<0) and (data_lagrange[ii]<svmThreshold):
   xB_normal.append(data_SVM[ii][0])
   yB_normal.append(data_SVM[ii][1])

  if (data_label[ii]>0) and (data_lagrange[ii]>svmThreshold):
   if (data_lagrange[ii]>0.95*svmSoft):
     xA_miss.append(data_SVM[ii][0])
     yA_miss.append(data_SVM[ii][1])
   else:
     xA_support.append(data_SVM[ii][0])
     yA_support.append(data_SVM[ii][1])

  if (data_label[ii]<0) and (data_lagrange[ii]>svmThreshold):
   if (data_lagrange[ii]>0.95*svmSoft):
     xB_miss.append(data_SVM[ii][0])
     yB_miss.append(data_SVM[ii][1])
   else:
     xB_support.append(data_SVM[ii][0])
     yB_support.append(data_SVM[ii][1])

print("# number of support vectors ",len(xA_support),len(xB_support))
print("# miss-classfied    vectors ",len(xA_miss),len(xB_miss))


#
# --- SVM plane and margins
#

svmPlane_x = [0.0, 0.0]
svmPlane_y = [0.0, 0.0]
svmMar_A_x = [0.0, 0.0]
svmMar_A_y = [0.0, 0.0]
svmMar_B_x = [0.0, 0.0]
svmMar_B_y = [0.0, 0.0]
r  = math.sqrt(data_ww[0]*data_ww[0]+data_ww[1]*data_ww[1])
rr = r*r
Len = 4.0
svmPlane_x[0] = (data_BB-0.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmPlane_y[0] = (data_BB-0.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmPlane_x[1] = (data_BB-0.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmPlane_y[1] = (data_BB-0.0)*data_ww[1]/rr + Len*data_ww[0]/r

svmMar_A_x[0] = (data_BB-1.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmMar_A_y[0] = (data_BB-1.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmMar_A_x[1] = (data_BB-1.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmMar_A_y[1] = (data_BB-1.0)*data_ww[1]/rr + Len*data_ww[0]/r

svmMar_B_x[0] = (data_BB+1.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmMar_B_y[0] = (data_BB+1.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmMar_B_x[1] = (data_BB+1.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmMar_B_y[1] = (data_BB+1.0)*data_ww[1]/rr + Len*data_ww[0]/r

#
# --- printing, including covariance ellipse
#

xEllipse_A, yEllipse_A = plotEllipseData(coVar_A)
xEllipse_B, yEllipse_B = plotEllipseData(coVar_B)  

Z_95 = math.sqrt(5.991)                  
xE_95_A = [Z_95*xx + mean_A[0] for xx in xEllipse_A]
yE_95_A = [Z_95*yy + mean_A[1] for yy in yEllipse_A]
xE_95_B = [Z_95*xx + mean_B[0] for xx in xEllipse_B]
yE_95_B = [Z_95*yy + mean_B[1] for yy in yEllipse_B]

x_connectMean = [ mean_A[0], mean_B[0] ]
y_connectMean = [ mean_A[1], mean_B[1] ]

plt.plot(xE_95_A, yE_95_A, "k", linewidth=2.0, label="95%")
plt.plot(xE_95_B, yE_95_B, "b", linewidth=2.0, label="95%")

plt.plot(xA_normal, yA_normal, "ok", markersize=8)
plt.plot(xB_normal, yB_normal, "ob", markersize=8)
plt.plot(xA_support, yA_support, "ok", markersize=10, 
         markeredgewidth=2, mfc='none')
plt.plot(xB_support, yB_support, "ob", markersize=10,
         markeredgewidth=2, mfc='none')
plt.plot(xA_miss, yA_miss, "ok", markersize=8, markeredgewidth=2,
         mfc=[1.0,0.84,0.0], label="miss")
plt.plot(xB_miss, yB_miss, "ob", markersize=8, markeredgewidth=2,
         mfc=[1.0,0.84,0.0], label="miss")

plt.plot(x_connectMean, y_connectMean, "r", linewidth=3.0)
plt.plot(x_connectMean, y_connectMean, "or", markersize=9.0)

plt.plot(svmPlane_x, svmPlane_y, color=[1.0,0.84,0.0], 
         linewidth=4.0, label="SVM plane")
plt.plot(svmMar_A_x, svmMar_A_y, color=[1.0,0.84,0.0], 
         linewidth=4.0, linestyle="dashed", label="SVM margin")
plt.plot(svmMar_B_x, svmMar_B_y, color=[1.0,0.84,0.0],
         linewidth=4.0, linestyle="dashed", label="")


plt.legend(loc="upper left")
#plt.axis('square')                           # square plot
plt.savefig('foo.svg')                       # export figure
plt.show()

#!/usr/bin/env python3

import math                        # math
import matplotlib.pyplot as plt    # plotting
import numpy as np
from numpy.linalg import inv       # inverse matrix
from numpy import linalg as LA     # linear algebra

# ***
# *** covariance matrix from data
# ***

def covarianceMatrix(data):
  '''normalized covariance matrix of input data'''
  nRow = len(data)                     # number of data points
  nCol = len(data[0])                  # dimension
  mean = [0.0 for iCol in range(nCol)]
  for thisRow in data:                 # loop over data points
    mean += thisRow*1.0/nRow           # vector-mean
#
  coVar = np.zeros((nCol,nCol))        # empty matrix
  for thisRow in data:
    dataPoint = thisRow - mean         # shifted data point
    coVar += np.outer(dataPoint,dataPoint)/nRow    # outer product
#
  if (1==2):
    print("\ninput data\n", data)
    print("nRow, nCol\n", nRow, nCol)
    print("mean\n", mean)
    print("\n covariance matrix: measured \n",coVar)
#
  return mean, coVar                   # returning both

# ***
# *** generate data for plotting an ellipse
# ***

def plotEllipseData(data):
  "generates ellipse from symmetric 2x2 slice of the input matrix"
#
  slice22 = [inner[:2] for inner in data[:2]]   # slice
  slice22[0][1] = slice22[1][0]                 # symmetrize
                                                # ^ should not be necessary
  eigenValues, eigenVectors = LA.eig(slice22)   # eigensystem
#
  if (eigenValues[0]<0.0) or (eigenValues[1]<0.0):
    print("# plotEllipseData: only positive eigenvalues (variance)")
    return
#
  if (1==2):
    print("\n# coVar     \n", coVar  )
    print("\n# slice22     ", slice22)
    print("\n# eigenValues ", eigenValues)
    print("\n# eigenVectors\n",eigenVectors)
#
  a = math.sqrt(eigenValues[0])
  b = math.sqrt(eigenValues[1])
  cTheta = eigenVectors[0][0]
  sTheta = eigenVectors[1][0]
  x = []
  y = []
  for i in range(nPoints:=101):          # walrus assignment
    tt = i*2.0*math.pi/(nPoints-1)       # full loop
    cc = math.cos(tt)
    ss = math.sin(tt)
    xx = a*cTheta*cc - b*sTheta*ss
    yy = a*sTheta*cc + b*cTheta*ss
    x.append(xx)
    y.append(yy)
#   print(xx,yy)
  return x, y

# ***
# *** generate 2D test data
# ***

def testData(angle, var1, var2, nData, startMean=[0.0,0.0]):
  "2D Gaussian for a given angle and main variances.\
   A = \sum_i \lambda_i |lambda_i><lambda_i|"
#
  eigen1 = [math.cos(angle),-math.sin(angle)]
  eigen2 = [math.sin(angle), math.cos(angle)]
  startCoVar  = var1*np.outer(eigen1,eigen1)
  startCoVar += var2*np.outer(eigen2,eigen2)
#  print("\n covariance matrix: data generation \n",startCoVar)
  return np.random.multivariate_normal(startMean, startCoVar, nData)

# *******
# *** SVM 
# *******

def SVM_ascent(svmData, svmL, svmSoft):
  "Simple SVM code"
  nIter        =  10000         # fixed number of update iterations
  epsilon      =  0.001         # update rate
  svmThreshold =  0.001         # for support vectors
  nData = len(svmL)             # total number of labeled data points
  svmA  = np.random.rand(nData) # random [0,1] Lagrange parameters
  svmW  = np.zeros(2)           # 2D w-vector

for iIter in range(nIter):                 # loop
    svmW = [0.0, 0.0]
    for ii in range(nData):   
      svmW += svmL[ii]*svmA[ii]*svmData[ii]
#   print(f'# {iIter:5d} {svmW[0]:6.2f} {svmW[1]:6.2f}')

for ii in range(nData):                  # updating L-parameters
      svmA[ii] += epsilon*(1.0-svmL[ii]*np.dot(svmData[ii],svmW))

factor = np.dot(svmA,svmL)*1.0/nData
    for ii in range(nData):                  # orthogonalization
      svmA[ii] = svmA[ii] - factor*svmL[ii]

for ii in range(nData): 
      svmA[ii] = max(0.0,svmA[ii])           # positiveness
      svmA[ii] = min(svmSoft,svmA[ii])       # soft bound
#
# --- iteration loop finished
# --- threshold per support vector
#
  svmB  = np.zeros(nData) 
  BB_mean   = 0.0
  BB_number = 0
  for ii in range(nData):                    # positiveness
    if (svmA[ii]>svmThreshold):
      svmB[ii] = np.dot(svmData[ii],svmW) - svmL[ii]
      BB_mean   += svmB[ii]
      BB_number += 1
  BB_mean = BB_mean/BB_number

print("# SVM data ")
# print(f'# {svmW[0]:6.2f} {svmW[1]:6.2f}')
  for ii in range(nData):                  
#   if (svmA[ii]>svmThreshold):
      print(f'{ii:5d} {svmA[ii]:10.6f} {svmL[ii]:5.1f} {svmB[ii]:10.4f} ',
            end="")
      print(f'{svmData[ii][0]:6.2f} {svmData[ii][1]:6.2f} ')
#
  return svmA, BB_mean, svmW     # Lagrange / threshold / w-vector

# ********
# *** main
# ********

dataMatrix_A  = testData( 0.2*math.pi, 1.0, 9.0, 20, [ 4.0,0.0])
dataMatrix_B  = testData(-0.3*math.pi, 1.0,12.0, 20, [-4.0,0.0])
mean_A, coVar_A = covarianceMatrix(dataMatrix_A)
mean_B, coVar_B = covarianceMatrix(dataMatrix_B)

#
# --- SVM preparation
#
dataLabel_A = [ 1.0 for _ in range(len(dataMatrix_A))]
dataLabel_B = [-1.0 for _ in range(len(dataMatrix_B))]
data_SVM   = np.vstack((dataMatrix_A,dataMatrix_B))    # vertical stacking
data_label = np.hstack((dataLabel_A,dataLabel_B))      # horizontal stacking
                                                       # do SVM
svmSoft = 15.0                                         # soft margin threshold
data_lagrange, data_BB, data_ww = SVM_ascent(data_SVM, data_label, svmSoft)

if (1==2):
  print()
  print(dataMatrix_A)
  print(dataLabel_A)
  print(dataMatrix_B)
  print(dataLabel_B)
  print(data_SVM)
  print(data_label)

#
# --- data points (normal, support)
#

xA_normal  = []
yA_normal  = []
xB_normal  = []
yB_normal  = []

xA_support = []
yA_support = []
xB_support = []
yB_support = []

xA_miss    = []
yA_miss    = []
xB_miss    = []
yB_miss    = []

svmThreshold = 0.001
for ii in range(len(data_label)):
  if (data_label[ii]>0) and (data_lagrange[ii]<svmThreshold):
   xA_normal.append(data_SVM[ii][0])
   yA_normal.append(data_SVM[ii][1])

if (data_label[ii]<0) and (data_lagrange[ii]<svmThreshold):
   xB_normal.append(data_SVM[ii][0])
   yB_normal.append(data_SVM[ii][1])

if (data_label[ii]>0) and (data_lagrange[ii]>svmThreshold):
   if (data_lagrange[ii]>0.95*svmSoft):
     xA_miss.append(data_SVM[ii][0])
     yA_miss.append(data_SVM[ii][1])
   else:
     xA_support.append(data_SVM[ii][0])
     yA_support.append(data_SVM[ii][1])

if (data_label[ii]<0) and (data_lagrange[ii]>svmThreshold):
   if (data_lagrange[ii]>0.95*svmSoft):
     xB_miss.append(data_SVM[ii][0])
     yB_miss.append(data_SVM[ii][1])
   else:
     xB_support.append(data_SVM[ii][0])
     yB_support.append(data_SVM[ii][1])

print("# number of support vectors ",len(xA_support),len(xB_support))
print("# miss-classfied    vectors ",len(xA_miss),len(xB_miss))

#
# --- SVM plane and margins
#

svmPlane_x = [0.0, 0.0]
svmPlane_y = [0.0, 0.0]
svmMar_A_x = [0.0, 0.0]
svmMar_A_y = [0.0, 0.0]
svmMar_B_x = [0.0, 0.0]
svmMar_B_y = [0.0, 0.0]
r  = math.sqrt(data_ww[0]*data_ww[0]+data_ww[1]*data_ww[1])
rr = r*r
Len = 4.0
svmPlane_x[0] = (data_BB-0.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmPlane_y[0] = (data_BB-0.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmPlane_x[1] = (data_BB-0.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmPlane_y[1] = (data_BB-0.0)*data_ww[1]/rr + Len*data_ww[0]/r

svmMar_A_x[0] = (data_BB-1.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmMar_A_y[0] = (data_BB-1.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmMar_A_x[1] = (data_BB-1.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmMar_A_y[1] = (data_BB-1.0)*data_ww[1]/rr + Len*data_ww[0]/r

svmMar_B_x[0] = (data_BB+1.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmMar_B_y[0] = (data_BB+1.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmMar_B_x[1] = (data_BB+1.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmMar_B_y[1] = (data_BB+1.0)*data_ww[1]/rr + Len*data_ww[0]/r

#
# --- printing, including covariance ellipse
#

xEllipse_A, yEllipse_A = plotEllipseData(coVar_A)
xEllipse_B, yEllipse_B = plotEllipseData(coVar_B)

Z_95 = math.sqrt(5.991)                  
xE_95_A = [Z_95*xx + mean_A[0] for xx in xEllipse_A]
yE_95_A = [Z_95*yy + mean_A[1] for yy in yEllipse_A]
xE_95_B = [Z_95*xx + mean_B[0] for xx in xEllipse_B]
yE_95_B = [Z_95*yy + mean_B[1] for yy in yEllipse_B]

x_connectMean = [ mean_A[0], mean_B[0] ]
y_connectMean = [ mean_A[1], mean_B[1] ]

plt.plot(xE_95_A, yE_95_A, "k", linewidth=2.0, label="95%")
plt.plot(xE_95_B, yE_95_B, "b", linewidth=2.0, label="95%")

plt.plot(xA_normal, yA_normal, "ok", markersize=8)
plt.plot(xB_normal, yB_normal, "ob", markersize=8)
plt.plot(xA_support, yA_support, "ok", markersize=10, 
         markeredgewidth=2, mfc='none')
plt.plot(xB_support, yB_support, "ob", markersize=10,
         markeredgewidth=2, mfc='none')
plt.plot(xA_miss, yA_miss, "ok", markersize=8, markeredgewidth=2,
         mfc=[1.0,0.84,0.0], label="miss")
plt.plot(xB_miss, yB_miss, "ob", markersize=8, markeredgewidth=2,
         mfc=[1.0,0.84,0.0], label="miss")

plt.plot(x_connectMean, y_connectMean, "r", linewidth=3.0)
plt.plot(x_connectMean, y_connectMean, "or", markersize=9.0)

plt.plot(svmPlane_x, svmPlane_y, color=[1.0,0.84,0.0], 
         linewidth=4.0, label="SVM plane")
plt.plot(svmMar_A_x, svmMar_A_y, color=[1.0,0.84,0.0], 
         linewidth=4.0, linestyle="dashed", label="SVM margin")
plt.plot(svmMar_B_x, svmMar_B_y, color=[1.0,0.84,0.0],
         linewidth=4.0, linestyle="dashed", label="")

plt.legend(loc="upper left")
#plt.axis('square')                           # square plot
plt.savefig('foo.svg')                       # export figure
plt.show()

SVM with Pytorch

first objective functions: minimize
:: $\ \mathbf{w}^2/2$
second objective functions: enforce margin condition
$z_i\le0,\qquad z_i = 1-l_i(\mathbf{w}\cdot\mathbf{x}_i-b)$
$\Rightarrow\ $ minimize $\ \max\big[0,z_i\big]$
$(1-C)$ / $C$ relative weights of objective functions
ReLU units: $\ \ \sigma(z)=\max(0,z)$.
direct minimization with Pytorch computational graph
:: backward pass evaluates gradients for
:: tensors with requires_grad=True

#!/usr/bin/env python3

import torch                       # pytorch; ML
import math                        # math
import matplotlib.pyplot as plt    # plotting
import numpy as np
from numpy.linalg import inv       # inverse matrix
from numpy import linalg as LA     # linear algebra


# ***
# *** covariance matrix from data
# ***

def covarianceMatrix(data):
  '''normalized covariance matrix of input data'''
  nRow = len(data)                     # number of data points
  nCol = len(data[0])                  # dimension
  mean = [0.0 for iCol in range(nCol)]
  for thisRow in data:                 # loop over data points
    mean += thisRow*1.0/nRow           # vector-mean
#
  coVar = np.zeros((nCol,nCol))        # empty matrix
  for thisRow in data:
    dataPoint = thisRow - mean         # shifted data point
    coVar += np.outer(dataPoint,dataPoint)/nRow    # outer product
#
  if (1==2):
    print("\ninput data\n", data)
    print("nRow, nCol\n", nRow, nCol)
    print("mean\n", mean)
    print("\n covariance matrix: measured \n",coVar)
#
  return mean, coVar                   # returning both

# ***
# *** generate data for plotting an ellipse
# ***

def plotEllipseData(data):
  "generates ellipse from symmetric 2x2 slice of the input matrix"
#
  slice22 = [inner[:2] for inner in data[:2]]   # slice
  slice22[0][1] = slice22[1][0]                 # symmetrize
                                                # ^ should not be necessary
  eigenValues, eigenVectors = LA.eig(slice22)   # eigensystem
#
  if (eigenValues[0]<0.0) or (eigenValues[1]<0.0):
    print("# plotEllipseData: only positive eigenvalues (variance)")
    return
#
  if (1==2):
    print("\n# coVar     \n", coVar  )
    print("\n# slice22     ", slice22)
    print("\n# eigenValues ", eigenValues)
    print("\n# eigenVectors\n",eigenVectors)
#
  a = math.sqrt(eigenValues[0])
  b = math.sqrt(eigenValues[1])
  cTheta = eigenVectors[0][0]
  sTheta = eigenVectors[1][0]
  x = []
  y = []
  for i in range(nPoints:=101):          # walrus assignment
    tt = i*2.0*math.pi/(nPoints-1)       # full loop
    cc = math.cos(tt)
    ss = math.sin(tt)
    xx = a*cTheta*cc - b*sTheta*ss
    yy = a*sTheta*cc + b*cTheta*ss
    x.append(xx)
    y.append(yy)
#   print(xx,yy)
  return x, y

# ***
# *** generate 2D test data
# ***

def testData(angle, var1, var2, nData, startMean=[0.0,0.0]):
  "2D Gaussian for a given angle and main variances.\
   A = \sum_i \lambda_i |lambda_i><lambda_i|"
#
  eigen1 = [math.cos(angle),-math.sin(angle)]
  eigen2 = [math.sin(angle), math.cos(angle)]
  startCoVar  = var1*np.outer(eigen1,eigen1)
  startCoVar += var2*np.outer(eigen2,eigen2)
#  print("\n covariance matrix: data generation \n",startCoVar)
  return np.random.multivariate_normal(startMean, startCoVar, nData)

# ******************
# *** SVM with torch
# ******************

def SVM_torch(Data, Label, C):
  """SVM; direct minimization of loss function, 
   Loss = (1-C)*|w|/2 
        + C ∑ max[0, 1 - y ( x*w - b ) ]^2"""
#
  nIter        = 10000          # fixed number of update iterations
  epsilon      =  0.01          # update rate

# adaptable variables  W, B
# on most machines double == float64
  W = torch.randn(2,requires_grad=True,dtype=torch.double)
  B = torch.randn(1,requires_grad=True,dtype=torch.double)

  for iIter in range(nIter):

# forward pass, with element-wise
# addition, multiplication, power, ReLU(x) = max(0,x)
    loss = (1.0-C)*0.5*torch.dot(W,W) \
         + C*torch.relu(1.0-Label*(torch.matmul(Data,W)-B[0])).pow(2).sum()

    loss.backward()                          # backward pass

    with torch.no_grad():                    # off computational graph
      W -= epsilon*W.grad                    # steepest descent
      B -= epsilon*B.grad                      
      W.grad = None                          # remove gradients
      B.grad = None
#
    if (iIter%100==0):
      print(f'{iIter:5d} {W[0]:8.4f} {W[1]:8.4f} {B[0]:8.4f} {0.0:8.4f} ')
#     print(f'{iIter:4d} {W.grad[0]:8.4f} {W.grad[1]:8.4f} {B.grad[0]:8.4f} ')
#
  return W.tolist(), B.tolist()              # returning non-tensors

# ********
# *** main
# ********

listData_A  = testData( 0.2*math.pi, 1.0, 9.0, 10, [ 4.0,0.0])
listData_B  = testData(-0.3*math.pi, 1.0,12.0, 10, [-4.0,0.0])
mean_A, coVar_A = covarianceMatrix(listData_A)
mean_B, coVar_B = covarianceMatrix(listData_B)

x_Ellipse_A, y_Ellipse_A = plotEllipseData(coVar_A)
x_Ellipse_B, y_Ellipse_B = plotEllipseData(coVar_B)  

#
# --- combine lists; lists to torch tensors
#
listLabel_A = [ 1.0 for _ in range(len(listData_A))]
listLabel_B = [-1.0 for _ in range(len(listData_B))]
list_allData  = np.vstack((listData_A,listData_B)) 
list_allLabel = np.hstack((listLabel_A,listLabel_B))  

tensor_allData  = torch.tensor(list_allData)     # fixed tensors
tensor_allLabel = torch.tensor(list_allLabel)    # non adaptable

torch_C = 0.9                                    # in [0,1]
data_ww, data_BB = SVM_torch(tensor_allData, tensor_allLabel, torch_C)


if (1==1):
  print("\n# list_allLabel  \n",list_allLabel)
  print("\n# data_ww        \n",data_ww      )
  print("\n# data_BB        \n",data_BB      )

#
# --- SVM plane and margins
#

svmPlane_x = [0.0, 0.0]
svmPlane_y = [0.0, 0.0]
svmMar_A_x = [0.0, 0.0]
svmMar_A_y = [0.0, 0.0]
svmMar_B_x = [0.0, 0.0]
svmMar_B_y = [0.0, 0.0]
r  = math.sqrt(data_ww[0]*data_ww[0]+data_ww[1]*data_ww[1])
rr = r*r
Len = 6.0
svmPlane_x[0] = (data_BB[0]-0.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmPlane_y[0] = (data_BB[0]-0.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmPlane_x[1] = (data_BB[0]-0.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmPlane_y[1] = (data_BB[0]-0.0)*data_ww[1]/rr + Len*data_ww[0]/r

svmMar_A_x[0] = (data_BB[0]-1.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmMar_A_y[0] = (data_BB[0]-1.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmMar_A_x[1] = (data_BB[0]-1.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmMar_A_y[1] = (data_BB[0]-1.0)*data_ww[1]/rr + Len*data_ww[0]/r

svmMar_B_x[0] = (data_BB[0]+1.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmMar_B_y[0] = (data_BB[0]+1.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmMar_B_x[1] = (data_BB[0]+1.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmMar_B_y[1] = (data_BB[0]+1.0)*data_ww[1]/rr + Len*data_ww[0]/r

#
# --- printing, including covariance ellipse
#

Z_95 = math.sqrt(5.991)                  
xE_95_A = [Z_95*xx + mean_A[0] for xx in x_Ellipse_A]
yE_95_A = [Z_95*yy + mean_A[1] for yy in y_Ellipse_A]
xE_95_B = [Z_95*xx + mean_B[0] for xx in x_Ellipse_B]
yE_95_B = [Z_95*yy + mean_B[1] for yy in y_Ellipse_B]

x_connectMean = [ mean_A[0], mean_B[0] ]
y_connectMean = [ mean_A[1], mean_B[1] ]

plt.plot(xE_95_A, yE_95_A, "k", linewidth=2.0, label="95%")
plt.plot(xE_95_B, yE_95_B, "b", linewidth=2.0, label="95%")

plt.plot(listData_A[:,0], listData_A[:,1], "ok", markersize=8)
plt.plot(listData_B[:,0], listData_B[:,1], "ob", markersize=8)

plt.plot(x_connectMean, y_connectMean, "r", linewidth=3.0)
plt.plot(x_connectMean, y_connectMean, "or", markersize=9.0)

plt.plot(svmPlane_x, svmPlane_y, color=[1.0,0.84,0.0], 
         linewidth=4.0, label="SVM plane")
plt.plot(svmMar_A_x, svmMar_A_y, color=[1.0,0.84,0.0], 
         linewidth=4.0, linestyle="dashed", label="SVM margin")
plt.plot(svmMar_B_x, svmMar_B_y, color=[1.0,0.84,0.0],
         linewidth=4.0, linestyle="dashed", label="")

plt.legend(loc="upper left")
#plt.axis('square')                           # square plot
plt.savefig('foo.svg')                       # export figure
plt.show()

#!/usr/bin/env python3

import torch                       # pytorch; ML
import math                        # math
import matplotlib.pyplot as plt    # plotting
import numpy as np
from numpy.linalg import inv       # inverse matrix
from numpy import linalg as LA     # linear algebra

# ***
# *** covariance matrix from data
# ***

def covarianceMatrix(data):
  '''normalized covariance matrix of input data'''
  nRow = len(data)                     # number of data points
  nCol = len(data[0])                  # dimension
  mean = [0.0 for iCol in range(nCol)]
  for thisRow in data:                 # loop over data points
    mean += thisRow*1.0/nRow           # vector-mean
#
  coVar = np.zeros((nCol,nCol))        # empty matrix
  for thisRow in data:
    dataPoint = thisRow - mean         # shifted data point
    coVar += np.outer(dataPoint,dataPoint)/nRow    # outer product
#
  if (1==2):
    print("\ninput data\n", data)
    print("nRow, nCol\n", nRow, nCol)
    print("mean\n", mean)
    print("\n covariance matrix: measured \n",coVar)
#
  return mean, coVar                   # returning both

# ***
# *** generate data for plotting an ellipse
# ***

def plotEllipseData(data):
  "generates ellipse from symmetric 2x2 slice of the input matrix"
#
  slice22 = [inner[:2] for inner in data[:2]]   # slice
  slice22[0][1] = slice22[1][0]                 # symmetrize
                                                # ^ should not be necessary
  eigenValues, eigenVectors = LA.eig(slice22)   # eigensystem
#
  if (eigenValues[0]<0.0) or (eigenValues[1]<0.0):
    print("# plotEllipseData: only positive eigenvalues (variance)")
    return
#
  if (1==2):
    print("\n# coVar     \n", coVar  )
    print("\n# slice22     ", slice22)
    print("\n# eigenValues ", eigenValues)
    print("\n# eigenVectors\n",eigenVectors)
#
  a = math.sqrt(eigenValues[0])
  b = math.sqrt(eigenValues[1])
  cTheta = eigenVectors[0][0]
  sTheta = eigenVectors[1][0]
  x = []
  y = []
  for i in range(nPoints:=101):          # walrus assignment
    tt = i*2.0*math.pi/(nPoints-1)       # full loop
    cc = math.cos(tt)
    ss = math.sin(tt)
    xx = a*cTheta*cc - b*sTheta*ss
    yy = a*sTheta*cc + b*cTheta*ss
    x.append(xx)
    y.append(yy)
#   print(xx,yy)
  return x, y

# ***
# *** generate 2D test data
# ***

def testData(angle, var1, var2, nData, startMean=[0.0,0.0]):
  "2D Gaussian for a given angle and main variances.\
   A = \sum_i \lambda_i |lambda_i><lambda_i|"
#
  eigen1 = [math.cos(angle),-math.sin(angle)]
  eigen2 = [math.sin(angle), math.cos(angle)]
  startCoVar  = var1*np.outer(eigen1,eigen1)
  startCoVar += var2*np.outer(eigen2,eigen2)
#  print("\n covariance matrix: data generation \n",startCoVar)
  return np.random.multivariate_normal(startMean, startCoVar, nData)

# ******************
# *** SVM with torch
# ******************

def SVM_torch(Data, Label, C):
  """SVM; direct minimization of loss function, 
   Loss = (1-C)*|w|/2 
        + C ∑ max[0, 1 - y ( x*w - b ) ]^2"""
#
  nIter        = 10000          # fixed number of update iterations
  epsilon      =  0.01          # update rate

# adaptable variables  W, B
# on most machines double == float64
  W = torch.randn(2,requires_grad=True,dtype=torch.double)
  B = torch.randn(1,requires_grad=True,dtype=torch.double)

for iIter in range(nIter):

# forward pass, with element-wise
# addition, multiplication, power, ReLU(x) = max(0,x)
    loss = (1.0-C)*0.5*torch.dot(W,W) \
         + C*torch.relu(1.0-Label*(torch.matmul(Data,W)-B[0])).pow(2).sum()

loss.backward()                          # backward pass

with torch.no_grad():                    # off computational graph
      W -= epsilon*W.grad                    # steepest descent
      B -= epsilon*B.grad                      
      W.grad = None                          # remove gradients
      B.grad = None
#
    if (iIter%100==0):
      print(f'{iIter:5d} {W[0]:8.4f} {W[1]:8.4f} {B[0]:8.4f} {0.0:8.4f} ')
#     print(f'{iIter:4d} {W.grad[0]:8.4f} {W.grad[1]:8.4f} {B.grad[0]:8.4f} ')
#
  return W.tolist(), B.tolist()              # returning non-tensors

# ********
# *** main
# ********

listData_A  = testData( 0.2*math.pi, 1.0, 9.0, 10, [ 4.0,0.0])
listData_B  = testData(-0.3*math.pi, 1.0,12.0, 10, [-4.0,0.0])
mean_A, coVar_A = covarianceMatrix(listData_A)
mean_B, coVar_B = covarianceMatrix(listData_B)

x_Ellipse_A, y_Ellipse_A = plotEllipseData(coVar_A)
x_Ellipse_B, y_Ellipse_B = plotEllipseData(coVar_B)

#
# --- combine lists; lists to torch tensors
#
listLabel_A = [ 1.0 for _ in range(len(listData_A))]
listLabel_B = [-1.0 for _ in range(len(listData_B))]
list_allData  = np.vstack((listData_A,listData_B)) 
list_allLabel = np.hstack((listLabel_A,listLabel_B))

tensor_allData  = torch.tensor(list_allData)     # fixed tensors
tensor_allLabel = torch.tensor(list_allLabel)    # non adaptable

torch_C = 0.9                                    # in [0,1]
data_ww, data_BB = SVM_torch(tensor_allData, tensor_allLabel, torch_C)

if (1==1):
  print("\n# list_allLabel  \n",list_allLabel)
  print("\n# data_ww        \n",data_ww      )
  print("\n# data_BB        \n",data_BB      )

#
# --- SVM plane and margins
#

svmPlane_x = [0.0, 0.0]
svmPlane_y = [0.0, 0.0]
svmMar_A_x = [0.0, 0.0]
svmMar_A_y = [0.0, 0.0]
svmMar_B_x = [0.0, 0.0]
svmMar_B_y = [0.0, 0.0]
r  = math.sqrt(data_ww[0]*data_ww[0]+data_ww[1]*data_ww[1])
rr = r*r
Len = 6.0
svmPlane_x[0] = (data_BB[0]-0.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmPlane_y[0] = (data_BB[0]-0.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmPlane_x[1] = (data_BB[0]-0.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmPlane_y[1] = (data_BB[0]-0.0)*data_ww[1]/rr + Len*data_ww[0]/r

svmMar_A_x[0] = (data_BB[0]-1.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmMar_A_y[0] = (data_BB[0]-1.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmMar_A_x[1] = (data_BB[0]-1.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmMar_A_y[1] = (data_BB[0]-1.0)*data_ww[1]/rr + Len*data_ww[0]/r

svmMar_B_x[0] = (data_BB[0]+1.0)*data_ww[0]/rr + Len*data_ww[1]/r
svmMar_B_y[0] = (data_BB[0]+1.0)*data_ww[1]/rr - Len*data_ww[0]/r
svmMar_B_x[1] = (data_BB[0]+1.0)*data_ww[0]/rr - Len*data_ww[1]/r
svmMar_B_y[1] = (data_BB[0]+1.0)*data_ww[1]/rr + Len*data_ww[0]/r

#
# --- printing, including covariance ellipse
#

Z_95 = math.sqrt(5.991)                  
xE_95_A = [Z_95*xx + mean_A[0] for xx in x_Ellipse_A]
yE_95_A = [Z_95*yy + mean_A[1] for yy in y_Ellipse_A]
xE_95_B = [Z_95*xx + mean_B[0] for xx in x_Ellipse_B]
yE_95_B = [Z_95*yy + mean_B[1] for yy in y_Ellipse_B]

x_connectMean = [ mean_A[0], mean_B[0] ]
y_connectMean = [ mean_A[1], mean_B[1] ]

plt.plot(xE_95_A, yE_95_A, "k", linewidth=2.0, label="95%")
plt.plot(xE_95_B, yE_95_B, "b", linewidth=2.0, label="95%")

plt.plot(listData_A[:,0], listData_A[:,1], "ok", markersize=8)
plt.plot(listData_B[:,0], listData_B[:,1], "ob", markersize=8)

plt.plot(x_connectMean, y_connectMean, "r", linewidth=3.0)
plt.plot(x_connectMean, y_connectMean, "or", markersize=9.0)

plt.plot(svmPlane_x, svmPlane_y, color=[1.0,0.84,0.0], 
         linewidth=4.0, label="SVM plane")
plt.plot(svmMar_A_x, svmMar_A_y, color=[1.0,0.84,0.0], 
         linewidth=4.0, linestyle="dashed", label="SVM margin")
plt.plot(svmMar_B_x, svmMar_B_y, color=[1.0,0.84,0.0],
         linewidth=4.0, linestyle="dashed", label="")

plt.legend(loc="upper left")
#plt.axis('square')                           # square plot
plt.savefig('foo.svg')                       # export figure
plt.show()

SVM regression

linear regression

data points $\ \ \big(\mathbf{x}_i,y_i\big)$
linear function $\ \ \mathbf{w}\cdot\mathbf{x}_i-b, \qquad R^d\,\to\, R$

$\qquad\quad \mathrm{min}\,\sum_i \big[y_i-(\mathbf{w}\cdot\mathbf{x}_i-b)\big]^2 $

least square fit

soft SVM

slack variables $\ \ \xi_i,\, \xi_i^*\ge 0$
conditions

$$ y_i-\mathbf{w}\cdot\mathbf{x}_i+b \le \epsilon+\xi_i, \qquad\quad \mathbf{w}\cdot\mathbf{x}_i-b-y_i \le \epsilon+\xi_i^* $$

miminize

$$ \frac{1}{2}\mathbf{w}^2 + C\,\sum_i\big(\xi_i+\xi_i^*\big) $$

Lagrange parameters $\ \ \Rightarrow \ \ $ dual form
: literature

SVM regression involves a confidence intervall $\ \ 2\epsilon$

non-linear classification

overfitting

given enough parameters,
one can fit an elephant

fitting $ \Leftrightarrow $ generalization
dilemma of ML
however:
double descent
SVM parsimonious
: kernel trick

non-linear transformation

linear classification on transformed data
$\qquad\quad \mathbf{x}\ \to \mathbf{z},\qquad R^d \ \to \ R^D $
high-dimensional $\ \ D$
: brain, neural nets, echo-state machine, $\dots$
examples

$$ \begin{array}{rcl c rcl c rcl} z_1 &=& x_1 &\quad& z_2&=& x_2 &\quad& z_3&=& x_1^2 \\ z_4 &=& \cos(x_7-x_8) &\quad& z_5&=& \tanh(x_5 x_6^2) &\quad& z_6&=& \sin(3 x_4^3) \\ \end{array} $$

$\dots$ explosively large $\ \ D$

classification

parameters $\ \ \mathbf{w}_D=\big(w_1,\,\dots,\, w_{D-1}\big)$
training data $ \ \ \mathbf{x}_i \ \to\ \mathbf{z}_i$

$$ \mathbf{w}_D\cdot\mathbf{z}_i-b \quad \to\quad \left\{ \begin{array}{rcl} >0 && \mbox{first class} \\ <0 && \mbox{second class} \end{array}\right. $$

high-dimenional dual

number of trainings data: $\ \ N$

$$ \underset{\mathbf{a}}{\mathrm{max}}\left[\, \sum_{i=0}^{N-1} a_i-\frac{1}{2}\sum_{i,j=0}^{N-1} a_i H_{ij} a_j \,\right]\,, \qquad\quad H_{ij} = l_i\big(\mathbf{z}_i\cdot\mathbf{z}_j\big)l_j $$

with slack variables

$$ \mathbf{w}_D=\sum_{j=0}^N a_jl_j\mathbf{z}_j, \qquad\quad 0\le a_i\le C, \quad\qquad \sum_{i=0}^{N-1} a_il_i=0 $$

classification variable

$$ \mathbf{w}_D\cdot\mathbf{z}_k-b = \sum_{j=0}^N a_jl_j\big(\mathbf{z}_j\cdot\mathbf{z}_k\big)-b $$

depends exclusively on the D-dimensional
scalar product $\ \ \mathbf{z}_j\cdot\mathbf{z}_k$
idem for offset $$ b=\mathbf{w}_D\cdot\mathbf{z}_s-l_s= \sum_{j=0}^N a_jl_j\big(\mathbf{z}_j\cdot\mathbf{z}_s\big)-l_s $$ for any $\ \ 0\le a_s\le S \ \ $ (support vector $\ \mathbf{z}_s$)

kernel trick

the large dimension D enters only implicitly
via the D-dimensional scalar products
$\ \ \mathbf{z}_i\cdot\mathbf{z}_j$

kernel $\ \ K(\mathbf{u},\mathbf{v})$

$$ \mathbf{z}_i\cdot\mathbf{z}_j\quad \to\quad K(\mathbf{x}_i,\mathbf{x}_j) $$

rewrite everything in terms of the kernel
Lagrangian

$$ \underset{\mathbf{a}}{\mathrm{max}}\left[\, \sum_{i=0}^{N-1} a_i-\frac{1}{2}\sum_{i,j=0}^{N-1} a_i l_i K(\mathbf{x}_i,\mathbf{x}_j) l_ja_j \,\right]\,, $$

classification argument $$ \mathbf{w}_D\cdot\mathbf{z}_k-b = \sum_{j=0}^N a_jl_j K(\mathbf{x}_j,\mathbf{x}_k) -b, \qquad\quad b= \sum_j a_jl_j K(\mathbf{x}_j,\mathbf{x}_s)-l_s $$ combined

$$ \mathbf{w}_D\cdot\mathbf{z}_k-b = \sum_{j=0}^N a_jl_j \Big[ K(\mathbf{x}_j,\mathbf{x}_k)-K(\mathbf{x}_j,\mathbf{x}_s) \Big] +l_s, $$

may be averaged over support vectors

kernel trick

select simple kernels

kernels

polynomial of degeree $\ \ \alpha$
mixed polynomial
Gaussian / radial basis functions (RBF)
sigmoidal (neuronal)

kernel optimization

kernel may depend on (many) parameters
tune kernel parameters using classification success as objective function
$\equiv$ meta-paramter optimization
: stochastic search / gradient descent, $\dots$
: ML methods

SVM kernel with Sklearn

sklearn

Scikit-learn

#!/usr/bin/env python3
# source
# https://scikit-learn.org/stable/auto_examples/exercises/plot_iris_exercise.html

import matplotlib.pyplot as plt
import numpy as np

from sklearn import datasets, svm

iris = datasets.load_iris()    # Iris (flower) dataset
X = iris.data
y = iris.target
#print(X)
#print(y)

# Iris has three classes; 0, 1, 2
# and four attributes;    0, 1, 2, 3    
X = X[y != 0, :2]   # attributes: 0, 1
y = y[y != 0]       # classes   : 1, 2

n_sample = len(X)

np.random.seed(0)
order = np.random.permutation(n_sample)
X = X[order]                             # shuffle
y = y[order].astype(float)               # dtype casting

X_train = X[: int(0.9 * n_sample)]       # first 90%
y_train = y[: int(0.9 * n_sample)]
X_test = X[int(0.9 * n_sample) :]        # last 10%
y_test = y[int(0.9 * n_sample) :]

# SVM analysis with several kernel types
for kernel in ("linear", "rbf", "poly"):
    clf = svm.SVC(kernel=kernel, gamma=10)  # SVM object
    clf.fit(X_train, y_train)

    plt.figure()
    plt.clf()
    plt.scatter(X[:, 0], X[:, 1], c=y, zorder=10, 
                cmap=plt.cm.Paired, edgecolor="k", s=20)

# Circle out the test data
    plt.scatter(X_test[:, 0], X_test[:, 1], s=80, 
                facecolors="none", zorder=10, edgecolor="k")

    plt.axis("tight")
    x_min = X[:, 0].min()
    x_max = X[:, 0].max()
    y_min = X[:, 1].min()
    y_max = X[:, 1].max()

    XX, YY = np.mgrid[x_min:x_max:200j, y_min:y_max:200j]
    Z = clf.decision_function(np.c_[XX.ravel(), YY.ravel()])

# Put the result into a color plot
    Z = Z.reshape(XX.shape)
    plt.pcolormesh(XX, YY, Z > 0, cmap=plt.cm.Paired)
    plt.contour(XX, YY, Z,
        colors=["k", "k", "k"],
        linestyles=["--", "-", "--"],
        levels=[-0.5, 0, 0.5])

    plt.title(kernel)
plt.show()

#!/usr/bin/env python3
# source
# https://scikit-learn.org/stable/auto_examples/exercises/plot_iris_exercise.html

import matplotlib.pyplot as plt
import numpy as np

from sklearn import datasets, svm

iris = datasets.load_iris()    # Iris (flower) dataset
X = iris.data
y = iris.target
#print(X)
#print(y)

# Iris has three classes; 0, 1, 2
# and four attributes;    0, 1, 2, 3    
X = X[y != 0, :2]   # attributes: 0, 1
y = y[y != 0]       # classes   : 1, 2

n_sample = len(X)

np.random.seed(0)
order = np.random.permutation(n_sample)
X = X[order]                             # shuffle
y = y[order].astype(float)               # dtype casting

X_train = X[: int(0.9 * n_sample)]       # first 90%
y_train = y[: int(0.9 * n_sample)]
X_test = X[int(0.9 * n_sample) :]        # last 10%
y_test = y[int(0.9 * n_sample) :]

# SVM analysis with several kernel types
for kernel in ("linear", "rbf", "poly"):
    clf = svm.SVC(kernel=kernel, gamma=10)  # SVM object
    clf.fit(X_train, y_train)

plt.figure()
    plt.clf()
    plt.scatter(X[:, 0], X[:, 1], c=y, zorder=10, 
                cmap=plt.cm.Paired, edgecolor="k", s=20)

# Circle out the test data
    plt.scatter(X_test[:, 0], X_test[:, 1], s=80, 
                facecolors="none", zorder=10, edgecolor="k")

plt.axis("tight")
    x_min = X[:, 0].min()
    x_max = X[:, 0].max()
    y_min = X[:, 1].min()
    y_max = X[:, 1].max()

XX, YY = np.mgrid[x_min:x_max:200j, y_min:y_max:200j]
    Z = clf.decision_function(np.c_[XX.ravel(), YY.ravel()])

# Put the result into a color plot
    Z = Z.reshape(XX.shape)
    plt.pcolormesh(XX, YY, Z > 0, cmap=plt.cm.Paired)
    plt.contour(XX, YY, Z,
        colors=["k", "k", "k"],
        linestyles=["--", "-", "--"],
        levels=[-0.5, 0, 0.5])

plt.title(kernel)
plt.show()

kernel, convolution and Green's function

kernel
later:
neural tangent kernels
convolution when $\ k(y,x)=k(y-x)$
later:
convolution nets

Green's function for the Poisson equation

density of walkers at $\ (x,t)$
diffusion constant $\ D$
starting state $\ \rho_0(y)=\rho(y,t\!=\!0))$
Green's function: propagation in time
proof via direct differentiation
diffusive transport $\ \langle x^2(t) \rangle = \sigma^2 = 2D\,t$
ballistic transport $\ \langle x(t) \rangle = v\,t$