Scientific Machine Learning:
Machine Learning
Machine learning is the set of tools with which a model is fitted to data: a hypothesis class, a loss function that measures the misfit, and an optimizer that reduces it. This notebook builds those tools from the simplest case, linear regression, through classification with the perceptron, logistic regression and support vector machines, and ends with the nonlinear models that lead to neural networks. Every later notebook in this course, physics-informed or not, is an instance of the same three ingredients.
Learning Objectives¶
After studying this notebook, you should be able to
- set up linear regression as a least-squares problem, solve it by the normal equation, and identify the three ingredients — a model, a loss function, an optimization algorithm — that every later method is assembled from;
- state the perceptron algorithm, implement it, and explain both why it converges on separable data and what it cannot do;
- derive logistic regression from the idea of using every point's distance rather than only the misclassified ones, read its output as a probability, and connect the cross-entropy loss to maximum likelihood;
- explain what a deep learning framework provides — automatic differentiation, hardware acceleration, and reusable components — and why hand-derived gradients do not scale;
- write a training loop in PyTorch, using
requires_grad,backward(), an update undertorch.no_grad(), andgrad.zero_(), and use it to solve both linear and logistic regression; - explain why a built-in loss such as
BCEWithLogitsLossis preferred to composing a sigmoid with a logarithm by hand.
Regression is a fundamental concept in machine learning, playing a crucial role in understanding and predicting continuous outcomes based on input features. Unlike classification, which assigns data points to discrete categories, regression models aim to establish a relationship between independent variables (features) and a continuous dependent variable (target).
At its core, regression seeks to capture patterns in data and predict numerical values by minimizing the error between the predicted values and the actual values. This makes regression indispensable for a wide range of applications.
Why We Begin with Linear Regression
Linear regression is the simplest supervised learning problem, and it deserves careful study for a reason that goes well beyond its own usefulness. Almost every method in this course, including the neural networks introduced later, is assembled from the same three ingredients:
A model: a family of candidate functions $f(x\,;\theta)$ indexed by parameters $\theta$. The model expresses the shape of relationship we are willing to consider before seeing the data.
A loss function: a single number measuring how badly a particular choice of $\theta$ explains the observed data. Learning is meaningless until we have committed to a definition of "badly".
An optimization algorithm: a procedure for finding the $\theta$ that makes the loss as small as possible.
Linear regression is the one case in which all three ingredients are simple enough to be carried through completely by hand. The model is a straight line, the loss is a sum of squared errors, and the optimization has a closed-form solution obtained from linear algebra alone. Once this pattern is clear, the later methods can be read as substitutions into the same template: replace the straight line with a neural network, replace the squared error with cross-entropy, and replace the closed-form solution with gradient descent. The template itself never changes.
Roadmap
Sections 1.1 to 1.3 set up the problem, recast it as a least squares problem, and solve it using linear algebra. Section 1.4 solves the same problem with scikit-learn in a few lines of code. Section 1.5 lists the topics that a complete treatment of regression would include, which are left for independent study.
from IPython.display import YouTubeVideo
YouTubeVideo('IRrGVQV8vZ8', width = '560', height = '315', extras = ['referrerpolicy="strict-origin-when-cross-origin"'])
Let us begin by considering a simple linear regression problem, where the goal is to predict a continuous output based on a single input feature.
Throughout this section, keep a concrete situation in mind. Suppose $x_i$ is the applied load on a specimen and $y_i$ is the measured elongation, and that we have performed $m$ tests. We want a rule that predicts the elongation for a load we have not tested. Two things are worth noticing immediately about how this differs from the way such a problem is usually posed in a mechanics course. We are not deriving the relationship from Hooke's law, and we are not assuming the measurements are exact. We are asking what relationship the data itself supports.
Given a dataset of input-output pairs:
$$\; \begin{cases}x_{i} \; \text{: inputs} \\y_{i} \; \text{: outputs}\end{cases} \qquad \text{for} \;i = 1, 2, \cdots, m$$
We can organize the inputs and outputs as column vectors:
$$x= \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{m} \end{bmatrix}, \qquad y= \begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{m} \end{bmatrix}$$
Here $m$ denotes the number of data points, and the subscript $i$ indexes a single measurement. Each pair $(x_i, y_i)$ is one experiment, one observation, or one row of a data table. Stacking them into vectors is not merely notational tidiness. It is what allows the entire dataset to be handled by a single matrix operation later, which is the reason linear algebra is the natural language for this subject.
Prediction Model
Let $\hat{y}_i$ be the predicted output for input $x_i$. In general, we assume there's a function $f(x_i; \theta)$ that maps the input to the predicted output, parameterized by $\theta$:
$$ \hat{y}_{i} = f(x_{i}\,; \theta) $$
For linear regression, we use a linear model to approximate the relationship between $x_i$ and $y_i$:
$$y_i \approx \hat{y}_{i} = \theta_{0} + \theta_{1}x_{i}$$
Here,
- $\theta_0$ is the intercept (bias term),
- $\theta_1$ is the slope (weight) of the input, and
- $\theta = \begin{bmatrix} \theta_0 \\ \theta_1 \end{bmatrix}$ are the model parameters we want to learn.
The distinction between $x$ and $\theta$ is worth stating explicitly, because it is a common source of confusion. The data $x_i$ and $y_i$ are given and fixed; we cannot change them. The parameters $\theta_0$ and $\theta_1$ are unknown and adjustable, and "learning" means nothing more than choosing numerical values for them. Once $\theta$ is fixed, the model is a specific line and prediction is simple arithmetic. The entire difficulty of machine learning lies in choosing $\theta$ well.
Note also that $\hat{y}_i$ carries a hat and $y_i$ does not. This distinction is used consistently throughout the notebook: $y_i$ is what was measured, and $\hat{y}_i$ is what the model claims. They will never be exactly equal, and the gap between them is precisely what we are about to minimize.
Objective: Minimizing Prediction Error
To find the best values of $\theta_0$ and $\theta_1$, we minimize the sum of squared errors between the predicted values $\hat{y}_i$ and the actual values $y_i$:
$$ \min\limits_{\theta_{0}, \theta_{1}}\sum\limits_{i = 1}^{m} \left(\hat{y}_{i} - y_{i} \right)^2$$
This is known as the least squares criterion and is a standard method for fitting linear models.
Why Squared Error?
The choice of the squared error is so standard that it is easy to accept without asking why, but three separate reasons support it.
Errors must not cancel. If we simply summed $\hat{y}_i - y_i$, a large positive error and a large negative error would cancel, and a terrible model could report zero total error. Squaring makes every contribution non-negative.
The objective is smooth. The squared error is differentiable everywhere, so the minimum can be found by setting a derivative to zero. The absolute error $\lvert \hat{y}_i - y_i \rvert$ also prevents cancellation, but it has a corner at zero and does not admit a clean closed-form solution.
It has a probabilistic justification. If the measurement noise is assumed to be independent and Gaussian with constant variance, then minimizing the sum of squared errors is exactly equivalent to maximizing the likelihood of the observed data. We will meet this argument again in Section 3.1, where the same reasoning applied to binary labels produces the cross-entropy loss instead.
One consequence deserves a warning. Because the error is squared, a point that is twice as far from the line contributes four times as much to the loss. Least squares is therefore sensitive to outliers, and a single badly corrupted measurement can visibly tilt the fitted line.
Vector Form of the Model
We can express the model for all data points at once using vectors and matrix operations:
Read the expression below from right to left, as a statement about vectors rather than about individual points. The predicted output vector $\hat{y}$ is built from just two fixed vectors: a vector of ones and the vector of inputs. The parameters $\theta_0$ and $\theta_1$ are the amounts of each that we mix together.
This is a genuinely useful way to think about regression. Whatever values of $\theta_0$ and $\theta_1$ we choose, the resulting $\hat{y}$ is confined to the set of all linear combinations of those two vectors, that is, to a two-dimensional plane inside $m$-dimensional space. The measured vector $y$ almost certainly does not lie on that plane. Fitting the model therefore amounts to finding the point on the plane that is closest to $y$, and we will see this geometric picture again in Section 1.2.
$$ y= \begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{m} \end{bmatrix} \approx \begin{bmatrix} \hat y_{1} \\ \hat y_{2} \\ \vdots \\ \hat y_{m} \end{bmatrix} = \theta_{0} \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} + \theta_{1} \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{m} \end{bmatrix}$$
For convenience, we define a function that maps inputs to feature vectors, denoted as $\phi$
The purpose of this rewriting is not merely convenience, and it is worth pausing on. Introducing the feature map $\phi$ separates the model into two stages: a transformation of the input, followed by a linear combination. The transformation $\phi$ absorbs everything specific to the problem, while the linear combination stays the same regardless of what $\phi$ does.
The payoff appears as soon as we want a model that is not a straight line. Choosing $\phi(x_i) = \begin{bmatrix} 1 & x_i & x_i^2 \end{bmatrix}^T$ gives a quadratic fit, and choosing $\phi(x_i) = \begin{bmatrix} 1 & \sin x_i & \cos x_i \end{bmatrix}^T$ gives a sinusoidal fit, yet in both cases the solution procedure derived below is used without a single modification. The model is nonlinear in $x$, but it remains linear in $\theta$, and it is linearity in $\theta$, not in $x$, that makes least squares work. This idea returns in Chapter 2 under the name of feature representation, and again when kernels are discussed.
$$\begin{array}{lcr}\begin{align*} \hat{y}_{i} & = \theta_0 + x_i \theta_1 = 1 \cdot \theta_0 + x_i \theta_1 \\ \\ & = \begin{bmatrix}1 & x_{i}\end{bmatrix}\begin{bmatrix}\theta_{0} \\ \theta_{1}\end{bmatrix} \\\\ & =\begin{bmatrix}1 \\ x_{i} \end{bmatrix}^{T}\begin{bmatrix}\theta_{0} \\ \theta_{1}\end{bmatrix} \\\\ & =\phi^{T}(x_{i})\theta \end{align*}\end{array} \begin{array}{lcr} \quad \quad \text{feature vector} \; \phi(x_{i}) = \begin{bmatrix}1 \\ x_{i}\end{bmatrix} \end{array}$$
We can then construct the matrix $\Phi$ by stacking the transposed feature vectors for all $m$ data points:
$$\Phi = \begin{bmatrix}1 & x_{1} \\ 1 & x_{2} \\ \vdots \\1 & x_{m} \end{bmatrix}=\begin{bmatrix}\phi^T(x_{1}) \\\phi^T(x_{2}) \\\vdots \\\phi^T(x_{m}) \end{bmatrix} \quad \implies \quad \hat{y} = \begin{bmatrix}\hat{y}_{1} \\\hat{y}_{2} \\\vdots \\\hat{y}_{m}\end{bmatrix}=\Phi\theta$$
Optimization Problem (or Least Squares)
We aim to minimize the total squared error:
$$\min\limits_{\theta_{0}, \theta_{1}}\sum\limits_{i = 1}^{m} (\hat{y}_{i} - y_{i})^2 =\color{red}{\min\limits_{\theta}\lVert\Phi\theta-y\rVert^2_2} \qquad \qquad \left(\text{same as} \; \min_{x} \lVert Ax-b \rVert_2^2 \right)$$
The solution to this least squares problem is:
This expression is called the normal equation, and it can be derived in one line. Expanding the objective gives
$$\lVert A\theta - b \rVert_2^2 = \theta^T A^T A \theta - 2 b^T A \theta + b^T b .$$
Setting the gradient with respect to $\theta$ to zero yields $2 A^T A \theta - 2 A^T b = 0$, that is, $A^T A \theta = A^T b$, from which the formula follows provided $A^T A$ is invertible. Since the objective is a convex quadratic function of $\theta$, this stationary point is the global minimum, not merely a local one.
$$ \theta^* = (A^TA)^{-1}A^T b $$
Why? (here, $A_i$ means the $i$th row of matrix $A$)
$$\sum_{i=1}^{m} \lvert \underbrace{ A_i x - b_i}_{\text{scalar}} \rvert^2 = \left\| \begin{bmatrix} A_1 x - b_1 \\ \vdots \\ A_m x - b_m \end{bmatrix} \right\|^2 = \left\| \begin{bmatrix} A_1 \\ \vdots \\ A_m \end{bmatrix} x - \begin{bmatrix} b_1 \\ \vdots \\ b_m \end{bmatrix} \right\|^2 = \lVert \underbrace{Ax-b}_{\text{vector}} \rVert_2^2$$
The geometric interpretation of the least squares solution is illustrated below.
From Input to Prediction: Feature Mapping and Linear Model
$$\begin{array}{lcr} \text{input} \\ x_{i} \end{array} \quad \rightarrow \quad \begin{array}{lcr} \text{feature vector} \\ \begin{bmatrix} 1 \\ x_{i} \end{bmatrix} \end{array} \quad \rightarrow \quad \begin{array}{lcr} \text{predicted output} \\ \hat{y}_{i} \end{array}$$
Matrix Form of the Linear System
We collect all feature vectors into a design matrix and model the predictions as:
$$ \underbrace{ \begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_m \end{bmatrix} }_{A = \Phi} \underbrace{ \begin{bmatrix} \theta_0 \\ \theta_1 \end{bmatrix} }_{\theta} = \underbrace{ \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_m \end{bmatrix} }_{y} $$
Below each column of $A$, we can identify the column vectors (features):
$$ \begin{array}{cccccc} \vec{A}_1 = \begin{bmatrix}1 \\ 1 \\ \vdots \\ 1\end{bmatrix}, & \quad \vec{A}_2 = \begin{bmatrix}x_1 \\ x_2 \\ \vdots \\ x_m\end{bmatrix} , &\quad \quad \theta = \begin{bmatrix} \theta_0 \\ \theta_1 \end{bmatrix}, & \quad \quad y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_m \end{bmatrix} \end{array} $$
This system is typically over-determined (more equations than unknowns), so we seek the best approximate solution via least squares projection.
What "Over-determined" Means Geometrically
With $m$ measurements and only two unknowns, we have $m$ equations in two unknowns. For $m > 2$ such a system generally has no exact solution, and this is not a defect of the data. It is the normal situation, and it reflects the fact that real measurements contain noise while a straight line does not.
The column space picture explains what least squares does about it. The product $A\theta$ is a linear combination of the two columns $\vec{A}_1$ and $\vec{A}_2$, so as $\theta$ ranges over all possible values, $A\theta$ traces out the plane spanned by these two vectors inside $\mathbb{R}^m$. Solving $A\theta = y$ exactly would require $y$ to lie in that plane, which it does not.
Least squares therefore asks a weaker and answerable question: which point in the plane is closest to $y$? The answer is the orthogonal projection of $y$ onto the column space, and the residual $y - A\theta^*$ is perpendicular to every column of $A$. Writing that perpendicularity condition as $A^T(y - A\theta^*) = 0$ gives $A^T A \theta^* = A^T y$, which is the normal equation again. The algebraic derivation and the geometric derivation lead to the same place.
$$A(= \Phi) = \begin{bmatrix} \vec{A}_1 & \vec{A}_2 \end{bmatrix}$$
Solution via Linear Algebra
The linear regression problem can be solved analytically using linear algebra, in what is known as the least squares solution. The optimal solution is given by the normal equation:
$$ \theta = (A^TA)^{-1}A^T y $$
Three remarks about this formula are worth keeping in mind.
When does the inverse exist? The matrix $A^TA$ is invertible if and only if the columns of $A$ are linearly independent. For simple linear regression this fails only in a degenerate case, namely when every input $x_i$ takes the same value, in which case the data contains no information about the slope. With many features, however, the condition fails whenever one feature is an exact linear combination of the others, and it becomes numerically fragile when one feature is nearly a linear combination of the others.
Do not compute the inverse explicitly. Although the formula is written with $(A^TA)^{-1}$, forming the inverse is both slower and numerically less stable than solving the linear system directly. In practice one uses
numpy.linalg.lstsqornumpy.linalg.pinv, which rely on matrix factorizations rather than explicit inversion.This closed form is a luxury. Linear regression is one of the few models whose optimum can be written down in a single formula. For logistic regression in Chapter 3, and for every neural network afterwards, no such formula exists and the minimum must be found iteratively by gradient descent. It is worth appreciating the closed form here precisely because it is the last time we will have one.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# data points in column vector [input, output]
x = np.array([0.1, 0.4, 0.7, 1.2, 1.3, 1.7, 2.2, 2.8, 3.0, 4.0, 4.3, 4.4, 4.9]).reshape(-1, 1)
y = np.array([0.5, 0.9, 1.1, 1.5, 1.5, 2.0, 2.2, 2.8, 2.7, 3.0, 3.5, 3.7, 3.9]).reshape(-1, 1)
plt.figure(figsize = (6, 4))
plt.plot(x, y, 'ko', alpha = 0.3)
plt.xlabel('X')
plt.ylabel('Y')
plt.axis('equal')
plt.grid(alpha = 0.3)
plt.xlim([0, 5])
plt.show()
m = y.shape[0]
# A = np.hstack([np.ones([m, 1]), x])
A = np.hstack([x**0, x])
A = np.asmatrix(A)
theta = (A.T*A).I*A.T*y
print('theta:\n', theta)
# to plot
plt.figure(figsize = (6, 4))
plt.xlabel('X')
plt.ylabel('Y')
plt.plot(x, y, 'ko', label = 'data', alpha = 0.3)
# to plot a straight line (fitted line)
xp = np.arange(0, 5, 0.01).reshape(-1, 1)
yp = theta[0,0] + theta[1,0]*xp
plt.plot(xp, yp, 'r', linewidth = 2, label = 'regression')
plt.legend()
plt.axis('equal')
plt.grid(alpha = 0.3)
plt.xlim([0, 5])
plt.show()
Scikit-Learn is a powerful and user-friendly machine learning library in Python. It provides efficient implementations of many standard algorithms, including linear regression.
Having solved the problem by hand, we now solve it again with a library. The point of doing both is not redundancy. Working through the linear algebra shows what is being computed, while the library shows what practice actually looks like, and knowing only one of the two is a weak position. When a fitted model behaves strangely, the ability to reason about $A^TA$ is what allows the cause to be diagnosed rather than guessed at.
Scikit-learn also introduces an interface convention that recurs throughout the Python machine learning ecosystem. A model object is created, fit estimates its parameters from data, and predict applies the fitted model to new inputs. Nearly every scikit-learn model follows this same pattern, so learning it once for linear regression transfers directly to the perceptron and to logistic regression later in this notebook.
- Simple and efficient tools for data mining and data analysis
- Accessible to everybody, and reusable in various contexts
- Built on NumPy, SciPy, and Matplotlib
- Open source
- https://scikit-learn.org/stable/index.html#
from sklearn import linear_model
reg = linear_model.LinearRegression()
reg.fit(x, y)
reg.coef_
reg.intercept_
# to plot
plt.figure(figsize = (6, 4))
plt.xlabel('X')
plt.ylabel('Y')
plt.plot(x, y, 'ko', label = 'data', alpha = 0.3)
# to plot a straight line (fitted line)
plt.plot(xp, reg.predict(xp), 'r', linewidth = 2, label = 'regression')
plt.legend()
plt.axis('equal')
plt.grid(alpha = 0.3)
plt.xlim([0, 5])
plt.show()
The following topics are important for a comprehensive understanding of regression analysis but are not discussed in detail here. Readers are encouraged to study them independently.
Two of these deserve a comment before the list, because they will reappear in every later chapter. Overfitting and regularization are not peculiarities of regression. They are the central tension of all supervised learning: a model flexible enough to fit the training data well is also flexible enough to fit its noise, and the remedy is always some form of restriction on that flexibility. Chapter 2 of the SciML material makes the same point in physical terms, where the restriction comes from governing equations rather than from a penalty term.
(1) Feature Engineering
- The performance of regression models depends critically on the quality of the input features. Preprocessing techniques such as feature scaling, categorical encoding, and the construction of interaction terms can substantially improve model performance.
(2) Multivariate Regression
- An extension of linear regression in which multiple output variables (dependent variables) are predicted simultaneously from one or more input variables.
(3) Nonlinear Regression
Models the relationship between independent variables and dependent variables using nonlinear functions.
In contrast to linear regression, nonlinear models can represent complex, non-linear relationships present in data.
(4) Overfitting
Overfitting occurs when a regression model fits noise and incidental patterns in the training data rather than the underlying data-generating process.
Consequently, the model exhibits good performance on the training data but poor generalization to unseen (test) data.
(5) Regularization
Regularization incorporates a penalty term into the loss function to control model complexity and mitigate overfitting.
By discouraging excessively large parameter values, regularization improves generalization performance.
from IPython.display import YouTubeVideo
YouTubeVideo('IRrGVQV8vZ8?si=_hkRnh2jEhJVDXsY&start=931', width = '560', height = '315', extras = ['referrerpolicy="strict-origin-when-cross-origin"'])
Classification is another core task in machine learning. Its objective is to assign input data to one of several predefined classes based on observed features. Unlike regression, where the target variable is continuous, classification problems involve discrete output variables. Classification plays a central role in many real-world applications, including spam detection, image recognition, and medical diagnosis.
We begin with binary classification problems in order to develop and understand the core concepts underlying classification methods.
Why Not Simply Use Regression?
A natural first reaction is to encode the two classes as $0$ and $1$ and fit a straight line by least squares, as in Chapter 1. This can be made to work in easy cases, but it fails for reasons that motivate everything in this chapter.
The output of a line is unbounded. Given an input far from the origin, the fitted line will confidently predict $3.7$ or $-2.1$, values that no label can take and that cannot be interpreted.
The squared error penalizes points for being correct. A point of class $1$ lying far on the correct side of the boundary produces a large residual, so least squares will tilt the boundary in order to reduce an "error" that is not an error at all.
What we actually need is a decision, not a numerical value. The quantity of interest is which side of a boundary a point falls on, and possibly how confident we are about it. Chapter 1 gave us the wrong objective for that question.
The methods below correct this in two stages. The perceptron discards the magnitude of the prediction and keeps only its sign. Logistic regression restores the magnitude, but passes it through a squashing function so that it becomes a probability rather than an unbounded number.
In a typical machine learning course, commonly studied classification algorithms include the perceptron, support vector machines (SVMs), and logistic regression. In this chapter, however, we omit a detailed treatment of SVMs. Instead, we focus on the perceptron and logistic regression. These two models are particularly important because they form the conceptual and mathematical foundation of modern artificial neural networks.
The reason for this choice will become concrete by the end of the chapter. A single perceptron is exactly one artificial neuron, and logistic regression is exactly one neuron with a sigmoid activation trained by cross-entropy. A deep network is these two objects stacked and composed. The SVM, by contrast, is a beautiful and important method, but it leads in a different direction and is not a building block of the architectures we will need later.
The Perceptron is one of the simplest types of artificial neural networks and is used for binary classification. It was invented by Frank Rosenblatt in 1958 and serves as the foundation for more advanced neural networks.
To understand how the perceptron works, let us consider a bank loan approval scenario where the bank decides whether to approve or reject a loan application based on specific criteria.
The example is worth taking seriously rather than treating as a formality, because it makes the structure of the model transparent. Each applicant is described by a list of numbers: annual income, credit score, years of employment, existing debt, and so on. The bank forms a single score by weighting each of these and adding them up, then compares that score against a threshold. Applicants above the threshold are approved.
Everything the perceptron does is contained in that sentence. The weights $\omega_i$ encode how much each attribute matters and in which direction, so a positive weight on income and a negative weight on existing debt would both be expected. The threshold sets how demanding the bank is. What makes this a machine learning problem rather than a policy decision is that we do not set these numbers ourselves. We infer them from records of past applications and their outcomes.
- For input $x = \begin{bmatrix}x_1\\ \vdots\\ x_d \end{bmatrix}\;$ 'attributes (or features) of a customer'
- weights $\omega = \begin{bmatrix}\omega_1\\ \vdots\\ \omega_d \end{bmatrix}$
Suppose a decision can be made based on a simple linear combination of features:
$$\begin{align*} \text{Approve credit if} \; & \sum\limits_{i=1}^{d}\omega_ix_i > \text{threshold}, \\ \text{Deny credit if} \; & \sum\limits_{i=1}^{d}\omega_ix_i < \text{threshold}. \end{align*}$$
$$g(x) = \text{sign} \left(\left( \sum\limits_{i=1}^{d}\omega_ix_i \right)- \text{threshold} \right) = \text{sign}\left(\left( \sum\limits_{i=1}^{d}\omega_ix_i \right)+ \omega_0\right)$$
where a sign function is:
The sign function is what converts a continuous score into a discrete decision. Note the consequence of applying it: all information about the size of the score is destroyed, and only its direction survives. An applicant who exceeds the threshold by a wide margin and one who barely exceeds it receive an identical output of $+1$. This is the defining limitation of the perceptron, and Section 3.1 begins by removing it.
$$ \text{sign}(x) = \begin{cases} 1, &\text{if }\; x > 0\\ 0, &\text{if }\; x = 0\\ -1, &\text{if }\; x < 0 \end{cases} $$
Introduce an artificial coordinate $x_0 = 1$:
- To simplify the perceptron formulation, we can include the bias term $\omega_0$ into the weight vector by adding an artificial coordinate $x_0 = 1$ to the input vector $x$. This eliminates the explicit bias term and simplifies the computation.
$$g(x) = \text{sign}\left( \sum\limits_{\color{red}{i=0}}^{d}\omega_ix_i \right)$$
In a vector form
$$g(x) = \text{sign}\left( \omega^T x \right)$$
Classification boundary: hyperplane $g(x) = 0 \implies \omega^T x = 0$
The word hyperplane simply means the straight generalization of a line to arbitrary dimension: a line in two dimensions, an ordinary plane in three, and a flat $(d-1)$-dimensional surface in $d$ dimensions. The listed properties follow directly from the equation $\omega^T x = 0$, and each is worth verifying rather than accepting.
Why $\omega$ is orthogonal to the hyperplane. If $x_a$ and $x_b$ both lie on the boundary, then $\omega^T x_a = 0$ and $\omega^T x_b = 0$, so $\omega^T (x_a - x_b) = 0$. The vector $x_a - x_b$ points along the boundary, and $\omega$ is perpendicular to it. Since $x_a$ and $x_b$ were arbitrary, $\omega$ is perpendicular to every direction lying within the boundary.
Why $\omega$ determines the two sides. The sign of $\omega^T x$ tells us whether $x$ lies on the side that $\omega$ points toward (positive) or the opposite side (negative). The vector $\omega$ therefore does two jobs at once: its direction orients the boundary, and it also labels which half-space is the positive class.
Why the augmented form is convenient. Without the artificial coordinate, the boundary is $\omega^T x + \omega_0 = 0$, an affine surface that need not pass through the origin. After appending $x_0 = 1$ and absorbing $\omega_0$ into the weight vector, the boundary becomes $\omega^T x = 0$, a hyperplane through the origin in the enlarged space. Nothing has been lost; the offset has simply been rewritten as one more weight so that a single formula covers both cases. This trick recurs constantly in machine learning implementations.
Separates a D-dimensional space into two half-spaces
Defined by an outward pointing normal vector $\omega$
Assume the hyperplane passes through origin, $\omega^T x = 0$ with $x_0 = 1$
$\omega$ is orthogonal to any vector lying on the hyperplane
$$ \begin{align*} \omega = \begin{bmatrix}\omega_1 \\ \omega_2 \end{bmatrix}, \quad x = \begin{bmatrix} x_1 \\ x_2\end{bmatrix} &\implies g(x) = \omega_0 + \omega_1 x_1 + \omega_2 x_2 = \omega_0 + \omega^T x \qquad \text{or}\\\\\\ \omega = \begin{bmatrix}\omega_0 \\ \omega_1 \\ \omega_2 \end{bmatrix}, \quad x = \begin{bmatrix} 1 \\ x_1 \\ x_2\end{bmatrix} &\implies g(x) = \omega_0 + \omega_1 x_1 + \omega_2 x_2 = \omega^T x \end{align*} $$
Given the vector $ \omega $, the decision boundary is defined by the hyperplane:
$$ \omega^T x = 0 $$
This hyperplane is orthogonal to the vector $ \omega $, meaning $ \omega $ determines the direction perpendicular to the boundary.
There is one more fact about $\omega^T x$ that will be needed in Chapter 3, so it is worth recording here. For a point $x$, the quantity
$$ \frac{\omega^T x}{\lVert \omega \rVert} $$
is the signed perpendicular distance from $x$ to the hyperplane. Its magnitude is the distance, and its sign indicates the side.
This means $\omega^T x$ is not merely a device for producing a sign. Up to the scale factor $\lVert \omega \rVert$, it measures how far a point lies from the decision boundary, and it is a distance to a surface rather than a distance between two data points. The perceptron throws this magnitude away by applying the sign function. Logistic regression keeps it, which is exactly the difference between the two methods.
Blue Points: Since all blue points are located on the same side of the hyperplane relative to $ \omega $, they satisfy:
$$ \omega^T x > 0 $$
This corresponds to the region where $ g(x) > 0 $, classifying these points as belonging to Class 1.
Green Points: Conversely, the green points are positioned on the opposite side of the hyperplane relative to $ \omega $, satisfying:
$$ \omega^T x < 0 $$
This corresponds to the region where $ g(x) < 0 $, classifying these points as belonging to Class 0.
In this linear classification framework, $ \omega $ not only defines the orientation of the boundary but also plays a crucial role in determining which points belong to each class based on their position relative to the hyperplane.
Learning a Hyperplane for Classification
(1) Goal:
The objective is to learn the hyperplane $g_{\omega}(x)=0$ using the given training data. This hyperplane serves as the decision boundary that separates different classes.
(2) How to find $\omega$:
The parameter vector $\omega$ defines the orientation and position of the hyperplane. The classification rule is determined as follows:
- For all data points belonging to Class 1:
$$g(x) > 0$$
- For all data points belonging to Class 0:
$$g(x) < 0$$
This formulation highlights the fundamental idea behind linear classifiers, where the decision boundary is defined by a linear function of the input features. By appropriately determining $\omega$, we can effectively separate the data points into their respective classes.
Notice what kind of problem this is. We are given the data and asked to find a vector $\omega$ such that a list of inequalities is satisfied simultaneously, one inequality per training point. If such an $\omega$ exists, the dataset is called linearly separable. If the two classes interlock in a way that no straight boundary can resolve, no $\omega$ satisfies all the inequalities and the problem as posed has no solution. Keep this possibility in mind, because it determines whether the algorithm of the next section terminates.
Note that in traditional methods, the vector $\omega$ is often determined by domain experts based on their prior knowledge and insights about the system. However, in machine learning, the objective is to estimate $\omega$ directly from the data without any inherent bias or preconception. This data-driven approach allows the model to adaptively learn the optimal decision boundary, potentially uncovering complex patterns that may not be apparent through manual design.
Key Insights
'Learning from Data' as a Paradigm Shift:
The concept of "learning from data" represents a significant philosophical shift for many engineers. Traditional engineering methods often rely on explicitly defined models and expert-derived parameters. In contrast, machine learning emphasizes empirical learning, where patterns and relationships are discovered directly from observed data.
Estimating $\omega$:
There are various techniques available to estimate $\omega$, each designed to suit different types of data and modeling requirements. We will explore several key algorithms for effectively estimating $\omega$. These methods form the foundation of machine learning approaches to classification.
2.2.1. Perceptron Algorithm¶
We will first walk through the perceptron algorithm and then explore the underlying principles that explain how and why it works.
The perceptron implements
$$h(x) = \text{sign}\left( \omega^Tx \right)$$
Given the training set
$$(x_1, y_1), (x_2, y_2), \cdots, (x_N, y_N) \quad \text{where } y_i \in \{-1,1\}$$
- pick a misclassified point
$$ \text{sign}\left(\omega^Tx_n \right) \neq y_n$$
- and update the weight vector
$$\omega \leftarrow \omega + y_nx_n$$
The two steps are repeated until no misclassified point remains, or until a prescribed maximum number of iterations is reached.
Read the update rule carefully, because it is remarkably economical. It contains no learning rate, no loss function, and no derivative. When a point of class $+1$ is missed, the misclassified input is added to the weight vector; when a point of class $-1$ is missed, it is subtracted. In words: rotate the weight vector slightly toward the point that was gotten wrong. Section 2.2.1 explains why this simple gesture is guaranteed to help.
Convergence
The perceptron has a genuine theoretical guarantee, known as the perceptron convergence theorem: if the training data is linearly separable, the algorithm finds a separating hyperplane in a finite number of updates. The bound on that number depends on the margin of the data and not on the number of training points.
The converse is equally important in practice. If the data is not linearly separable, the algorithm never terminates. It will cycle indefinitely, since at every stage some point remains misclassified. This is why implementations always impose a maximum iteration count, and why a perceptron that fails to converge is telling you something about your data rather than about your code.
Why Perceptron Updates Work ?
- Let us look at a misclassified positive example ($y_n = +1$)
- perceptron (wrongly) thinks $\omega_{old}^T x_n < 0$
- Updates would be
$$ \begin{align*}\omega_{new} &= \omega_{old} + y_n x_n = \omega_{old} + x_n \\ \\ \omega_{new}^T x_n &= (\omega_{old} + x_n)^T x_n = \omega_{old}^T x_n + x_n^T x_n \geq \omega_{old}^T x_n \end{align*}$$
- Thus $\omega_{new}^T x_n$ is less negative than $\omega_{old}^T x_n$
The increase equals $x_n^T x_n = \lVert x_n \rVert^2$, which is strictly positive for any nonzero input. The update therefore always moves the score for that point in the correct direction, by an amount determined by the magnitude of the input.
- The negative case works identically. For a misclassified example with $y_n = -1$, the perceptron wrongly believes $\omega_{old}^T x_n > 0$, and the update is $\omega_{new} = \omega_{old} - x_n$. Then
$$\omega_{new}^T x_n = \omega_{old}^T x_n - \lVert x_n \rVert^2 \leq \omega_{old}^T x_n ,$$
so the score becomes less positive, again the desired direction. Both cases are summarized by the single statement $y_n \omega_{new}^T x_n \geq y_n \omega_{old}^T x_n$.
Two cautions are in order. The update improves the score for the point currently under consideration, but it may break points that were previously classified correctly. Progress is not monotone in the number of correct classifications, which is why the convergence theorem requires a more careful argument than the one line above. In addition, the size of the update scales with $\lVert x_n \rVert^2$, so features measured on very different scales will produce updates of very different magnitudes. This is one of the standard reasons for normalizing input features before training.
Note:
(1) Perceptron Update Rule (discrete version):
$$\omega \leftarrow \omega \pm 1 \cdot x_n$$
(2) Gradient Descent Update Rule (continuous version):
$$\omega \leftarrow \omega - \alpha \nabla_{\omega} f$$
Understanding this similarity provides insight into the foundational mechanisms of various learning algorithms, illustrating how they iteratively adjust model parameters to enhance performance.
The parallel is closer than it may appear, and the next section makes it exact. Both rules have the same shape: start from the current parameters, and add a correction term computed from a single training example. The difference is only in where the correction comes from. The perceptron obtains it from a hand-designed rule justified by the geometric argument above, while gradient descent obtains it by differentiating an explicitly stated loss function.
The second approach is the one that generalizes. Once a loss is written down, the same optimization machinery applies to any differentiable model, which is why every method after this section is formulated as the minimization of a loss rather than as a bespoke update rule.
Summary
The perceptron is a simple yet powerful model for binary classification tasks such as bank loan approval, classifying applicants on the basis of features such as credit score, income, and debt. While it performs well on linearly separable data, it has clear limitations on more complex datasets. Understanding how the perceptron learns and updates its weights is fundamental to understanding modern neural networks.
Before moving on, it is worth consolidating what has been established.
- The model is $h(x) = \text{sign}(\omega^T x)$, a linear score followed by a hard threshold.
- The decision boundary is the hyperplane $\omega^T x = 0$, oriented by the vector $\omega$.
- Learning consists of repeatedly selecting a misclassified point and applying $\omega \leftarrow \omega + y_n x_n$.
- If the data is linearly separable the procedure terminates; otherwise it does not.
- The hyperplane it returns is some separating hyperplane, not the best one.
The last two points are the open issues. Section 2.2.3 addresses the second, and Chapter 3 addresses the first by replacing the hard threshold with a smooth function.
This diagram symbolically illustrates the Perceptron, a fundamental model in machine learning. The Perceptron operates as follows:
(1) Input Stage
The model takes multiple inputs, denoted as $x_0, x_1, \cdots, x_d$, each associated with a corresponding weight $\omega_0, \omega_1, \cdots, \omega_d$.
These weights are key parameters that the model adjusts during the learning process.
(2) Weighted Sum Calculation
- Each input is multiplied by its respective weight, and the results are summed:
$$ z = \omega_0 x_0 + \omega_1 x_1 + \cdots + \omega_d x_d $$
(3) Activation Function
- The summation result $ z $ is passed through the sign function, which determines the predicted output:
$$ \hat{y} = \text{sign}(z) $$
- The sign function outputs either $ +1 $ or $ -1 $, representing the two possible classes in binary classification.
(4) Weight Update (Learning Process)
- The predicted output $ \hat{y} $ is compared with the true label $ y $.
- If the prediction is incorrect, the weights are updated as follows:
$$ \omega \leftarrow \omega + y_n x_n $$
- This iterative adjustment is the foundation of the Perceptron Learning Algorithm, ensuring that the decision boundary improves over time based on the training data. This is a clear example of learning from data.
Important Insight for Deep Learning
This diagram is also significant because the Perceptron can be viewed as a simple neuron in an Artificial Neural Network (ANN). When we study deep learning later, this structure will serve as the fundamental building block of more complex models. Understanding this visualization now will provide valuable intuition when exploring neural networks and their hierarchical architecture.
The mapping between the two vocabularies is exact, and it is worth memorizing now so that later chapters read as familiar material rather than as new material:
| Perceptron | Neural network |
|---|---|
| input $x_0, \cdots, x_d$ | input layer |
| weights $\omega_0, \cdots, \omega_d$ | connection weights |
| weighted sum $z = \omega^T x$ | pre-activation |
| $\text{sign}(z)$ | activation function |
| update $\omega \leftarrow \omega + y_n x_n$ | one step of training |
A deep network changes exactly two things about this picture. The sign function is replaced by a smooth activation, such as the sigmoid of Chapter 3 or the ReLU mentioned in the next section, so that derivatives exist. And many such units are stacked in layers, so that the output of one becomes the input of the next. Nothing else about the diagram changes, which is why the perceptron is the right place to start.
Perceptron loss function
Rather than explicitly checking whether each point is misclassified, we can express the same rule through a loss function that automatically accumulates contributions only from the misclassified points.
The loss for an individual sample is defined as:
$$ \mathscr{L}(\omega) = \max \left\{ 0, -y_n \cdot \left(\omega^T x_n \right)\right\} $$
- Loss $\mathscr{L}(\omega) = 0$ on examples where perceptron is correct, i.e., $y_n \cdot \left(\omega^T x_n \right) > 0$
- Loss $\mathscr{L}(\omega) > 0$ on examples where perceptron misclassified, i.e., $y_n \cdot \left(\omega^T x_n \right) < 0$
Note:
$\text{sign}\left(\omega^T x_n \right) \neq y_n$ is equivalent to $ y_n \cdot \left(\omega^T x_n \right) < 0$
$\text{ReLU}(z) = \max(0, z)$: the Rectified Linear Unit (ReLU) will be revisited later in the discussion, as it plays a crucial role not only in perceptron loss formulation but also in modern deep learning architectures where it is widely used as an activation function.
This function returns zero when the point is correctly classified and returns a positive value proportional to the margin violation when misclassified.
Compact Expression for Total Perceptron Loss
The total loss across all samples in the dataset $D$ is given by:
$$ \mathscr{L}(\omega) = \sum_{n =1}^{m} \max \left\{ 0, -y_n \cdot \left(\omega^T x_n \right)\right\} $$
This summation aggregates the loss over all training samples and only penalizes incorrectly classified points.
Reading the Loss Function
The quantity $y_n \cdot (\omega^T x_n)$ appearing inside is worth a name of its own. It is called the margin of the point, and it combines correctness and confidence into one number: it is positive when the prediction is right, negative when the prediction is wrong, and large in magnitude when the point lies far from the boundary. The perceptron loss is then simply $\max\{0, -\text{margin}\}$, which is zero for correct points and grows linearly with how badly a wrong point is wrong.
Recovering the Update Rule
The reason for introducing the loss is that the perceptron update now follows from calculus rather than from geometry. For a misclassified point the loss is $-y_n \omega^T x_n$, whose gradient with respect to $\omega$ is
$$\nabla_{\omega} \mathscr{L} = -y_n x_n .$$
For a correctly classified point the loss is identically zero and the gradient vanishes. Substituting into the gradient descent rule $\omega \leftarrow \omega - \alpha \nabla_{\omega}\mathscr{L}$ with a step size of $\alpha = 1$ gives
$$\omega \leftarrow \omega + y_n x_n ,$$
which is precisely the perceptron update introduced earlier. The hand-designed rule and the gradient of the loss are the same object seen from two directions. This is the first instance of a pattern that will hold for every remaining model in the course: specify a loss, differentiate it, and the learning rule is whatever comes out.
One technical remark. The function $\max\{0, z\}$ is not differentiable at $z = 0$, so strictly speaking the derivative above is a subgradient rather than a gradient. In practice this causes no difficulty, since the non-differentiable point is a single value that optimization almost never lands on exactly. The same remark applies to the ReLU activation in deep learning, which is the same function.
2.2.2. Perceptron in Python¶
By adding an artificial coordinate $x_0 = 1$, we simplify the perceptron implementation and avoid the need to handle the bias separately. This technique is commonly used in machine learning to streamline calculations and make the formulation more uniform.
Pay attention to the shape conventions below, since they are used consistently in the code that follows and mismatched shapes are the most common source of errors in practice. Each row of the matrix $x$ is one training example, and each column is one feature. With $m$ examples and two original features, $x$ is $m \times 3$ after the leading column of ones is prepended, $\omega$ is $3 \times 1$, and the product $x\omega$ is $m \times 1$, giving one score per example. The superscript $(i)$ indexes the example, while the subscript indexes the feature.
$$g(x) = \omega_0 + \omega^Tx = \omega_0 + \omega_1x_1 + \omega_2x_2 = 0$$
$$ \begin{align*} \omega &= \begin{bmatrix} \omega_0 \\ \omega_1 \\ \omega_2\end{bmatrix}\\ \\ x &= \begin{bmatrix} \left(x^{(1)}\right)^T \\ \left(x^{(2)}\right)^T \\ \left(x^{(3)}\right)^T\\ \vdots \\ \left(x^{(m)}\right)^T \end{bmatrix} = \begin{bmatrix} 1 & x_1^{(1)} & x_2^{(1)} \\ 1 & x_1^{(2)} & x_2^{(2)} \\ 1 & x_1^{(3)} & x_2^{(3)}\\\vdots & \vdots & \vdots \\ 1 & x_1^{(m)} & x_2^{(m)}\end{bmatrix} \\ \\ y &= \begin{bmatrix}y^{(1)} \\ y^{(2)} \\ y^{(3)}\\ \vdots \\ y^{(m)} \end{bmatrix} \end{align*}$$
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
#training data gerneration
m = 100
x1 = 8*np.random.rand(m, 1)
x2 = 7*np.random.rand(m, 1) - 4
g = 0.8*x1 + x2 - 3
C1 = np.where(g >= 1)
C0 = np.where(g < -1)
print(C1)
C1 = np.where(g >= 1)[0]
C0 = np.where(g < -1)[0]
print(C1.shape)
print(C0.shape)
plt.figure(figsize = (6, 4))
plt.plot(x1[C1], x2[C1], 'ro', alpha = 0.4, label = 'C1')
plt.plot(x1[C0], x2[C0], 'bo', alpha = 0.4, label = 'C0')
plt.title('Linearly Separable Classes')
plt.legend(loc = 1)
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.show()
$$ \begin{align*} x &= \begin{bmatrix} \left(x^{(1)}\right)^T \\ \left(x^{(2)}\right)^T \\ \left(x^{(3)}\right)^T\\ \vdots \\ \left(x^{(m)}\right)^T \end{bmatrix} = \begin{bmatrix} 1 & x_1^{(1)} & x_2^{(1)} \\ 1 & x_1^{(2)} & x_2^{(2)} \\ 1 & x_1^{(3)} & x_2^{(3)}\\\vdots & \vdots & \vdots \\ 1 & x_1^{(m)} & x_2^{(m)}\end{bmatrix} \\\\ y &= \begin{bmatrix}y^{(1)} \\ y^{(2)} \\ y^{(3)}\\ \vdots \\ y^{(m)} \end{bmatrix} \end{align*}$$
X1 = np.hstack([np.ones([C1.shape[0],1]), x1[C1], x2[C1]])
X0 = np.hstack([np.ones([C0.shape[0],1]), x1[C0], x2[C0]])
X = np.vstack([X1, X0])
y = np.vstack([np.ones([C1.shape[0],1]), -np.ones([C0.shape[0],1])])
X = np.asmatrix(X)
y = np.asmatrix(y)
$$\omega = \begin{bmatrix} \omega_0 \\ \omega_1 \\ \omega_2\end{bmatrix}$$
$$\omega \leftarrow \omega + yx$$ where $(x, y)$ is a misclassified training point
The implementation follows the algorithm literally. Initialize $\omega$, then loop: compute the predictions, find a point whose prediction disagrees with its label, apply the update, and repeat. Because every update depends on a single training example, the perceptron is an example of an online learning algorithm, one that can process data one point at a time rather than requiring the whole dataset at once.
w = np.ones([3,1])
w = np.asmatrix(w)
n_iter = y.shape[0]
for k in range(n_iter):
for i in range(n_iter):
if y[i,0] != np.sign(X[i,:]*w)[0,0]:
w += y[i,0]*X[i,:].T
print(w)
$$ \begin{align*} g(x) &= \omega_0 + \omega^Tx = \omega_0 + \omega_1x_1 + \omega_2x_2 = 0 \\\\ \implies x_2 &= -\frac{\omega_1}{\omega_2} x_1 - \frac{\omega_0}{\omega_2} \end{align*} $$
The rearrangement is purely for plotting. The learned model is the vector $\omega$, and the boundary is defined implicitly by $\omega^T x = 0$. To draw it on a two-dimensional scatter plot we solve that implicit equation for $x_2$ as an explicit function of $x_1$, which is the familiar slope-intercept form of a line. The slope is $-\omega_1/\omega_2$ and the intercept is $-\omega_0/\omega_2$. Note that this form breaks down when $\omega_2 = 0$, which corresponds to a vertical boundary, a reminder that the implicit form $\omega^T x = 0$ is the more general description.
x1p = np.linspace(0,8,100).reshape(-1,1)
x2p = - w[1,0]/w[2,0]*x1p - w[0,0]/w[2,0]
plt.figure(figsize = (6, 4))
plt.plot(x1[C1], x2[C1], 'ro', alpha = 0.4, label = 'C1')
plt.plot(x1[C0], x2[C0], 'bo', alpha = 0.4, label = 'C0')
plt.plot(x1p, x2p, c = 'k', linewidth = 3, label = 'perceptron')
plt.xlim([0, 8])
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.legend(loc = 1)
plt.show()
Perceptron using Scikit-Learn
In scikit-learn, the Perceptron includes the bias term by default, so you do not need to add an artificial coordinate manually. This makes it more convenient for implementing and training perceptron models on datasets.
Compare the matrix below with the one used in the manual implementation: the leading column of ones is absent. This is a small but frequent source of confusion when moving between hand-written code and library code, so it is worth stating the rule explicitly. When you form the design matrix yourself, you must supply the constant column. When a library manages the intercept internally, supplying it yourself creates a redundant feature.
Note also that the fitted parameters are then reported in two separate places, coef_ for the weights and intercept_ for the bias, rather than in a single augmented vector.
$$ \begin{align*} x &= \begin{bmatrix} \left(x^{(1)}\right)^T \\ \left(x^{(2)}\right)^T \\ \left(x^{(3)}\right)^T\\ \vdots \\ \left(x^{(m)}\right)^T \end{bmatrix} = \begin{bmatrix} x_1^{(1)} & x_2^{(1)} \\ x_1^{(2)} & x_2^{(2)} \\ x_1^{(3)} & x_2^{(3)}\\ \vdots & \vdots \\ x_1^{(m)} & x_2^{(m)}\end{bmatrix} \\\\ y &= \begin{bmatrix}y^{(1)} \\ y^{(2)} \\ y^{(3)}\\ \vdots \\ y^{(m)} \end{bmatrix} \end{align*}$$
X1 = np.hstack([x1[C1], x2[C1]])
X0 = np.hstack([x1[C0], x2[C0]])
X = np.vstack([X1, X0])
y = np.vstack([np.ones([C1.shape[0],1]), -np.ones([C0.shape[0],1])])
from sklearn import linear_model
clf = linear_model.Perceptron(tol = 1e-3)
clf.fit(X, np.ravel(y))
clf.predict([[3, -2]])
clf.predict([[6, 2]])
clf.coef_
clf.intercept_
$$ \begin{align*} g(x) &= \omega_0 + \omega^Tx = \omega_0 + \omega_1x_1 + \omega_2x_2 = 0 \\\\ \implies x_2 &= -\frac{\omega_1}{\omega_2} x_1 - \frac{\omega_0}{\omega_2} \end{align*} $$
w0 = clf.intercept_[0]
w1 = clf.coef_[0,0]
w2 = clf.coef_[0,1]
x1p = np.linspace(0,8,100).reshape(-1,1)
x2p = - w1/w2*x1p - w0/w2
plt.figure(figsize = (6, 4))
plt.plot(x1[C1], x2[C1], 'ro', alpha = 0.4, label = 'C1')
plt.plot(x1[C0], x2[C0], 'bo', alpha = 0.4, label = 'C0')
plt.plot(x1p, x2p, c = 'k', linewidth = 4, label = 'perceptron')
plt.xlim([0, 8])
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.legend(loc = 1)
plt.show()
2.2.3. The Best Hyperplane Separator?¶
Limitations
The perceptron algorithm is guaranteed to find a separating hyperplane if the data is linearly separable. However, there are typically many possible hyperplanes that can separate the data. The perceptron algorithm does not guarantee the optimal hyperplane; it merely returns one of the many admissible solutions.
Improvements
While the perceptron identifies one valid separating hyperplane, identifying the best hyperplane requires an optimization framework. This concept is crucial in understanding more advanced models like Support Vector Machines (SVM) and Logistic Regression, which are designed to find optimal decision boundaries with improved performance and generalization.
The issue is easy to see and is worth thinking through before reading on. Suppose the two classes are well separated, with a visible empty corridor between them. Any hyperplane drawn inside that corridor classifies every training point correctly, and the perceptron will stop at whichever one it happens to reach first, which depends on the initial weights and on the order in which misclassified points were selected. A boundary hugging the edge of one class and a boundary running down the middle of the corridor are indistinguishable to the algorithm, since both achieve zero training error.
They are not equally good, however. Training error is not the quantity we care about. We care about performance on data we have not yet seen, and a boundary that passes very close to the training points of one class will misclassify a new point from that class that happens to fall slightly on the other side. The boundary running down the middle of the corridor leaves room for such variation on both sides.
Making this intuition precise requires stating an objective that distinguishes among hyperplanes that are all correct. The SVM does so by maximizing the width of the corridor, called the margin. Logistic regression, which we take up next, does so by a different route: it assigns a probability to every point and maximizes the likelihood of the observed labels, with the effect that a boundary comfortably separating the data is preferred over one that barely does.
2.2.4. Historical Notes¶
In the late 1960s, Marvin Minsky and Seymour Papert published the influential book "Perceptrons" (1969), where they mathematically demonstrated the limitations of the perceptron. One key result from their work was that single-layer perceptrons cannot solve the XOR problem (exclusive OR), highlighting a significant limitation of early neural networks.
The XOR example is small enough to verify by hand, and doing so is instructive. Consider four points with the labels of the exclusive-or function:
| $x_1$ | $x_2$ | $y$ |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 1 | 1 |
| 1 | 0 | 1 |
| 1 | 1 | 0 |
The two points labelled $1$ lie on one diagonal of the square and the two labelled $0$ lie on the other. No single straight line can place both members of one diagonal on one side and both members of the other diagonal on the opposite side. The dataset is not linearly separable, and by the discussion in Section 2.2.1 the perceptron algorithm will therefore never terminate on it.
This is a limitation of the model, not of the training procedure, and no cleverer update rule can repair it. What repairs it is a change of representation. If a third feature $x_1 x_2$ is added, the four points move into three dimensions and become linearly separable there. A hidden layer in a neural network does precisely this, except that it learns the additional features from data instead of requiring a person to invent them.
Stagnation in Neural Network Research:
- The book led to a significant decline in interest and funding for neural networks, as it was seen as proof that neural networks had severe limitations.
- This period is often referred to as the "AI Winter".
Rediscovery with Multi-Layer Networks:
- In the 1980s, researchers like Geoffrey Hinton showed that multi-layer perceptrons (MLPs) with backpropagation could learn non-linear functions like XOR.
- Adding hidden layers and non-linear activation functions allows neural networks to learn complex, non-linear patterns.
Logistic regression is fundamentally a classification algorithm, even though the term "regression" appears in its name.
The previous chapter ended with two unresolved difficulties. The perceptron discards the magnitude of $\omega^T x$ and keeps only its sign, so it cannot distinguish a confident prediction from a marginal one. And because every separating hyperplane achieves zero training error, the algorithm has no basis for preferring one over another.
Logistic regression resolves both with a single change. Instead of thresholding the score, it passes the score through a smooth squashing function whose output lies between $0$ and $1$ and can be read as a probability. Distance information is retained, every training point contributes to the objective, and the resulting optimization problem is convex and therefore has a unique solution.
Section 3.1 motivates the construction geometrically, arriving at the sigmoid function and at the maximum likelihood objective. Section 3.2 applies the method using scikit-learn. Section 3.3 revisits the resulting loss function from the perspective of entropy, connecting it to a concept familiar from physics.
from IPython.display import YouTubeVideo
YouTubeVideo('IRrGVQV8vZ8?si=8H3cmkAUyNIu2NDb&start=2750', width = '560', height = '315', extras = ['referrerpolicy="strict-origin-when-cross-origin"'])
Perceptron: Utilizes only the sign of the data to determine class labels, relying solely on binary outcomes without considering the distance of points from the decision boundary.
Logistic Regression: To better utilize the distance information of all data points, we turn to logistic regression, which considers the positions of all points in the dataset. This broader perspective enhances the model's ability to capture patterns and improve generalization.
Recall from Section 2.2 that $\omega^T x$, divided by $\lVert \omega \rVert$, is the signed perpendicular distance from a point to the decision boundary. The perceptron computes this quantity and then throws away everything except its sign. The question driving this section is what could be done with the discarded information.
Let us start with the case of only two data points. Which linear classification boundary would be considered better?
On the left: The classification boundary is positioned near the center of the data points.
On the right: The classification boundary is biased toward one of the data points.
The boundary on the right is positioned too close to one of the data points, leaving minimal space (margin) on that side. This smaller margin increases the risk of misclassifying new data points that may fall within this narrow buffer zone. The boundary that maintains a balanced margin between the data points is preferred for its improved stability and performance.
Basic idea:
- Identify the decision boundary $g(x)=\omega^T x =0$ that maximizes the product of the distances $\lvert h(x) \rvert$ from the hyperplane to the nearest points of each class:
$$\prod_i \lvert h_i \rvert$$
Why?
- By the Arithmetic Mean-Geometric Mean (AM-GM) Inequality, for non-negative real numbers $\lvert h_1 \rvert$ and $\lvert h_2 \rvert$:
$$ \frac{\lvert h_1 \rvert + \lvert h_2 \rvert}{2} \geq \sqrt{\lvert h_1 \rvert \lvert h_2 \rvert} $$
- Equality holds if and only if $\lvert h_1 \rvert = \lvert h_2 \rvert$. Therefore, maximizing $\prod_i \lvert h_i \rvert$ encourages the distances $\lvert h_i \rvert$ to be equal, effectively positioning the hyperplane equidistantly between the two classes.
Roughly speaking, this optimization of $\max \prod_i \lvert h_i \rvert$ tends to position a hyperplane in the middle of two classes
It is worth being explicit about why the product is used rather than the sum. Maximizing $\sum_i \lvert h_i \rvert$ would be satisfied by pushing the boundary far away from one point while allowing another to sit almost on top of it, since a large term can compensate for a small one. A product cannot be rescued this way: if any single $\lvert h_i \rvert$ approaches zero, the entire product approaches zero regardless of how large the other factors are. The product therefore penalizes any point that comes close to the boundary, and the AM-GM inequality makes this precise by showing that, for a fixed sum, the product is largest when all the factors are equal.
The final line records that maximizing over $\omega^T x$ rather than over the properly normalized distance changes nothing essential, since the two differ only by the factor $1/\lVert \omega \rVert$ that is common to all points.
$$h = \frac{g(x)}{\lVert \omega \rVert} = \frac{\omega^T x}{\lVert \omega \rVert} \sim \omega^T x$$
Considering all data points may seem more effective than an approach such as the support vector machine, which relies only on the support vectors near the margin. However, some data points may be located far from the decision boundary. These extreme points (or outliers) can negatively impact the model’s performance by overly influencing the decision boundary due to the multiplication of all $\lvert h_i \rvert$.
To mitigate this negative impact, we apply a squeezing function that compresses the range of extreme values, mapping $(-\infty, +\infty)$ to $(0,1)$. This transformation ensures that points far from the decision boundary have minimal impact on the model's outcome.
By squeezing the values in this way, we create a more robust model that effectively emphasizes the most informative points near the decision boundary while reducing the adverse effects of outliers.
It is worth being clear about what the squashing function must accomplish, because these requirements determine the choice uniquely up to minor variations.
- It must be bounded, so that a point very far from the boundary cannot dominate the objective.
- It must be monotonically increasing, so that a larger distance still corresponds to greater confidence and the ordering of points is preserved.
- It must be smooth, so that the objective can be differentiated and optimized by gradient descent.
- Its output should lie in $(0, 1)$, so that it admits an interpretation as a probability.
The logistic function introduced next satisfies all four.
Here, let us examine the sigmoid function, also known as the logistic function, denoted as $\sigma(z)$, which often serves as a squeezing function.
$$ \sigma(z) = \frac{1}{1+e^{-z}} \implies \sigma \left(\omega^T x \right) = \frac{1}{1+e^{-\omega^T x}}$$
- The logistic function always returns a value between 0 and 1
- Crosses 0.5 at the origin, then flattens out
- Monotonic: the ordering of the distances is preserved, so the resulting optimization problem is essentially unchanged
- Continuous and differentiable: suitable for gradient descent optimization
- Reduces the influence of points lying far from the decision boundary
Two Properties Worth Verifying
The following identities are used repeatedly, and both can be checked in a line of algebra.
$$ 1 - \sigma(z) = \sigma(-z) $$
This symmetry says that the probability assigned to the negative class at $z$ equals the probability assigned to the positive class at $-z$, so the two classes are treated on an equal footing.
$$ \frac{d\sigma}{dz} = \sigma(z)\left(1 - \sigma(z)\right) $$
The derivative of the sigmoid is expressible in terms of the sigmoid itself, which makes gradient computations unusually clean and is the reason the formula in Section 3.1 turns out as simple as it does. Note also that this derivative is largest at $z = 0$ and decays toward zero as $\lvert z \rvert$ grows. Points near the boundary exert the most influence on learning, while points far from it exert almost none, which is exactly the behavior the previous section asked for.
A word of caution for later sections: this same decay is the origin of the vanishing gradient problem in deep networks, where sigmoid activations stacked in many layers multiply small derivatives together until learning stalls. It is one of the main reasons ReLU replaced the sigmoid in hidden layers of modern architectures.
# plot a sigmoid function
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
z = np.linspace(-4, 4, 100)
s = 1/(1 + np.exp(-z))
plt.figure(figsize = (6, 3))
plt.plot(z, s, linewidth = 3)
plt.xlim([-4, 4])
plt.axis('equal')
plt.grid(alpha = 0.3)
plt.show()
The output of the logistic function is bounded between 0 and 1, making it suitable for binary classification tasks.
- For logistic regression, the output label is $y_i \in \{0, 1\}$. Note that this convention differs from the perceptron, where the labels were $y_i \in \{-1, +1\}$
The logistic function compresses large positive and negative values:
- Large positive values are mapped close to 1.
- Large negative values are mapped close to 0.
One significant advantage of the logistic function is that its output can be interpreted as a probability.
- Often we do not care about predicting the label $y$
- Rather, we want to predict the label probabilities $P\left(y \mid x\,;\omega\right)$
$$P\left(y = 1 \mid x\,;\omega\right) = \frac{1}{1+e^{-\omega^T x}} \;\; \in \; [0,1]$$
- the probability that the label is $1$
$$P\left(y = 1 \mid x\,;\omega\right)$$
- the probability that the label is $0$
$$P\left(y = 0 \mid x\,;\omega\right) = 1 - P\left(y = 1 \mid x\,;\omega\right)$$
This reinterpretation is the conceptual turning point of the section, so it is worth stating plainly what has changed. We began with a geometric quantity, the distance from a point to a hyperplane. We compressed it with the sigmoid to keep outliers from dominating. The compressed value happens to lie in $(0,1)$ and to increase monotonically with the distance, so nothing prevents us from reading it as a probability.
Having done so, we gain access to an entire framework. Probabilities can be combined by multiplication across independent observations, they can be compared using likelihood, and parameters can be estimated by maximum likelihood estimation, a standard tool of statistics. The next several cells simply carry out that standard procedure. The geometric motivation has served its purpose and can now be set aside.
The predicted probability is also useful in its own right, independently of the classification decision. A medical screening model reporting a probability of $0.55$ and one reporting $0.99$ both predict the positive class, yet the two situations call for very different responses. The perceptron cannot distinguish them.
The Role of Distance in Classification
The perceptron, the support vector machine (SVM), and logistic regression are three fundamental approaches to linear classification. The primary distinction among them lies in how each defines and utilizes the distance from the decision boundary to the data points. Although the SVM is not treated in detail in this section, it is included here for comparison. The figures below illustrate these differences.
- Perceptron: Considers only the sign of the data points, completely ignoring the distance.
- SVM: Utilizes the distance, but only through the margin, that is, the distance from the decision boundary to the closest points (the support vectors).
- Logistic Regression: Considers the distance of all points, assigning probabilities accordingly, with points farther from the boundary contributing less significantly.
The three methods can be arranged on a single scale according to how much of the distance information they retain. The perceptron keeps none of it, using only the sign. The SVM keeps the distance of a few critical points and ignores the rest entirely. Logistic regression keeps all of the distances but weights them unequally, letting influence decay smoothly as a point moves away from the boundary.
Seen this way, the SVM and logistic regression are two different answers to the same question, namely how to prevent distant points from dominating. The SVM answers by truncation, discarding them outright. Logistic regression answers by attenuation, retaining them with diminishing weight.
Note:
Although the term "regression" appears in logistic regression, it is, in fact, a classification algorithm. This terminology may seem counterintuitive at first. However, upon closer examination, logistic regression can be understood as a regression-based approach in which the distances between data points and the linear decision boundary are transformed via a logistic (sigmoid) function. This transformation maps the outputs to a probability distribution within the range [0, 1], thereby enabling binary classification. Notably, the transformed distance is continuous, which aligns with the regression-like nature of the model.
Goal: we need to fit $\omega$ to our data
For a single data point $(x,y)$ with parameters $\omega$
$$ \begin{align*} P\left(y = 1 \mid x\,;\omega\right) &= h_{\omega}(x) = \sigma \left(\omega^T x \right)\\ P\left(y = 0 \mid x\,;\omega\right) &= 1 - h_{\omega}(x) = 1- \sigma \left(\omega^T x \right) \end{align*} $$
It can be compactly written as (since $y$ is either $0$ or $1$)
$$P\left(y \mid x\,;\omega\right) = \left(h_{\omega}(x) \right)^y \left(1 - h_{\omega}(x)\right)^{1-y}$$
The compact form is a standard trick and is easier to read once its mechanism is noticed. The exponents act as switches. When $y = 1$ the second factor becomes $(1 - h_\omega(x))^0 = 1$ and the expression reduces to $h_\omega(x)$; when $y = 0$ the first factor becomes $1$ and the expression reduces to $1 - h_\omega(x)$. A single formula thus covers both cases, which is what allows the product over all training points to be written down cleanly below.
For $m$ training data points, the likelihood function of the parameters:
$$ \begin{align*} \mathscr{L}(\omega) &= P\left(y^{(1)}, \cdots, y^{(m)} \mid x^{(1)}, \cdots, x^{(m)}\,;\omega\right)\\ &= \prod\limits_{i=1}^{m}P\left(y^{(i)} \mid x^{(i)}\,;\omega\right)\\ &= \prod\limits_{i=1}^{m}\left(h_{\omega}\left(x^{(i)}\right) \right)^{y^{(i)}} \left(1 - h_{\omega}\left(x^{(i)}\right)\right)^{1-y^{(i)}} \qquad \left(\sim \prod_i \lvert h_i \rvert \right) \end{align*} $$
The second line uses the assumption that the training examples are independent given the parameters, which is what permits a joint probability to be written as a product of individual probabilities.
The logic of maximum likelihood is worth stating in words, since the notation can obscure a simple idea. The expression $\mathscr{L}(\omega)$ is the probability that the model, with parameters $\omega$, would have produced exactly the labels we observed. Different values of $\omega$ make the observed data more or less probable. Maximum likelihood estimation selects the $\omega$ under which what actually happened was most likely to happen.
Note the parenthetical remark at the end of the line. The likelihood is a product over all data points, precisely the structure that the geometric argument of Section 3.1 arrived at when it proposed maximizing $\prod_i \lvert h_i \rvert$. The two lines of reasoning, one geometric and one probabilistic, converge on the same objective.
It would be easier to work on the log likelihood.
- Taking the logarithm converts the product into a sum
- Probabilities close to 0 can cause numerical instability when directly multiplied.
- The log transformation helps keep values in a manageable range.
$$\ell(\omega) = \log \mathscr{L}(\omega) = \sum_{i=1}^{m} y^{(i)} \log h_{\omega} \left(x^{(i)} \right) + \left(1-y^{(i)} \right) \log \left(1-h_{\omega} \left(x^{(i)} \right) \right)$$
- This is known as cross-entropy, a commonly used loss function in classification problems.
Two further reasons for taking the logarithm are worth adding to those listed above. The logarithm is monotonically increasing, so the value of $\omega$ that maximizes $\ell(\omega)$ is the same one that maximizes $\mathscr{L}(\omega)$; nothing is lost by the transformation. And the resulting sum is far easier to differentiate than a product of $m$ factors would be.
The structure of the summand is worth reading directly. For a point with $y^{(i)} = 1$, only the first term survives and the contribution is $\log h_\omega(x^{(i)})$. This is near zero when the model assigns high probability to the correct class, and large and negative when the model assigns a probability near zero. In other words, the loss punishes confident mistakes severely, and grows without bound as the assigned probability of the true class approaches zero. This is a deliberate and desirable property.
Then, the logistic regression problem can be solved as a (convex) optimization problem as
$$\hat{\omega} = \arg\max_{\omega} \ell(\omega)$$
The word convex in the preceding sentence carries real weight and should not be passed over. It means the objective has no local maxima other than the global one, so any optimization method that keeps improving will eventually reach the best possible solution, and the answer does not depend on where the search was started.
This is a stronger guarantee than we will have available for any later model in this course. Neural network objectives are not convex, and their training is a search over a landscape with many local minima. Logistic regression is therefore a useful reference point: it is the last model whose optimization we can fully trust.
Note finally that no closed-form solution exists here, unlike in Chapter 1. Setting the derivative of $\ell(\omega)$ to zero produces a system of transcendental equations with no algebraic solution, so the maximum must be found iteratively.
To apply the gradient descent method, we need the derivative of the log-likelihood.
Carrying out the differentiation and using the identity $\sigma'(z) = \sigma(z)(1 - \sigma(z))$ from Section 3.1, a substantial amount of algebra collapses and the result is remarkably simple:
$$ \frac{\partial \ell(\omega)}{\partial \omega} = \sum_{i=1}^{m} \left( y^{(i)} - h_{\omega}\left(x^{(i)}\right) \right) x^{(i)} $$
The gradient ascent update is therefore
$$ \omega \leftarrow \omega + \alpha \sum_{i=1}^{m} \left( y^{(i)} - h_{\omega}\left(x^{(i)}\right) \right) x^{(i)} $$
where $\alpha$ is the learning rate.
This deserves a moment of attention, because its form is familiar. The correction applied to $\omega$ is the input vector scaled by the prediction error $y^{(i)} - h_\omega(x^{(i)})$. A point the model already predicts correctly and confidently contributes almost nothing, since its error is near zero. A point the model gets badly wrong contributes an error approaching $\pm 1$ and moves the weights substantially.
Compare this with the perceptron update $\omega \leftarrow \omega + y_n x_n$. The structure is identical, an input vector scaled by a measure of error, and the difference is only that the perceptron's error is the discrete quantity $y_n$ evaluated at misclassified points, whereas logistic regression uses a continuous error evaluated at every point. This same form, "error times input", will appear once more as the backpropagation rule for a single layer of a neural network. It is not a coincidence, and recognizing it is a good check that the material has been understood.
We can implement logistic regression from scratch in Python, but instead, we will use scikit-learn for convenience and efficiency.
As with the perceptron, scikit-learn handles the intercept internally, so the design matrix below contains no column of ones. Two further points about the library implementation are worth knowing.
Scikit-learn applies $L_2$ regularization by default, adding a penalty on $\lVert \omega \rVert^2$ to the objective. The strength is controlled by the parameter
C, where a smallerCmeans stronger regularization. This is why the coefficients obtained from the library may differ slightly from those obtained by maximizing the unpenalized log-likelihood by hand.The method
predictreturns hard class labels obtained by thresholding at $0.5$, whilepredict_probareturns the underlying probabilities. The probabilities are usually the more informative output, and the threshold of $0.5$ is a convention that should be reconsidered whenever the two types of error carry different costs.
$$ \begin{align*} \omega &= \begin{bmatrix} \omega_1 \\ \omega_2\end{bmatrix}, \qquad \omega_0, \qquad x = \begin{bmatrix} x_1 \\ x_2\end{bmatrix}\\ \\ X &= \begin{bmatrix} \left(x^{(1)}\right)^T \\ \left(x^{(2)}\right)^T \\ \left(x^{(3)}\right)^T \\ \vdots\end{bmatrix} = \begin{bmatrix} x_1^{(1)} & x_2^{(1)} \\ x_1^{(2)} & x_2^{(2)} \\ x_1^{(3)} & x_2^{(3)} \\ \vdots & \vdots \\\end{bmatrix}\\ \\ y & = \begin{bmatrix} y^{(1)}\\ y^{(2)} \\y^{(3)} \\ \vdots \end{bmatrix} \end{align*} $$
# datat generation
m = 100
w0 = -6
w = np.array([[2], [1]])
X = np.hstack([4*np.random.rand(m,1), 4*np.random.rand(m,1)])
w = np.asmatrix(w)
X = np.asmatrix(X)
y = 1/(1 + np.exp(-w0-X*w)) > 0.5
C1 = np.where(y == True)[0]
C0 = np.where(y == False)[0]
y = np.empty([m,1])
y[C1] = 1
y[C0] = 0
plt.figure(figsize = (6, 4))
plt.plot(X[C1,0], X[C1,1], 'ro', alpha = 0.3, label = 'C1')
plt.plot(X[C0,0], X[C0,1], 'bo', alpha = 0.3, label = 'C0')
plt.xlabel(r'$x_1$', fontsize = 15)
plt.ylabel(r'$x_2$', fontsize = 15)
plt.legend(loc = 1, fontsize = 12)
plt.axis('equal')
plt.xlim([0,4])
plt.ylim([0,4])
plt.show()
from sklearn import linear_model
clf = linear_model.LogisticRegression(solver = 'lbfgs')
clf.fit(np.asarray(X), np.ravel(y))
clf.coef_
clf.intercept_
w0 = clf.intercept_[0]
w1 = clf.coef_[0,0]
w2 = clf.coef_[0,1]
xp = np.linspace(0,4,100).reshape(-1,1)
yp = - w1/w2*xp - w0/w2
plt.figure(figsize = (6, 4))
plt.plot(X[C1,0], X[C1,1], 'ro', alpha = 0.3, label = 'C1')
plt.plot(X[C0,0], X[C0,1], 'bo', alpha = 0.3, label = 'C0')
plt.plot(xp, yp, 'g', linewidth = 4, label = 'Logistic Regression')
plt.title('Logistic Regression')
plt.xlabel(r'$x_1$', fontsize = 15)
plt.ylabel(r'$x_2$', fontsize = 15)
plt.legend(loc = 1, fontsize = 12)
plt.axis('equal')
plt.xlim([0,4])
plt.ylim([0,4])
plt.show()
You might have seen the concept of entropy in physics. In fact, the concept of entropy is closely related to the cross-entropy that we just encountered in logistic regression.
Entropy in Physics In statistical physics, entropy measures the uncertainty or disorder of a system:
$$ S = -k_B \sum_{i} p_i \log p_i, $$
where:
- $ S $ is the Boltzmann entropy.
- $ k_B $ is Boltzmann's constant.
- $ p_i $ is the probability of the system being in state $ i $.
Entropy quantifies how much "randomness" or "information" is needed to describe the system's configuration.
The parallel is not a loose analogy. Shannon arrived at the same functional form independently, in the context of communication, by asking how many bits are needed on average to encode the outcome of a random variable. In information theory the entropy is written
$$ H(p) = - \sum_i p_i \log p_i , $$
which is the Boltzmann expression with $k_B$ set to one. Entropy is minimized, and equal to zero, when one outcome is certain, and it is maximized when all outcomes are equally likely. In both fields the interpretation is the same: entropy measures how much we do not know.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
p = np.linspace(0.01, 0.99, 100)
S = -p*np.log(p) - (1-p)*np.log(1-p)
plt.figure(figsize = (4, 3))
plt.plot(p, S, linewidth = 3)
plt.xlabel('p')
plt.axis('equal')
plt.grid(alpha = 0.3)
plt.show()
Cross-Entropy in Logistic Regression In logistic regression, cross-entropy measures how well the predicted probability distribution matches the true distribution:
$$ \mathcal{L} = - \sum_{i} \left[\, y_i \log \hat{y}_i + (1 - y_i) \log (1 - \hat{y}_i) \,\right], $$
where:
- $ y_i $ is the true label (1 or 0).
- $ \hat{y}_i $ is the predicted probability for class 1.
- The loss is high when the model predicts probabilities far from the true class.
Connecting the Concepts
- In physics, Boltzmann entropy measures the disorder of a system based on probabilities of different states.
- In machine learning, cross-entropy loss measures the "disorder" or "uncertainty" when the predicted probabilities do not align with the true labels.
- Improving classification performance therefore corresponds to reducing the uncertainty between the predicted and the true distributions.
Why It Is Called Cross-Entropy
The prefix "cross" indicates that two distributions are involved rather than one. Ordinary entropy $H(p) = -\sum_i p_i \log p_i$ measures the uncertainty of a single distribution $p$. Cross-entropy
$$ H(p, q) = - \sum_i p_i \log q_i $$
measures the average number of bits needed to encode outcomes drawn from the true distribution $p$ when the encoding was designed for a different distribution $q$. It is minimized, and equal to $H(p)$, exactly when $q = p$. The excess, $H(p,q) - H(p)$, is the Kullback-Leibler divergence, a measure of how far the predicted distribution is from the true one.
In our setting $p$ is the true label, a distribution placing all of its mass on one class, and $q$ is the model's predicted probability. Minimizing the cross-entropy therefore means adjusting the model until its predicted distribution matches the observed labels as closely as possible.
The Connection to Section 3.1
Note that nothing new has been introduced here. The cross-entropy is simply the negative log-likelihood derived earlier, renamed. Maximizing the likelihood of the observed labels and minimizing the cross-entropy between the true and predicted distributions are the same operation described in two vocabularies, one from statistics and one from information theory. The name matters mainly because deep learning frameworks use the second one: what PyTorch calls BCEWithLogitsLoss is the objective we derived by hand in Section 3.1.
Everything up to this point could be implemented with NumPy alone, because the models were small and their gradients could be derived by hand. That approach does not scale. A modern network has millions of parameters and a loss function whose derivative no one would attempt to write out, so a different tool is needed.
Deep learning frameworks solve three problems at once.
Automatic differentiation. The framework records the operations performed during the forward computation and applies the chain rule mechanically to obtain exact gradients. Derivatives no longer have to be derived by hand, and they are not approximated by finite differences either.
Hardware acceleration. The same code runs on a GPU, where the matrix operations that dominate training are executed in parallel across thousands of cores.
Reusable components. Layers, activation functions, loss functions, and optimizers are provided as standard building blocks, so that specifying a model becomes a matter of composition rather than implementation.
The three most widely used frameworks are introduced below. PyTorch is used for the examples in this notebook.
TensorFlow
- Platform: Linux, Mac OS, Windows
- Written in: C++, Python
- Interface: Python, C/C++, Java, Go, R
- https://www.tensorflow.org/
Keras
PyTorch is an open-source software library for deep learning.
It provides a framework for performing computation efficiently, and it can make use of the GPU (Graphics Processing Unit) to accelerate computation further. PyTorch is controlled through a simple Python API.
PyTorch is one of the most widely used libraries for implementing machine learning and other algorithms that involve a large number of mathematical operations. It was developed by Meta AI, is now governed by the PyTorch Foundation, and is among the most popular machine learning libraries on GitHub.
PyTorch gets its name from "tensors", which are arrays of arbitrary dimensionality, and from Torch, the scientific computing library it descends from.
A tensor in PyTorch is essentially a NumPy array that can live on a GPU and that remembers how it was computed. Training and inference in a neural network involve the propagation of matrix computations through many nodes of a computational graph.
The computational graph is the idea that makes automatic differentiation possible, and it is worth understanding rather than treating as an implementation detail. When an expression such as $(\omega - 4)^2$ is evaluated, PyTorch does not merely compute a number. It records the sequence of elementary operations that produced it: a subtraction, then a squaring. Each of these has a known derivative, so the chain rule can be applied backwards through the recorded sequence to obtain the derivative of the result with respect to any input.
Two consequences follow, and both matter in practice. The gradients obtained are exact, not finite-difference approximations, so no step size has to be chosen and no truncation error is incurred. And the cost of computing the gradient with respect to all parameters is comparable to the cost of a single forward evaluation, regardless of how many parameters there are. This second property is what makes training a network with millions of parameters feasible at all.
The graph in PyTorch is built as the forward computation runs, and is consumed when backward() is called. Nothing has to be declared in advance, which is why an ordinary Python for loop or if statement may appear inside a model.
Side note:
In solid mechanics, a tensor is a multi-dimensional array of numbers that describes the physical state of a material. Tensors are used to describe physical quantities like stress and strain, which have magnitude and two or more directions.
While the specific applications of tensors in PyTorch and solid mechanics differ, the underlying concept, a mathematical object capable of representing multi-dimensional relationships, remains the same in both cases.
Here we demonstrate how to implement gradient descent in PyTorch as a way of becoming familiar with the framework.
This example illustrates how PyTorch handles automatic differentiation and gradient descent updates, which is useful when building deep learning models and optimizing more complex loss functions.
The example below is deliberately trivial. The function $(\omega - 4)^2$ is a parabola whose minimum is obviously at $\omega = 4$, and no numerical method is needed to find it. That is exactly the point: because the answer is known in advance, any discrepancy in the output is attributable to the procedure rather than to the problem.
Watch for four elements in the code, since every training loop in deep learning contains the same four, no matter how complicated the model.
- The variable to be optimized is created with
requires_grad = True, which marks it as something PyTorch should track and update. Ordinary tensors are not tracked. - Calling
backward()on the loss applies the chain rule backwards through the recorded operations and deposits the result inw.grad. - The update
w -= LR*w.gradis performed insidetorch.no_grad(), so that the update itself is not recorded as part of the next graph. w.grad.zero_()clears the gradient, because PyTorch accumulates gradients rather than overwriting them. Forgetting this line is the single most common mistake in a hand-written training loop.
Observe also how the learning rate behaves. Too small, and the iterations creep toward the answer; too large, and the updates overshoot and may diverge. This parameter has no automatic setting and remains one of the most consequential choices in training any model.
$$\min_{\omega}\;\;(\omega - 4)^2$$
import numpy as np
import matplotlib.pyplot as plt
import torch
%matplotlib inline
w = torch.tensor(0.0, requires_grad = True)
LR = 0.05
# Training
cost_record = []
for i in range(50):
cost = w*w - 8*w + 16
cost.backward() # d(cost)/dw, stored in w.grad
cost_record.append(cost.item())
with torch.no_grad():
w -= LR*w.grad
w.grad.zero_() # gradients accumulate, so clear them
print('\n optimal w =', w.item())
print('\n')
plt.figure(figsize = (6, 4))
plt.plot(cost_record)
plt.xlabel('iteration')
plt.ylabel('cost')
plt.show()
We now return to the problem of Chapter 1, but solve it by a completely different route. There the optimum was obtained in one step from the normal equation $\theta = (A^TA)^{-1}A^Ty$. Here it is obtained by taking many small steps downhill on the loss surface.
For linear regression this is a step backwards in efficiency, and it should be recognized as such. Its value is that the iterative approach continues to work when the closed form does not exist, which is the situation for every model beyond this chapter. Solving a problem whose exact answer we already know is the safest way to confirm that the iterative machinery is correct, and it also allows the two answers to be compared directly.
$$\hat{y} = \omega x + b$$
- Given $x$ and $y$
- Want to estimate $\omega$ and $b$
# data generation
# data points in column vector [input, output]
train_x = np.array([0.1, 0.4, 0.7, 1.2, 1.3, 1.7, 2.2, 2.8, 3.0, 4.0, 4.3, 4.4, 4.9]).reshape(-1, 1)
train_y = np.array([0.5, 0.9, 1.1, 1.5, 1.5, 2.0, 2.2, 2.8, 2.7, 3.0, 3.5, 3.7, 3.9]).reshape(-1, 1)
m = train_x.shape[0]
plt.figure(figsize = (6, 4))
plt.plot(train_x, train_y, 'ko', alpha = 0.3)
plt.title('Data')
plt.xlabel('X')
plt.ylabel('Y')
plt.axis('equal')
plt.grid(alpha = 0.3)
plt.xlim([0, 5])
plt.show()
- Given $(x_i, y_i)$ for $i=1,\cdots, m$
$$ \hat{y}_{i} = \omega x_{i} + b \; \quad \text{ such that }\quad \min\limits_{\omega, b}\sum\limits_{i = 1}^{m} (\hat{y}_{i} - y_{i})^2$$
LR = 0.001
n_iter = 1000
train_x_t = torch.tensor(train_x, dtype = torch.float32)
train_y_t = torch.tensor(train_y, dtype = torch.float32)
w = torch.zeros((1, 1), requires_grad = True)
b = torch.zeros((1, 1), requires_grad = True)
loss_record = []
for i in range(n_iter):
cost = torch.mean((w*train_x_t + b - train_y_t)**2)
cost.backward()
loss_record.append(cost.item())
with torch.no_grad():
w -= LR*w.grad
b -= LR*b.grad
w.grad.zero_()
b.grad.zero_()
w_val = w.detach().numpy()
b_val = b.detach().numpy()
print('\n optimal w =', w_val)
print('\n optimal b =', b_val)
print('\n')
plt.figure(figsize = (6, 4))
plt.plot(loss_record)
plt.xlabel('iteration')
plt.ylabel('loss')
plt.show()
xp = np.arange(0, 5, 0.01).reshape(-1, 1)
yp = w_val*xp + b_val
plt.figure(figsize = (6, 4))
plt.plot(train_x, train_y, 'ko', alpha = 0.3)
plt.plot(xp, yp, 'r')
plt.title('Data')
plt.xlabel('X')
plt.ylabel('Y')
plt.axis('equal')
plt.grid(alpha = 0.3)
plt.xlim([0, 5])
plt.show()
This is the same model derived in Chapter 3, now expressed in PyTorch. Comparing the two implementations is instructive: the mathematics is identical, but the gradient formula derived by hand in Section 3.1 no longer appears anywhere in the code. It is obtained automatically by differentiating the loss.
This is the pattern that carries forward into deep learning. From this point on, specifying a model means stating the architecture and the loss function, and the learning rule follows from them without further human intervention.
$$ \begin{align*} \omega &= \begin{bmatrix} \omega_0 \\ \omega_1 \\ \omega_2\end{bmatrix}, \qquad x = \begin{bmatrix} 1 \\ x_1 \\ x_2\end{bmatrix}\\ \\ X &= \begin{bmatrix} \left(x^{(1)}\right)^T \\ \left(x^{(2)}\right)^T \\ \left(x^{(3)}\right)^T \\ \vdots\end{bmatrix} = \begin{bmatrix} 1 & x_1^{(1)} & x_2^{(1)} \\ 1 & x_1^{(2)} & x_2^{(2)} \\ 1 & x_1^{(3)} & x_2^{(3)} \\ \vdots & \vdots & \vdots \\\end{bmatrix}, \quad y = \begin{bmatrix} y^{(1)}\\ y^{(2)} \\y^{(3)} \\ \vdots \end{bmatrix} \end{align*} $$
# datat generation
m = 1000
true_w = np.array([[-6], [2], [1]])
train_X = np.hstack([np.ones([m,1]), 5*np.random.rand(m,1), 4*np.random.rand(m,1)])
true_w = np.asmatrix(true_w)
train_X = np.asmatrix(train_X)
train_y = 1/(1 + np.exp(-train_X*true_w)) > 0.5
C1 = np.where(train_y == True)[0]
C0 = np.where(train_y == False)[0]
train_y = np.empty([m,1])
train_y[C1] = 1
train_y[C0] = 0
plt.figure(figsize = (6, 4))
plt.plot(train_X[C1,1], train_X[C1,2], 'ro', alpha = 0.3, label = 'C1')
plt.plot(train_X[C0,1], train_X[C0,2], 'bo', alpha = 0.3, label = 'C0')
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.legend(loc = 1)
plt.axis('equal')
plt.ylim([0,4])
plt.show()
$$ \begin{align*} \ell(\omega) = \log \mathscr{L}(\omega) &= \sum_{i=1}^{m} y^{(i)} \log h_{\omega} \left(x^{(i)} \right) + \left(1-y^{(i)} \right) \log \left(1-h_{\omega} \left(x^{(i)} \right) \right)\\\\ &\Rightarrow \frac{1}{m} \sum_{i=1}^{m} y^{(i)} \log h_{\omega} \left(x^{(i)} \right) + \left(1-y^{(i)} \right) \log \left(1-h_{\omega} \left(x^{(i)} \right) \right) \end{align*} $$
In practice, optimization problems in machine learning are commonly formulated as the minimization of an objective function (or loss function) rather than the maximization.
Two changes are made in the expressions below, and each has a practical reason.
Division by $m$. Replacing the sum by an average makes the loss independent of the number of training examples. Without it, a suitable learning rate for a dataset of $100$ points would be far too large for one of $10{,}000$ points, since the gradient would scale with the dataset size.
Negation. Optimizers in every deep learning framework are written to minimize, so a maximization problem is converted by flipping the sign. Maximizing the log-likelihood and minimizing the negative log-likelihood, which is the cross-entropy of Section 3.3, are the same problem.
$$ \begin{align*} \hat{\omega} &= \arg\max_{\omega} \ell(\omega)\\ & = \arg \min_{\omega} -\ell(\omega) \end{align*} $$
LR = 0.1
n_iter = 3000
train_x_t = torch.tensor(np.asarray(train_X), dtype = torch.float32)
train_y_t = torch.tensor(np.asarray(train_y), dtype = torch.float32)
w = torch.zeros((3, 1), requires_grad = True)
loss_record = []
for i in range(n_iter):
y_pred = torch.sigmoid(torch.matmul(train_x_t, w))
loss = - train_y_t*torch.log(y_pred) - (1 - train_y_t)*torch.log(1 - y_pred)
loss = torch.mean(loss)
loss.backward()
loss_record.append(loss.item())
with torch.no_grad():
w -= LR*w.grad
w.grad.zero_()
w_hat = w.detach().numpy()
print(w_hat)
print('\n')
plt.figure(figsize = (6, 4))
plt.plot(loss_record)
plt.xlabel('iteration')
plt.ylabel('loss')
plt.show()
xp = np.arange(0, 4, 0.01).reshape(-1, 1)
yp = - w_hat[1,0]/w_hat[2,0]*xp - w_hat[0,0]/w_hat[2,0]
plt.figure(figsize = (6, 4))
plt.plot(train_X[C1,1], train_X[C1,2], 'ro', alpha = 0.3, label = 'C1')
plt.plot(train_X[C0,1], train_X[C0,2], 'bo', alpha = 0.3, label = 'C0')
plt.plot(xp, yp, 'g', linewidth = 3, label = 'Logistic Regression')
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.legend(loc = 1)
plt.axis('equal')
plt.ylim([0,4])
plt.show()
Instead of manually defining the cross-entropy loss function, PyTorch's built-in loss functions can be utilized for greater convenience and efficiency.
The reason for preferring the built-in function is numerical, not merely a matter of convenience, and it is worth understanding because the same issue arises in every framework.
Computing the sigmoid first and then taking its logarithm invites disaster at the extremes. For a large negative input the sigmoid underflows to exactly $0$ in floating point, and $\log 0$ is $-\infty$, which propagates through the rest of the computation as nan. The built-in functions avoid forming the intermediate probability at all, combining the sigmoid and the logarithm into a single algebraically equivalent expression that remains stable for arguments of any magnitude.
This is the origin of the term logits appearing in the function names. A logit is the raw score $\omega^T x$ before the sigmoid is applied. These functions expect the logits as input, not the probabilities, and passing an already-squashed probability into them is a common and quiet error: the code runs, but the sigmoid is applied twice and the model trains poorly for no visible reason.
PyTorch embedded functions
- torch.nn.BCEWithLogitsLoss for binary classification
- torch.nn.CrossEntropyLoss for multiclass classification
LR = 0.1
n_iter = 3000
train_x_t = torch.tensor(np.asarray(train_X), dtype = torch.float32)
train_y_t = torch.tensor(np.asarray(train_y), dtype = torch.float32)
w = torch.zeros((3, 1), requires_grad = True)
loss_fn = torch.nn.BCEWithLogitsLoss() # sigmoid and log in one stable step
loss_record = []
for i in range(n_iter):
y_pred = torch.matmul(train_x_t, w) # the logits, before the sigmoid
loss = loss_fn(y_pred, train_y_t)
loss.backward()
loss_record.append(loss.item())
with torch.no_grad():
w -= LR*w.grad
w.grad.zero_()
w_hat = w.detach().numpy()
print(w_hat)
print('\n')
xp = np.arange(0, 4, 0.01).reshape(-1, 1)
yp = - w_hat[1,0]/w_hat[2,0]*xp - w_hat[0,0]/w_hat[2,0]
plt.figure(figsize = (6, 4))
plt.plot(train_X[C1,1], train_X[C1,2], 'ro', alpha = 0.3, label = 'C1')
plt.plot(train_X[C0,1], train_X[C0,2], 'bo', alpha = 0.3, label = 'C0')
plt.plot(xp, yp, 'g', linewidth = 3, label = 'Logistic Regression')
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.legend(loc = 1)
plt.axis('equal')
plt.ylim([0,4])
plt.show()
Key takeaways
- Every method in this notebook is made of the same three parts: a model $f(x\,;\theta)$, a loss that scores how badly a given $\theta$ explains the data, and an optimization algorithm that reduces it. Learning means nothing until the second part has been committed to.
- Linear regression is least squares. Written as $\lVert A\theta - y \rVert^2$ it has the closed-form solution $\theta = (A^TA)^{-1}A^Ty$, which is exact and available in one step — a luxury that disappears for every later model.
- The perceptron classifies by the sign of $\omega^Tx$ and updates only on the points it gets wrong. It converges whenever the data are linearly separable, but it stops at the first separating hyperplane it finds and has nothing to say when no such hyperplane exists.
- Logistic regression replaces the hard sign by the sigmoid, so that every point contributes according to its distance from the boundary. Its output is a probability, its loss is the cross-entropy, and minimizing that loss is maximum likelihood.
- A deep learning framework supplies automatic differentiation, GPU execution, and standard building blocks. Automatic differentiation is the essential one: it records the forward computation and applies the chain rule backwards, so gradients are exact rather than finite-difference approximations, and cost about as much as one forward pass no matter how many parameters there are.
- In PyTorch a training loop is four lines of bookkeeping — mark the parameters with
requires_grad, callbackward()on the loss, update undertorch.no_grad(), and clear the gradients withzero_(), since PyTorch accumulates them. Solving linear and logistic regression this way gives the same answers obtained analytically, which is the point: the iterative machinery is verified on problems whose answers are already known. - Built-in losses such as
BCEWithLogitsLosstake the logits, not the probabilities, and combine the sigmoid with the logarithm into one numerically stable expression. Passing an already-squashed probability into them is a quiet and common error.
Looking ahead. From this point on, specifying a model means stating an architecture and a loss function; the learning rule follows from them automatically. That is what makes it possible to leave linear boundaries behind. The next notebook stacks the logistic units of Chapter 3 into layers, so that the features a linear model must be given by hand are instead learned from the data, and shows how the same training loop written here scales to a network with thousands of parameters.