Linear Models and Gradient Descent from Scratch

01A line, a loss, and what the loss assumes

Linear regression predicts a number as a weighted sum of features plus a bias: ŷ = wᵀx + b. Stacking the n examples as rows of a design matrix X turns all predictions into one matrix–vector product Xw.^[1] The fit is scored by the mean squared error, and the choice of that loss is not arbitrary: it is the negative log-likelihood of a model where the target is the linear prediction plus Gaussian noise.^[1]

The quantity inside the parentheses is the residual. Squaring it penalizes large misses far more than small ones, which is both the strength of the loss and its sensitivity to outliers.^[1]

Because the loss is a quadratic bowl, calculus finds its bottom directly. Setting the gradient to zero gives the normal equations and the solution w* = (XᵀX)⁻¹Xᵀy, unique whenever the columns of X are linearly independent.^[1][6] Solving them directly becomes numerically fragile when features are nearly collinear: XᵀX is then close to singular and the fitted weights can take large magnitudes. Adding a regularization term, as in §5, keeps the matrix non-singular.^[6]

02Gradient descent and the learning rate

Few models past linear regression have a closed-form optimum, so the method that matters is gradient descent: evaluate the gradient of the loss at the current parameters and step against it, scaled by a learning rate η.^[1]

For the squared loss, ∇L = (2/n) Xᵀ(Xw − y): the residuals pushed back through the features. The gradient of the logistic loss in §4 has the same shape.^[6]

The learning rate has a ceiling. In one dimension a quadratic loss ½λw² has gradient λw, so one step maps w to (1 − ηλ)w. That shrinks only when |1 − ηλ| < 1, which means η < 2/λ; past it the iterate overshoots and grows without bound.^[2] In several dimensions the same test applies along each eigenvector of the loss's Hessian, so the step must be small enough for the steepest direction, and the flattest direction then crawls. The ratio of the two curvatures is the condition number, and it is set by the data: a length measured in millimetres next to one in kilometres gives curvatures about 10¹² apart, a bowl whose axes differ a million to one. That is why inputs are standardized to zero mean and unit variance before fitting, and why preconditioning, which in its diagonal form is a per-coordinate learning rate, improves plain descent.^[2][8]

03Minibatches trade noise for speed

The full gradient averages over every training example, which costs a pass over the dataset per step. Stochastic gradient descent uses a random minibatch instead. The minibatch gradient is an unbiased estimate of the full one, so on average the steps point the right way; what changes is the noise around that average.^[3]

The noise follows the §2 rule for averages on the probability page: the standard deviation of a B-example gradient falls like 1/√B. A batch of 4 is half as noisy as a batch of 1 and costs four times as much per step; the batch sizes used in practice, often 32 to 256, balance that against memory and hardware throughput.^[1] The noise does not vanish as training proceeds, so a constant learning rate leaves the iterate jittering around the optimum, and the convergence proofs for convex losses shrink the learning rate as training proceeds.^[3]

04Logistic regression: the same gradient, a different likelihood

For a yes-or-no label, the linear score z = wᵀx + b is passed through the sigmoid to produce a probability. Maximum likelihood for a Bernoulli label gives the cross-entropy loss from the probability page, and its gradient is (p − y)x: the same error-times-input form as linear regression, because the derivative of the sigmoid cancels against the derivative of the log.^[6] That is logistic regression, a classifier despite the name.

The decision boundary is where z = 0. Points far from it on the correct side contribute almost nothing to the gradient, because p − y is near zero for them; misclassified points and points near the boundary carry most of it.

05Weight decay and overfitting

A model that fits its training data perfectly can predict new data badly. The gap between training error and error on held-out data is overfitting; the measure that matters is its error on data it did not train on, which is why a validation set exists.^[5] A linear model can overfit when it has too many features for its data, including polynomial features built from one input.

Weight decay adds (λ/2)‖w‖² to the loss. The gradient of the penalty is λw, so each descent step first shrinks the weights by (1 − ηλ) and then applies the data gradient; λ dials the model's effective complexity continuously.^[4] In closed form it is ridge regression, solving (XᵀX + αI)w = Xᵀy; with the mean squared error as the data term the two penalties match when α = nλ/2. The added diagonal also fixes the collinearity problem of §1.^[7][6] The bias is conventionally left unpenalized.^[4]

06Checking a gradient numerically

