Logistic Regression from Scratch

Build a logistic regression classifier with PCA, regularization, and k-fold cross-validation in NumPy.

Introduction

30 Second Summary

Two lines of code. That's all sklearn needs to train an entire machine learning model:

from sklearn.linear_model import LogisticRegression

model = LogisticRegression().fit(X, y)

But what's actually happening inside that .fit() call? Gradient descent, loss functions, regularization, matrix math. In this project, you crack open the black box and build every piece yourself with nothing but NumPy.

What You'll Build

Picture running your script, watching it reduce 8 noisy features down to 2 principal components, train a regularized classifier with proper cross-validation, then print precision and recall scores that prove your model handles an imbalanced dataset correctly.

By the end of this project, you'll have:

  • A PCA implementation that computes eigenvectors from a covariance matrix and reduces 8 features to 2 principal components you can visualize.
  • A numerically stable logistic regression with L2 regularization, where you deliberately trigger NaN from naive math, then fix it with the log-sum-exp trick.
  • A k-fold cross-validation pipeline that properly separates training and validation data, tunes hyperparameters, and reports precision, recall, and F1 on a held-out test set.
  • Secret Mission: Skip PCA entirely and watch the curse of dimensionality cause overfitting, proving why dimensionality reduction matters.

Are there any prerequisites?

You should be comfortable with Python basics (functions, loops, arrays) and have a rough memory of what a derivative is from calculus. You'll also need NumPy and Matplotlib installed already.

Before We Start

Before We Start...

Before we dive in, let's lock in what you're about to build and why it matters.

Set Up Your Project

You're about to build a logistic regression classifier from raw math. Every function you write will trace directly back to a formula you understand.

To get started, you need a clean Python file with the right libraries imported. In this step, you'll create the project file and confirm that NumPy and Matplotlib load without errors.

In this step, get ready to:
  • Create your project file on the Desktop.
  • Import NumPy and Matplotlib.
  • Verify both libraries load correctly.
Create your project file

You need a single Python file that will hold your entire logistic regression implementation. You'll save it on the Desktop so it's easy to find and run.

  • In your code editor, press Cmd+N (macOS) or Ctrl+N (Windows) to create a new file.
  • Press Cmd+S (macOS) or Ctrl+S (Windows) to open the save dialog.
  • Navigate to your Desktop, name the file logistic_regression.py, and click Save.
Add the setup code
  • Add the following code to logistic_regression.py to import both libraries and print their versions:
import numpy as np
import matplotlib.pyplot as plt

# --- Setup verification ---
print(f"Setup complete. NumPy version: {np.__version__}, Matplotlib version: {plt.matplotlib.__version__}")

What does this code do?

  • import numpy as np loads NumPy and gives it the short alias np. You'll use this for all your math operations.
  • import matplotlib.pyplot as plt loads the plotting module so you can visualize data and training progress.
  • The print statement confirms both libraries are installed and accessible by printing their version numbers.
  • Save the file with Cmd+S (macOS) or Ctrl+S (Windows).

Now you'll run the script to confirm everything works.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type Terminal and press Enter to open it.
  • Move to your Desktop by running this command:
cd ~/Desktop

Before you run the script, what do you expect to see printed in the terminal?

  • Run the script by entering this command:
python logistic_regression.py

You should see output like this:

Setup complete. NumPy version: 2.5.1, Matplotlib version: 3.11.1

What should I see?

Your exact version numbers might differ slightly (for example 2.5.0 instead of 2.5.1). The important thing is that the message prints without any errors. If it prints, both libraries loaded successfully.

🙋‍♀️ Seeing ModuleNotFoundError?

  • Make sure you're running the same Python installation where you installed NumPy and Matplotlib. Try python3 linear_regression.py if python alone doesn't work.
  • Check that your terminal's working directory is the Desktop. Run ls and confirm you see linear_regression.py in the output.
  • If the module is genuinely missing, install it with pip install numpy matplotlib and try again.

Help me fix a ModuleNotFoundError when running my linear_regression.py script.

✔️ Awesome, I've got everything!

Great. Double check you've saved linear_regression.py and move on to the next step.

