Accompanying notebook for the recipe:
Christian Bauckhage: "NumPy / SciPy Recipes for Data Science: Kernel Least Squares Optimization (1)",
Technical Report, March 2015 (Download from ResearchGate)
Abstract of the paper: In this note, we show that least squares optimization is amenable to the kernel trick. This provides great flexibility in model fitting and we consider examples that illustrate this. In particular, we will discuss how kernel functions can implicitly introduce non-linearity into the least squares method.
import numpy as np
import numpy.linalg as la
import numpy.random as rnd
import scipy.spatial as spt
import scipy.sparse.linalg as sla
import matplotlib.pyplot as plt
def create_data_2D(n, w, xmin=-4, xmax=8):
x = rnd.random(n) * (xmax - xmin) + xmin
X = np.vander(x, len(w))
y = np.dot(X, w) + rnd.randn(n) * 0.5
return x, y
n = 50
wTrain = np.array([0.1, -0.8, 0.0, 11.5])
xTrain, yTrain = create_data_2D(n, wTrain)
plt.scatter(xTrain, yTrain)
plt.show()
XTrain = np.vander(xTrain, 2)
w = sla.lsmr(XTrain, yTrain, damp=1.)[0]
N = 2 * n
xTest = np.linspace(-4, 8, N)
XTest = np.vander(xTest, 2)
yTest = XTest.dot(w)
plt.scatter(xTrain, yTrain)
plt.plot(xTest, yTest, '-', color='orange')
plt.show()
def polyKernelMat(X, p):
return (X.dot(X.T) + 1.)**p
def polyKernelVec(x, X, p):
return (x.dot(X.T) + 1.)**p
p = 3
K = polyKernelMat(XTrain, p)
KI = la.inv(K + 1. * np.identity(n))
KIy = KI.dot(yTrain)
xTest = np.linspace(-4, 8, N)
XTest = np.vander(xTest, 2)
yTest = np.zeros(N)
for i in range(N):
k = polyKernelVec(XTest[i,:], XTrain, p)
yTest[i] = k.dot(KIy)
plt.scatter(xTrain, yTrain)
plt.plot(xTest, yTest, '-', color='orange')
plt.show()
def squaredEDM(X):
V = spt.distance.pdist(X, 'sqeuclidean')
D = spt.distance.squareform(V)
return D
def gaussKernelMat(X, s):
D = squaredEDM(X)
K = np.exp(-0.5/s**2 * D)
return K
def gaussKernelVec(x, X, s):
d = np.sum((X-x)**2, axis=1)
k = np.exp(-0.5/s**2 * d)
return k
s = 2.5
K = gaussKernelMat(XTrain, s)
KI = la.inv(K + 1. * np.identity(n))
KIy = KI.dot(yTrain)
xTest = np.linspace(-4, 8, N)
XTest = np.vander(xTest, 2)
yTest = np.zeros(N)
for i in range(N):
k = gaussKernelVec(XTest[i,:], XTrain, s)
yTest[i] = k.dot(KIy)
plt.scatter(xTrain, yTrain)
plt.plot(xTest, yTest, '-', color='orange')
plt.show()
© C. Bauckhage and O. Cremers Licensed under a CC BY-NC 4.0 . |
Acknowledgments: This material was prepared within the project P3ML which is funded by the Ministry of Education and Research of Germany (BMBF) under grant number 01/S17064. The authors gratefully acknowledge this support. |