An analytic gradient with a wrong sign or a dropped factor still produces a loss that goes down for a while, which makes the bug hard to see. The check is a finite difference: perturb one parameter by ±h, evaluate the loss twice, and compare (L(w + h) − L(w − h)) / 2h with the analytic component. The centered form is more accurate than a one-sided difference, and h around 10⁻⁵ is a common choice.^[9]

07Both models, checked against closed forms

These complete programs fit a line by gradient descent and compare it with the closed-form least-squares solution, then fit a logistic classifier by gradient descent and verify its analytic gradient against centered finite differences. The names match between the two listings.

#include <cassert>
#include <cmath>
#include <cstddef>
#include <vector>
struct Example { double feature; double target; };
// Mean squared error of y ≈ weight * feature + bias, and its gradient.
struct Fit { double loss, gradientWeight, gradientBias; };
Fit squaredLoss(const std::vector<Example>& data, double weight, double bias) {
    Fit fit{0, 0, 0};
    for (const Example& example : data) {
        const double residual = weight * example.feature + bias - example.target;
        fit.loss += residual * residual / data.size();
        fit.gradientWeight += 2 * residual * example.feature / data.size(); // Residual pushed back through the feature.
        fit.gradientBias += 2 * residual / data.size();
    }
    return fit;
}
double sigmoid(double logit) { return logit >= 0 ? 1 / (1 + std::exp(-logit)) : std::exp(logit) / (1 + std::exp(logit)); }
// Mean cross-entropy of a Bernoulli label under p = sigmoid(weight * feature + bias).
Fit logisticLoss(const std::vector<Example>& data, double weight, double bias) {
    Fit fit{0, 0, 0};
    for (const Example& example : data) {
        const double logit = weight * example.feature + bias, probability = sigmoid(logit);
        // log(1 + e^-|z|) + max(0, ∓z) equals −log p or −log(1 − p) without overflow.
        const double softplus = std::log1p(std::exp(-std::abs(logit)));
        fit.loss += (example.target == 1 ? softplus + std::fmax(0.0, -logit) : softplus + std::fmax(0.0, logit)) / data.size();
        fit.gradientWeight += (probability - example.target) * example.feature / data.size(); // (p − y) x
        fit.gradientBias += (probability - example.target) / data.size();
    }
    return fit;
}
int main() {
    std::vector<Example> regression;
    for (int index = 0; index < 20; ++index) {
        const double feature = -2.0 + 0.2 * index;
        regression.push_back({feature, 1.5 * feature + 0.5 + 0.3 * std::sin(7.0 * index)}); // Deterministic "noise".
    }
    // Closed form for one feature: weight = cov(x, y) / var(x), bias = mean(y) − weight * mean(x).
    double meanFeature = 0, meanTarget = 0;
    for (const Example& example : regression) { meanFeature += example.feature / 20; meanTarget += example.target / 20; }
    double covariance = 0, variance = 0;
    for (const Example& example : regression) {
        covariance += (example.feature - meanFeature) * (example.target - meanTarget);
        variance += (example.feature - meanFeature) * (example.feature - meanFeature);
    }
    const double closedWeight = covariance / variance, closedBias = meanTarget - closedWeight * meanFeature;
    double weight = 0, bias = 0;
    for (int step = 0; step < 2000; ++step) { // Plain gradient descent at a safe learning rate.
        const Fit fit = squaredLoss(regression, weight, bias);
        weight -= 0.1 * fit.gradientWeight; bias -= 0.1 * fit.gradientBias;
    }
    assert(std::abs(weight - closedWeight) < 1e-6 && std::abs(bias - closedBias) < 1e-6); // Descent reaches the normal-equations optimum.
    std::vector<Example> classification;
    for (int index = 0; index < 20; ++index) {
        const double feature = -3.0 + 0.3 * index + 0.8 * std::sin(5.0 * index);
        classification.push_back({feature, feature + 0.6 * std::cos(3.0 * index) > 0 ? 1.0 : 0.0}); // Overlapping classes.
    }
    const double probeWeight = 0.7, probeBias = -0.2, step = 1e-5;
    const Fit analytic = logisticLoss(classification, probeWeight, probeBias);
    const double numericalWeight = (logisticLoss(classification, probeWeight + step, probeBias).loss - logisticLoss(classification, probeWeight - step, probeBias).loss) / (2 * step);
    const double numericalBias = (logisticLoss(classification, probeWeight, probeBias + step).loss - logisticLoss(classification, probeWeight, probeBias - step).loss) / (2 * step);
    assert(std::abs(numericalWeight - analytic.gradientWeight) < 1e-8 && std::abs(numericalBias - analytic.gradientBias) < 1e-8); // Gradient check.
    double classifierWeight = 0, classifierBias = 0;
    const double initialLoss = logisticLoss(classification, 0, 0).loss;
    for (int iteration = 0; iteration < 500; ++iteration) {
        const Fit fit = logisticLoss(classification, classifierWeight, classifierBias);
        classifierWeight -= 0.5 * fit.gradientWeight; classifierBias -= 0.5 * fit.gradientBias;
    }
    assert(std::abs(initialLoss - std::log(2.0)) < 1e-12);                               // All-zero weights predict 0.5: loss log 2.
    assert(logisticLoss(classification, classifierWeight, classifierBias).loss < 0.5 * initialLoss); // Training halved it.
}