ⓧ I'd like to double check the full code

Here's the complete linear_regression.py so far:

import numpy as np
import matplotlib.pyplot as plt

# --- Setup verification ---
print(f"Setup complete. NumPy version: {np.__version__}, Matplotlib version: {plt.matplotlib.__version__}")

Your environment is confirmed working. Next up, you'll learn why the full ML pipeline matters before diving into code.

Generate Data and Understand the Pipeline

Your environment works. Now you need data to classify, and a mental model of how a real ML pipeline fits together.

In production, you never jump straight to training. You explore the data, engineer better features, check for class imbalance, and split carefully. Skipping any of these steps is how models fail silently in the real world.

In this step, get ready to:
  • Generate an imbalanced synthetic dataset with correlated features.
  • Understand why EDA, feature engineering, and class imbalance matter.
  • Learn why accuracy alone is a trap for imbalanced data.
Generate a synthetic imbalanced dataset

You will create a 2-class dataset with 8 features. Some features are correlated (making PCA useful later), and the class split is 70/30 (making accuracy misleading).

  • Add the data generation code to logistic_regression.py below the setup print statement:
# --- Generate imbalanced dataset with correlated features ---
np.random.seed(42)
n_samples = 300

# Class 0 (majority: 210 samples) and Class 1 (minority: 90 samples)
n_class0 = 210
n_class1 = 90

# Generate 4 base features per class with different means
X_class0 = np.random.randn(n_class0, 4) + np.array([0, 0, 0, 0])
X_class1 = np.random.randn(n_class1, 4) + np.array([1.5, 1.2, 0.8, 1.0])

# Stack into one dataset
X_base = np.vstack([X_class0, X_class1])
y = np.array([0] * n_class0 + [1] * n_class1)

# Create 4 correlated features (linear combos + noise)
X_corr = np.column_stack([
    X_base[:, 0] * 0.8 + X_base[:, 1] * 0.3 + np.random.randn(n_samples) * 0.2,
    X_base[:, 1] * 0.7 + X_base[:, 2] * 0.4 + np.random.randn(n_samples) * 0.2,
    X_base[:, 2] * 0.6 + X_base[:, 3] * 0.5 + np.random.randn(n_samples) * 0.3,
    np.random.randn(n_samples) * 2  # Pure noise feature
])

# Final dataset: 8 features total
X = np.hstack([X_base, X_corr])
print(f"X shape: {X.shape}, y shape: {y.shape}")
print(f"Class distribution: {np.sum(y==0)} negative, {np.sum(y==1)} positive")

What does this code do?

  • It creates 4 base features with different means per class. The class 1 samples are shifted, giving the classifier a signal to learn.
  • It adds 4 more features that are linear combinations of the originals. These correlated features are exactly what PCA will compress later.
  • The 70/30 class split means a model that always predicts class 0 gets 70% accuracy while catching zero positive cases.
  • Save logistic_regression.py and run it:
python logistic_regression.py

You should see X shape: (300, 8) and Class distribution: 210 negative, 90 positive printed.

Why the full pipeline matters

What is EDA and why does it matter?

Exploratory Data Analysis means looking at distributions, correlations, and outliers before modelling. A 5-minute histogram can reveal that one feature has 50% missing values, saving you hours of debugging a broken model.

In our dataset, you would check feature distributions per class and notice that features 5-7 are highly correlated with features 1-3. That correlation is redundant information that PCA will compress.

💡 What is feature engineering?

Feature engineering means creating better inputs from raw data. Combining a user's login count and session length into an "engagement score" gives the model a single informative feature instead of two noisy ones.

Our correlated features are an example of poor raw features. PCA is one form of automated feature engineering: it creates new axes (principal components) that capture the maximum variance in fewer dimensions.

💡 Why is class imbalance dangerous?

With 210 negatives and 90 positives, a model that always predicts "negative" scores 70% accuracy. That looks good on paper but catches zero fraud cases. This is why precision, recall, and F1 score exist.

Precision asks: of everything you flagged positive, how many actually were? Recall asks: of all real positives, how many did you catch? F1 balances both into a single number.

