In this blog post, you’ll implement the perceptron algorithm using numerical programming and demonstrate its use on synthetic data sets.
\[ \newcommand{\R}{\mathbb{R}} \newcommand{\vx}{\mathbf{x}} \newcommand{\vy}{\mathbf{y}} \newcommand{\mX}{\mathbf{X}} \newcommand{\vw}{\mathbf{w}} \newcommand{\bracket}[1]{\langle #1 \rangle} \newcommand{\paren}[1]{\left( #1 \right)} \]
In the perceptron algorithm, we have a data set described by:
We assume that the labels in our data can be (approximately) separated by a linear predictor. The strongest version of this assumption is that our data are linearly separable. This means that there exists a vector \(\vw \in p\) and a scalar \(b\) such that, for every \(i\),
\[ \begin{aligned} y_i &= \begin{cases} 0 & \langle \vw, \vx_i \rangle < b \\ 1 & \langle \vw, \vx_i \rangle \geq b \end{cases} \\ &= \mathbb{1}(\langle \vw, \vx_i \rangle \geq b) \end{aligned} \]
Here’s an example of what our data might look like. For visualization purposes, we are going to have \(p = 2\) features.
import numpy as np
import pandas as pd
import seaborn as sns
from matplotlib import pyplot as plt
from sklearn.datasets import make_blobs
np.random.seed(12345)
n = 100
p_features = 3
X, y = make_blobs(n_samples = 100, n_features = p_features - 1, centers = [(-1.7, -1.7), (1.7, 1.7)])
fig = plt.scatter(X[:,0], X[:,1], c = y)
xlab = plt.xlabel("Feature 1")
ylab = plt.ylabel("Feature 2")
You may be able to visualize a line that separates all the purple “0” labels from the yellow “1” labels. Our end-goal is a Python class that will allow us to find such a separating line, when it exists.
You are going to implement a Perceptron class in a script file called perceptron.py. (you should place it next to your .ipynb notebook in which you write your blog post.) You’ll import it like this:
from perceptron import PerceptronBy the time you have successfully implemented your Perceptron class, you should be able to replicate the following results on your blog post.
First, you should be able to import your class
from perceptron import PerceptronThen you should be able to create an instance of the class and fit it to data.
p = Perceptron()
p.fit(X, y, max_steps = 1000)This should result in p having an instance variable w of weights and an instance variable history of scores:
p.warray([2.10557404, 3.1165449 , 0.25079936])print(p.history[-10:]) #just the last few values[0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 1.0]You can visualize how the score evolved over time:
fig = plt.plot(p.history)
xlab = plt.xlabel("Iteration")
ylab = plt.ylabel("Accuracy")
In this particular example, the algorithm was able to achieve perfect, 100% classification.
You can also visualize the line that the algorithm finds to separate the data.
def draw_line(w, x_min, x_max):
x = np.linspace(x_min, x_max, 101)
y = -(w[0]*x + w[2])/w[1]
plt.plot(x, y, color = "black")
fig = plt.scatter(X[:,0], X[:,1], c = y)
fig = draw_line(p.w, -2, 2)
xlab = plt.xlabel("Feature 1")
ylab = plt.ylabel("Feature 2")
As we know from the previous investigations, the score on the data is:
p.score(X, y)1.0This blog posts as three main components: implement the perceptron algorithm, complete several experiments, and write careful prose to describe your findings.
Your class should have the following methods:
Perceptron.fit(X, y) is the primary method. This method has no return value. If p is a Perceptron, then after p.fit(X, y) is called, p should have an instance variable of weights called w. This w is the vector \(\tilde{\vw} = (\vw, -b)\) in the classifier above. Additionally, p should have an instance variable called p.history which is a list of the evolution of the score over the training period (see Perceptron.score(X, y) below.)Perceptron.predict(X) should return a vector \(\hat{\vy} \in \{0,1\}^n\) of predicted labels. These are the model’s predictions for the labels on the data.Perceptron.score(X, y) should return the accuracy of the perceptron as a number between 0 and 1, with 1 corresponding to perfect classification.Feel free to add any other methods or functions that you find helpful while implementing.
fit()To implement fit(), it’s convenient to consider a modified version of \(\mX\): \(\tilde{\mX} = [\mX, \mathbf{1}]\), where \(\mathbf{1} \in \R^n\) is a column-vector of \(1\)s. The reason this is handy is that if we also define \(\tilde{\vw} = (\vw, -b)\), then we can write our classification rule as
\[ \hat{y}_i = \mathbb{1}(\langle \tilde{\vw}, \tilde{\vx}_i\rangle \geq 0)\;. \]
This is mathematically convenient and makes it much easier for us to code up our algorithms.
With these definitions, the perceptron algorithm proceeds as follows:
In this expression, \(\tilde{y}_i = 2y_i - 1\) is a convenient version of \(y_i\) that takes values \(-1\) and \(1\) instead of \(0\) and \(1\).
This update is performed until either a user-specified maximum number of steps is reached or until the score (accuracy) reaches 1.0.
Note that in an iteration in which \(\tilde{y}_i \langle \tilde{\vw}^{(t)}, \tilde{\vx}_i\rangle \geq 0\), nothing happens. Take a moment to check that this occurs when the current weight vector \(\tilde{\vw}^{(t)}\) correctly classifies the tuple \((\vx_i, y_i)\).
You should be able to replicate the demo in Section 2 with your source code. Feel free to use that demo as a test case – your source code may be in good shape when you are able to fully replicate the results.
np.random.seed() immediately after importing your packages.An excellent solution will have exactly one for-loop, of the form:
for _ in range(max_steps):
# perform the perceptron update and log the score in self.historyThat is, you should not do any loops over the data! Use vectorized numpy operations and matrix-vector multiplication.
You should also not use if statements to perform comparisons between numbers.
For a hint on how you can avoid doing this, you can reflect on the following two code snippets:
print((1 < 2)*2)
print((1 > 2)*2)Please include informative docstrings for Perceptron.fit(), Perceptron.predict(), and Perceptron.score().
A concise solution should likely be no more than 60 lines of compact Python code (excluding comments and docstrings).
Please perform experiments (with visualizations) that illustrate the following:
In crafting your blog post, please include the following components:
perceptron.py) on your GitHub repo.What is the runtime complexity of a single iteration of the perceptron algorithm update as described by Equation 1? Assume that the relevant operations are addition and multiplication. Does the runtime complexity depend on the number of data points \(n\)? What about the number of features \(p\)?
You only need to consider this question in the context of a single update. The question of how many updates are required to converge is a trickier one that you don’t have to discuss in your blog post.
The code
X_ = np.append(X, np.ones((X.shape[0], 1)), 1)can be used to form the matrix \(\tilde{\mX} = [\mathbf{X}, \mathbf{1}]\).
© Phil Chodrow, 2023