struct Example { feature: f64, target: f64 }
// Mean squared error of y ≈ weight * feature + bias, and its gradient.
struct Fit { loss: f64, gradient_weight: f64, gradient_bias: f64 }
fn squared_loss(data: &[Example], weight: f64, bias: f64) -> Fit {
    let mut fit = Fit { loss: 0.0, gradient_weight: 0.0, gradient_bias: 0.0 };
    let count = data.len() as f64;
    for example in data {
        let residual = weight * example.feature + bias - example.target;
        fit.loss += residual * residual / count;
        fit.gradient_weight += 2.0 * residual * example.feature / count; // Residual pushed back through the feature.
        fit.gradient_bias += 2.0 * residual / count;
    }
    fit
}
fn sigmoid(logit: f64) -> f64 { if logit >= 0.0 { 1.0 / (1.0 + (-logit).exp()) } else { logit.exp() / (1.0 + logit.exp()) } }
// Mean cross-entropy of a Bernoulli label under p = sigmoid(weight * feature + bias).
fn logistic_loss(data: &[Example], weight: f64, bias: f64) -> Fit {
    let mut fit = Fit { loss: 0.0, gradient_weight: 0.0, gradient_bias: 0.0 };
    let count = data.len() as f64;
    for example in data {
        let logit = weight * example.feature + bias;
        let probability = sigmoid(logit);
        // log(1 + e^-|z|) + max(0, ∓z) equals −log p or −log(1 − p) without overflow.
        let softplus = (-logit.abs()).exp().ln_1p();
        fit.loss += if example.target == 1.0 { softplus + (-logit).max(0.0) } else { softplus + logit.max(0.0) } / count;
        fit.gradient_weight += (probability - example.target) * example.feature / count; // (p − y) x
        fit.gradient_bias += (probability - example.target) / count;
    }
    fit
}
fn main() {
    let regression: Vec<Example> = (0..20).map(|index| {
        let feature = -2.0 + 0.2 * index as f64;
        Example { feature, target: 1.5 * feature + 0.5 + 0.3 * (7.0 * index as f64).sin() } // Deterministic "noise".
    }).collect();
    // Closed form for one feature: weight = cov(x, y) / var(x), bias = mean(y) − weight * mean(x).
    let mean_feature = regression.iter().map(|example| example.feature).sum::<f64>() / 20.0;
    let mean_target = regression.iter().map(|example| example.target).sum::<f64>() / 20.0;
    let covariance: f64 = regression.iter().map(|example| (example.feature - mean_feature) * (example.target - mean_target)).sum();
    let variance: f64 = regression.iter().map(|example| (example.feature - mean_feature).powi(2)).sum();
    let closed_weight = covariance / variance;
    let closed_bias = mean_target - closed_weight * mean_feature;
    let (mut weight, mut bias) = (0.0, 0.0);
    for _ in 0..2000 { // Plain gradient descent at a safe learning rate.
        let fit = squared_loss(&regression, weight, bias);
        weight -= 0.1 * fit.gradient_weight;
        bias -= 0.1 * fit.gradient_bias;
    }
    assert!((weight - closed_weight).abs() < 1e-6 && (bias - closed_bias).abs() < 1e-6); // Descent reaches the normal-equations optimum.
    let classification: Vec<Example> = (0..20).map(|index| {
        let feature = -3.0 + 0.3 * index as f64 + 0.8 * (5.0 * index as f64).sin();
        Example { feature, target: if feature + 0.6 * (3.0 * index as f64).cos() > 0.0 { 1.0 } else { 0.0 } } // Overlapping classes.
    }).collect();
    let (probe_weight, probe_bias, step) = (0.7, -0.2, 1e-5);
    let analytic = logistic_loss(&classification, probe_weight, probe_bias);
    let numerical_weight = (logistic_loss(&classification, probe_weight + step, probe_bias).loss - logistic_loss(&classification, probe_weight - step, probe_bias).loss) / (2.0 * step);
    let numerical_bias = (logistic_loss(&classification, probe_weight, probe_bias + step).loss - logistic_loss(&classification, probe_weight, probe_bias - step).loss) / (2.0 * step);
    assert!((numerical_weight - analytic.gradient_weight).abs() < 1e-8 && (numerical_bias - analytic.gradient_bias).abs() < 1e-8); // Gradient check.
    let (mut classifier_weight, mut classifier_bias) = (0.0, 0.0);
    let initial_loss = logistic_loss(&classification, 0.0, 0.0).loss;
    for _ in 0..500 {
        let fit = logistic_loss(&classification, classifier_weight, classifier_bias);
        classifier_weight -= 0.5 * fit.gradient_weight;
        classifier_bias -= 0.5 * fit.gradient_bias;
    }
    assert!((initial_loss - 2.0_f64.ln()).abs() < 1e-12);                                      // All-zero weights predict 0.5: loss log 2.
    assert!(logistic_loss(&classification, classifier_weight, classifier_bias).loss < 0.5 * initial_loss); // Training halved it.
}