💡 What are boolean masks and why use them?

A boolean mask is a True/False array that selects specific rows. When you write X[y == 1], NumPy returns only the rows where y is 1. You will use masks to split data into folds, separate classes, and filter predictions throughout this project.

You have data and context. Next up, you will reduce those 8 features down to their essential components using PCA.

Implement PCA from Scratch

You have 8 features, but several are just noisy copies of each other. Training on redundant information slows convergence and adds noise.

Principal Component Analysis finds the directions of maximum variance in your data and projects it onto fewer dimensions. You will implement it from the covariance matrix and eigendecomposition, keeping only the top-k eigenvectors.

In this step, get ready to:
  • Standardize the features to zero mean and unit variance.
  • Compute the covariance matrix and its eigenvectors.
  • Project the data onto the top 2 principal components and visualize it.
Standardize the features

PCA is sensitive to scale. A feature measured in thousands will dominate one measured in decimals. Standardizing puts all features on the same footing.

  • Add the standardization code below your data generation in logistic_regression.py:
# --- PCA from scratch ---
# Step 1: Standardize features (zero mean, unit variance)
X_mean = np.mean(X, axis=0)
X_std = np.std(X, axis=0)
X_standardized = (X - X_mean) / X_std

print(f"Mean after standardization: {X_standardized.mean(axis=0).round(2)}")
print(f"Std after standardization: {X_standardized.std(axis=0).round(2)}")

What does this code do?

  • X_mean computes the mean of each column (feature). Subtracting it centers the data at zero.
  • Dividing by X_std scales each feature to unit variance. Now no feature dominates due to scale alone.
  • The print statements confirm all means are ~0 and all standard deviations are ~1.
  • Save logistic_regression.py and run it:
python logistic_regression.py

You should see arrays of zeros for the mean and ones for the standard deviation (or very close to them due to floating point).

Compute eigenvectors and eigenvalues

The covariance matrix captures how features move together. Its eigenvectors point in the directions of maximum variance, and eigenvalues tell you how much variance each direction explains.

  • Add the eigendecomposition code below your standardization:
# Step 2: Compute covariance matrix and eigendecomposition
cov_matrix = np.cov(X_standardized, rowvar=False)
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)

# Sort by eigenvalue descending (largest variance first)
sorted_indices = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[sorted_indices]
eigenvectors = eigenvectors[:, sorted_indices]

# Show how much variance each component explains
explained_variance_ratio = eigenvalues / eigenvalues.sum()
print(f"Explained variance per component: {explained_variance_ratio.round(3)}")
print(f"Top 2 components explain: {explained_variance_ratio[:2].sum():.1%} of total variance")

What does this code do?

  • np.cov(X_standardized, rowvar=False) computes the 8x8 covariance matrix. Each entry (i, j) measures how features i and j co-vary.
  • np.linalg.eigh returns eigenvalues and eigenvectors for a symmetric matrix. Eigenvectors are the principal component directions.
  • Sorting by eigenvalue descending ensures the first component captures the most variance. The top 2 components should explain most of it because our features were built from just 4 independent sources.
  • Save and run:
python logistic_regression.py

You should see that the top 2 components explain roughly 50-70% of the total variance. The correlated features collapse into fewer components.

Project onto top-k components and visualize
  • Add the projection and plotting code:
# Step 3: Project data onto top 2 principal components
k = 2  # Number of components to keep
W_pca = eigenvectors[:, :k]  # Top-k eigenvectors as columns
X_pca = X_standardized @ W_pca  # Project: (300, 8) @ (8, 2) = (300, 2)

print(f"Reduced X shape: {X_pca.shape}")

# Visualize the 2D projection colored by class
plt.figure(figsize=(8, 6))
plt.scatter(X_pca[y == 0, 0], X_pca[y == 0, 1], alpha=0.5, label="Class 0")
plt.scatter(X_pca[y == 1, 0], X_pca[y == 1, 1], alpha=0.5, label="Class 1")
plt.xlabel("Principal Component 1")
plt.ylabel("Principal Component 2")
plt.title("PCA: 8 Features Reduced to 2 Components")
plt.legend()
plt.show()

What does the projection do?

  • W_pca is shape (8, 2). Each column is an eigenvector pointing in a high-variance direction.
  • X_standardized @ W_pca projects each 8-dimensional sample onto those 2 directions. The result is a 2D representation that preserves the most variance.
  • The boolean mask y == 0 selects only class 0 rows for plotting. This is how you separate classes for visualization without loops.
  • Save logistic_regression.py and run:
python logistic_regression.py

A scatter plot opens showing two partially overlapping clusters. The classes are separable but not perfectly. Close the plot window to continue.

Clusters completely overlapping?

  • Make sure you sorted eigenvalues in descending order with [::-1]. Without sorting, you might project onto the lowest-variance directions.
  • Confirm rowvar=False in np.cov. Without it, NumPy treats rows as variables instead of samples.

Help me debug my PCA projection.

You compressed 8 features into 2 while keeping most of the signal. Next, you will build the logistic regression classifier that learns to separate these clusters.

Build Logistic Regression with Numerical Stability

Your data is now 2D and ready for classification. The next challenge is building the model that draws a decision boundary between the clusters.

Logistic regression uses a sigmoid function to squash predictions between 0 and 1. The problem is that naive implementations hit log(0) and produce NaN. You will deliberately trigger this, then fix it with the log-sum-exp trick.

In this step, get ready to:
  • Implement sigmoid and binary cross-entropy loss (the naive way that breaks).
  • Add L2 regularization to prevent overfitting.
  • Fix the numerical instability with a stable implementation.
Implement the naive sigmoid and loss

The sigmoid function maps any real number to (0, 1). Binary cross-entropy measures how wrong those predictions are. Together they form the logistic regression loss.

  • Add the naive implementation below your PCA code in logistic_regression.py:
# --- Logistic Regression ---
# Naive sigmoid (will cause numerical issues)
def sigmoid_naive(z):
    return 1 / (1 + np.exp(-z))

# Naive binary cross-entropy loss
def bce_loss_naive(y_true, y_pred):
    n = len(y_true)
    loss = -(1/n) * np.sum(
        y_true * np.log(y_pred) + (1 - y_true) * np.log(1 - y_pred)
    )
    return loss

# Test with extreme values to see the problem
z_extreme = np.array([500, -500, 0])
print(f"Naive sigmoid on extreme values: {sigmoid_naive(z_extreme)}")
preds_broken = sigmoid_naive(z_extreme)
print(f"Naive BCE with sigmoid=1.0: {bce_loss_naive(np.array([1, 0, 1]), preds_broken)}")

Before you run this, do you think it will produce valid numbers or break? What happens when np.log(0) is computed?

  • Save and run:
python logistic_regression.py

You will see a RuntimeWarning about overflow in exp and the BCE loss returns nan or inf. This is the designed failure.

Why does this break?

When z = 500, np.exp(-500) underflows to 0.0, making sigmoid return exactly 1.0. Then np.log(1 - 1.0) = np.log(0) = negative infinity.

This is numerical instability. The math is correct on paper, but floating-point representation breaks at the extremes. Every production ML library handles this internally.

Fix it with numerically stable implementations
  • Replace the naive functions with stable versions:
# Numerically stable sigmoid using conditional logic
def sigmoid(z):
    # For large positive z: 1 / (1 + exp(-z)) is safe
    # For large negative z: exp(z) / (1 + exp(z)) avoids overflow
    pos_mask = z >= 0
    neg_mask = ~pos_mask
    result = np.zeros_like(z, dtype=float)
    result[pos_mask] = 1 / (1 + np.exp(-z[pos_mask]))
    ez = np.exp(z[neg_mask])
    result[neg_mask] = ez / (1 + ez)
    return result

# Numerically stable BCE using log-sum-exp trick
def bce_loss(y_true, z, w, lam):
    n = len(y_true)
    # log(sigmoid(z)) = z - log(1 + exp(z)) [stable form]
    # log(1-sigmoid(z)) = -log(1 + exp(z)) [stable form]
    log1pexp = np.logaddexp(0, z)  # log(1 + exp(z)) without overflow
    loss = (1/n) * np.sum(-y_true * z + log1pexp)
    # L2 regularization penalty
    reg = (lam / (2 * n)) * np.sum(w ** 2)
    return loss + reg