08Where linear models go wrong

Report the loss on held-out data next to the training loss, and treat a widening gap between the two as the signal to regularize or simplify rather than to train longer.^[5]

Weights on collinear features are not interpretable individually. If two features are nearly proportional, the normal equations can trade weight between them almost freely, and small changes in the data swing both; only their combined effect is pinned down.^[6]

A learning rate that works on standardized inputs can diverge on raw ones, because the Hessian scales with the square of the feature scale. When a loss goes to infinity on the first steps, check the input scaling before anything else.

09What's next

Neural Networks keeps every piece of this loop and replaces the linear model with a stack of linear layers separated by nonlinearities. The gradient is no longer one line of algebra, which is why that page derives backpropagation, and why the page after it builds the machinery that computes any gradient automatically.

10Sources

The cited texts define the models and the update rules. Widget readouts report values computed in the page from the seeded datasets.

1. Aston Zhang, Zachary C. Lipton, Mu Li, Alexander J. Smola, 2023. Dive into Deep Learning, Section 3.1: Linear Regression, Cambridge University Press. The model and design matrix, squared loss, the analytic solution and its uniqueness condition, minibatch SGD and typical batch sizes, and squared loss as Gaussian maximum likelihood.
2. Aston Zhang et al., 2023. Dive into Deep Learning, Section 12.3: Gradient Descent. Learning rates too small or too large, divergence by overshoot, and preconditioning for mismatched feature scales.
3. Aston Zhang et al., 2023. Dive into Deep Learning, Section 12.4: Stochastic Gradient Descent. The stochastic gradient as an unbiased estimate, noise around the optimum, and decaying learning rates for convergence.
4. Aston Zhang et al., 2023. Dive into Deep Learning, Section 3.7: Weight Decay. The (λ/2)‖w‖² penalty, the (1 − ηλ) shrink in the update, and the unpenalized bias.
5. Aston Zhang et al., 2023. Dive into Deep Learning, Section 3.6: Generalization. Training versus generalization error, overfitting and underfitting, and the validation set.
6. Christopher M. Bishop, 2006. Pattern Recognition and Machine Learning, Springer. The normal equations (3.15), numerical difficulty with near-collinear features and regularization as its remedy (Section 3.1.1), the sequential update (3.22), the logistic cross-entropy error (4.90), its gradient (4.91), and overfitting of maximum likelihood on separable data (Section 4.3.2).
7. scikit-learn developers. Linear Models. The Ridge objective ‖Xw − y‖² + α‖w‖² and the regularized logistic regression objective.
8. scikit-learn developers. Preprocessing data. Standardization to zero mean and unit variance, and that linear models benefit from it.
9. Stanford CS231n course staff. Optimization: Stochastic Gradient Descent. The numerical gradient, the centered-difference formula, and h ≈ 10⁻⁵ as a working step size.