# Verify it handles extremes
print(f"Stable sigmoid on extreme values: {sigmoid(z_extreme)}")
print(f"Stable BCE: {bce_loss(np.array([1, 0, 1]), z_extreme, np.zeros(2), 0.0)}")

What makes this stable?

  • The boolean mask pos_mask = z >= 0 splits computation into two safe branches. Positive z uses the standard formula. Negative z uses the equivalent exp(z)/(1+exp(z)) form that never overflows.
  • np.logaddexp(0, z) computes log(1 + exp(z)) without ever computing exp(z) directly. NumPy handles the extreme cases internally.
  • The lam parameter is the L2 regularization strength. It penalizes large weights, preventing the model from overfitting to noise.
  • Save and run:
python logistic_regression.py

Now you should see valid sigmoid values [1.0, 0.0, 0.5] (approximately) and a finite BCE loss instead of NaN.

Implement gradient descent with L2 regularization
  • Add the gradient descent function:
# Gradient descent for logistic regression with L2 regularization
def train_logistic(X, y, lr, lam, iterations):
    n, d = X.shape
    w = np.zeros(d)
    b = 0.0
    losses = []

    for i in range(iterations):
        z = X @ w + b
        preds = sigmoid(z)
        errors = preds - y

        # Gradients with L2 regularization on weights
        dw = (1/n) * (X.T @ errors) + (lam/n) * w
        db = (1/n) * np.sum(errors)

        w -= lr * dw
        b -= lr * db

        if i % 100 == 0:
            loss = bce_loss(y, z, w, lam)
            losses.append(loss)
            if i % 500 == 0:
                print(f"Iter {i:4d} | Loss: {loss:.4f}")

    return w, b, losses

What does this code do?

  • The gradient dw has two parts: the data gradient X.T @ errors (how wrong are the predictions) plus the regularization term (lam/n) * w (penalize large weights).
  • The hyperparameter lr controls step size. Too large and you overshoot. Too small and training takes forever. You will tune this in the next step.
  • The lam hyperparameter controls regularization strength. Higher values mean simpler models that generalize better but may underfit.
  • Add a quick training test below the function to verify it works:
# Quick test: train on the PCA-reduced data
w_test, b_test, losses_test = train_logistic(X_pca, y, lr=0.1, lam=0.01, iterations=1000)
print(f"\nFinal weights: {w_test.round(4)}, bias: {b_test:.4f}")
print(f"Final loss: {losses_test[-1]:.4f}")
  • Save and run:
python logistic_regression.py

You should see the loss decreasing every 500 iterations, ending around 0.4-0.5. The model is learning to separate the classes.

Loss not decreasing?

  • Check that errors = preds - y subtracts in the right order. Reversing it flips the gradient and the model diverges.
  • Confirm you are passing X_pca (the 2-column PCA result), not X (the original 8-column data). Shape mismatch causes incorrect gradients.

Help me debug my logistic regression gradient descent.

Your model trains and converges. But how do you know if lr=0.1 and lam=0.01 are the right choices? Next, you will use k-fold cross-validation to find out.

Train with K-Fold Cross-Validation

Your model trains, but you picked hyperparameters by guessing. A model that looks great on training data might be memorizing noise instead of learning the pattern.

K-fold cross-validation splits data into k parts, trains on k-1, and validates on the remaining fold. Rotating which fold is held out gives you a robust estimate of real-world performance, and you will use it to tune learning rate and regularization strength.

In this step, get ready to:
  • Split data into train and test sets, holding out a final test set.
  • Implement k-fold cross-validation from scratch.
  • Tune hyperparameters and evaluate with precision, recall, and F1.
Split into train and test sets

The test set is sacred. You never train on it, never tune on it. It exists only to give an unbiased final score after all decisions are made.

  • Add the train/test split below your training test in logistic_regression.py:
# --- K-Fold Cross-Validation ---
# Shuffle the data first (permutation)
np.random.seed(42)
perm = np.random.permutation(len(y))
X_shuffled = X_pca[perm]
y_shuffled = y[perm]

# Hold out 20% as final test set
split_idx = int(0.8 * len(y))
X_train_full, X_test = X_shuffled[:split_idx], X_shuffled[split_idx:]
y_train_full, y_test = y_shuffled[:split_idx], y_shuffled[split_idx:]

print(f"Train set: {X_train_full.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")
print(f"Test class distribution: {np.sum(y_test==0)} neg, {np.sum(y_test==1)} pos")

Why shuffle with permutation?

Your data was created with all class 0 samples first, then all class 1. Without shuffling, the test set would be entirely class 1. np.random.permutation creates a random index order that mixes both classes throughout the splits.

  • Save and run:
python logistic_regression.py

You should see 240 training samples and 60 test samples with both classes represented in the test set.

Implement k-fold cross-validation
  • Add the k-fold implementation:
# K-fold cross-validation from scratch
def k_fold_cv(X, y, k, lr, lam, iterations):
    n = len(y)
    fold_size = n // k
    fold_scores = []

    for fold in range(k):
        # Create boolean mask for validation fold
        val_mask = np.zeros(n, dtype=bool)
        val_mask[fold * fold_size:(fold + 1) * fold_size] = True
        train_mask = ~val_mask

        X_train, X_val = X[train_mask], X[val_mask]
        y_train, y_val = y[train_mask], y[val_mask]

        # Train on this fold
        w, b, _ = train_logistic(X_train, y_train, lr, lam, iterations)

        # Predict on validation fold
        z_val = X_val @ w + b
        preds_val = (sigmoid(z_val) >= 0.5).astype(int)

        # Compute precision, recall, F1
        tp = np.sum((preds_val == 1) & (y_val == 1))
        fp = np.sum((preds_val == 1) & (y_val == 0))
        fn = np.sum((preds_val == 0) & (y_val == 1))

        precision = tp / (tp + fp) if (tp + fp) > 0 else 0.0
        recall = tp / (tp + fn) if (tp + fn) > 0 else 0.0
        f1 = 2 * precision * recall / (precision + recall) if (precision + recall) > 0 else 0.0
        fold_scores.append(f1)

    return np.mean(fold_scores), fold_scores

What does this code do?

  • The boolean mask val_mask selects one fold for validation. ~val_mask (bitwise NOT) gives you the training set. No data leaks between folds.
  • Precision, recall, and F1 are computed using boolean operations: (preds_val == 1) & (y_val == 1) counts true positives without any loops.
  • The function returns the mean F1 across all k folds. Higher means the model generalizes better at those hyperparameters.
Tune hyperparameters with cross-validation
  • Add the hyperparameter search:
# Hyperparameter search
learning_rates = [0.01, 0.1, 0.5]
reg_strengths = [0.001, 0.01, 0.1, 1.0]

print("\n--- Hyperparameter Search (5-fold CV) ---")
best_f1, best_lr, best_lam = 0, 0, 0

for lr in learning_rates:
    for lam in reg_strengths:
        mean_f1, _ = k_fold_cv(X_train_full, y_train_full, k=5, lr=lr, lam=lam, iterations=1000)
        if mean_f1 > best_f1:
            best_f1 = mean_f1
            best_lr = lr
            best_lam = lam

print(f"\nBest: lr={best_lr}, lam={best_lam}, CV F1={best_f1:.3f}")
  • Save and run:
python logistic_regression.py

The search tests 12 combinations. You should see the best F1 score printed, likely with lr=0.1 or lr=0.5 and a moderate regularization value.

Final evaluation on the held-out test set
  • Add the final evaluation using the best hyperparameters:
# Final model: train on full training set, evaluate on test set
w_final, b_final, _ = train_logistic(
    X_train_full, y_train_full, lr=best_lr, lam=best_lam, iterations=2000
)

# Predict on test set
z_test = X_test @ w_final + b_final
preds_test = (sigmoid(z_test) >= 0.5).astype(int)

# Final metrics
tp = np.sum((preds_test == 1) & (y_test == 1))
fp = np.sum((preds_test == 1) & (y_test == 0))
fn = np.sum((preds_test == 0) & (y_test == 1))
tn = np.sum((preds_test == 0) & (y_test == 0))

accuracy = (tp + tn) / len(y_test)
precision = tp / (tp + fp) if (tp + fp) > 0 else 0.0
recall = tp / (tp + fn) if (tp + fn) > 0 else 0.0
f1 = 2 * precision * recall / (precision + recall) if (precision + recall) > 0 else 0.0

print(f"\n--- Final Test Results ---")
print(f"Accuracy:  {accuracy:.3f}")
print(f"Precision: {precision:.3f}")
print(f"Recall:    {recall:.3f}")
print(f"F1 Score:  {f1:.3f}")

Why these metrics matter together

Accuracy tells you overall correctness but hides class imbalance problems. Precision tells you how trustworthy your positive predictions are. Recall tells you how many real positives you caught.

A spam filter with 99% precision but 10% recall is useless: it only catches 1 in 10 spam emails. F1 balances both into a single number you can optimize.

  • Save and run the final script:
python logistic_regression.py

You should see accuracy around 80-85%, with precision, recall, and F1 all reported. The F1 should be close to your cross-validation estimate, confirming the model generalizes.

F1 score is 0.0?

  • If F1 is zero, the model is predicting all zeros (never predicting class 1). Check that best_lr is not too small. A very low learning rate with only 1000 iterations may not converge.
  • Confirm your data was shuffled before splitting. If np.random.permutation is missing, all class 1 samples may be in the test set with none in training.

Help me debug my k-fold cross-validation.

Secret mission

Explore the Curse of Dimensionality

Your PCA step crushed 8 features down to 2. But was that really necessary? What if you skipped PCA entirely and fed all 8 features straight into logistic regression? In this mission, you'll find out exactly what goes wrong when dimensions grow unchecked.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs. Decide whether to keep your files for future reference, or delete them to tidy up your Desktop.

Resources you used:

  • linear_regression.py on your Desktop.
  • linear_regression_multi.py on your Desktop (Secret Mission).

Keep everything running

No action needed. Choose this if you plan to keep experimenting with your gradient descent implementation or want the code as a reference for future ML projects.

  • Your Python files remain on your Desktop, ready to run anytime.
  • NumPy and Matplotlib stay installed in your Python environment for other projects.

Pause - I'll come back to this later

Since everything is local files, there is nothing to pause. Your files sit on disk and consume negligible space.

  • Leave your files where they are. They take up only a few kilobytes.
  • When you return, run them again with the same commands you used during the project.

Delete - I don't want to use this again

Remove all project files from your Desktop to start fresh.

  • Delete the project files by running this command in your terminal:

macOS

rm ~/Desktop/linear_regression.py ~/Desktop/linear_regression_multi.py

Windows

Remove-Item ~\Desktop\linear_regression.py, ~\Desktop\linear_regression_multi.py
  • Confirm the files are gone by running:

macOS

ls ~/Desktop/linear_regression*.py

Windows

Get-ChildItem ~\Desktop\linear_regression*.py

You should see a "No such file" error (macOS) or no output (Windows), confirming both files are deleted.

Nice Work!

Nice Work!

You just built a complete logistic regression classifier from pure math and NumPy. No scikit-learn, no magic one-liners. Every sigmoid, every gradient update, every cross-validation fold came from code you wrote.

You've learned how to:

  • Implement PCA from scratch using eigendecomposition of a covariance matrix to reduce 8 noisy features down to 2 principal components.
  • Build a numerically stable logistic regression with the log-sum-exp trick and L2 regularization, after deliberately triggering NaN to understand why naive math fails.
  • Evaluate with k-fold cross-validation that properly isolates training from validation data, tunes hyperparameters, and reports precision, recall, and F1.
  • Secret Mission: Explored the curse of dimensionality by comparing full-dimension vs PCA-reduced model performance, observing overfitting emerge as feature count grows.

Ready to quiz yourself?