Lecture 08: Bias-Variance and Logistic Regression (1) – CS 189, Fall 2026

Reference: Bishop & Bishop, Deep Learning: Foundations and Concepts, §4.3 (Part 0), Ch. 5 §5.1 and §5.4.1-5.4.3 (Parts 1-6).

In [1]:
import numpy as np
import pandas as pd

# Part 0 (bias-variance)
import matplotlib.pyplot as plt
from numpy.polynomial import legendre as L

# Parts 1-6 (logistic regression)
import plotly.express as px
import plotly.graph_objects as go
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.preprocessing import StandardScaler
from scipy.optimize import minimize

import warnings
from sklearn.exceptions import ConvergenceWarning
warnings.filterwarnings("ignore", category=ConvergenceWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", message="Setting penalty=None")

import plotly.io as pio
pio.renderers.default = "notebook_connected"

C_DATA, C_FIT, C_ALT, C_RESID, C_SPAN = "#003262", "#FDB515", "#C4820E", "#D55E00", "#00553A"
plt.rcParams.update({"figure.figsize": (8, 5), "axes.grid": True, "grid.alpha": 0.3,
                     "font.size": 13, "axes.titlesize": 15, "axes.labelsize": 14})

Part 0. The Bias-Variance Decomposition (carried over from Lecture 07)¶

Take the expected squared error at an input $x$, averaged over training sets $\mathcal{D}$:

$$ \mathbb{E}_{\mathcal{D}}\big[(t - y(x;\mathcal{D}))^2\big] = \underbrace{\sigma^2}_{\text{noise}} + \underbrace{\big(\mathbb{E}_{\mathcal{D}}[y(x;\mathcal{D})] - f(x)\big)^2}_{\text{bias}^2} + \underbrace{\mathbb{E}_{\mathcal{D}}\big[(y(x;\mathcal{D}) - \mathbb{E}_{\mathcal{D}}[y(x;\mathcal{D})])^2\big]}_{\text{variance}} $$
  • Noise is irreducible. No model removes it.
  • Bias is how far the average fitted function sits from the truth: underfitting.
  • Variance is how much the fitted function moves when the training set changes: overfitting.

Every term is an expectation over training sets, so we can estimate all three by simulation: draw many training sets, fit each, and look at the spread of the fitted functions.

In [2]:
SIGMA_BV, NSETS = 0.30, 400
f_true = lambda x: np.sin(5 * x)
xg_bv = np.linspace(-0.9, 0.9, 120)

def Phi_leg(x, deg):
    "Legendre basis: well conditioned, so what we measure is variance and not rounding error."
    return L.legvander(np.asarray(x).ravel(), deg)

def decompose(n_train, degrees, seed=189):
    rng = np.random.default_rng(seed)
    out = []
    for d in degrees:
        preds = np.empty((NSETS, len(xg_bv)))
        for s in range(NSETS):
            xt = rng.uniform(-0.9, 0.9, n_train)
            tt = f_true(xt) + SIGMA_BV * rng.standard_normal(n_train)
            w = np.linalg.lstsq(Phi_leg(xt, d), tt, rcond=None)[0]
            preds[s] = Phi_leg(xg_bv, d) @ w
        bias2 = float(np.mean((preds.mean(0) - f_true(xg_bv)) ** 2))
        var = float(np.mean(preds.var(0)))
        out.append((d, bias2, var, SIGMA_BV ** 2, bias2 + var + SIGMA_BV ** 2, preds))
    return out

degrees = list(range(0, 11))
res100 = decompose(100, degrees)

print(f"N_train = 100,  {NSETS} training sets per degree\n")
print(f"{'deg':>4} {'bias^2':>9} {'variance':>9} {'noise':>8} {'total':>9}")
for d, b2, v, nz, tot, _ in res100:
    print(f"{d:>4} {b2:>9.4f} {v:>9.4f} {nz:>8.4f} {tot:>9.4f}")
best = min(res100, key=lambda r: r[4])
print(f"\nlowest expected error at degree {best[0]}  (total {best[4]:.4f})")

N_train = 100,  400 training sets per degree

 deg    bias^2  variance    noise     total
   0    0.4811    0.0053   0.0900    0.5764
   1    0.4813    0.0112   0.0900    0.5824
   2    0.4813    0.0191   0.0900    0.5904
   3    0.0726    0.0080   0.0900    0.1706
   4    0.0728    0.0120   0.0900    0.1748
   5    0.0017    0.0067   0.0900    0.0984
   6    0.0018    0.0084   0.0900    0.1001
   7    0.0000    0.0097   0.0900    0.0998
   8    0.0000    0.0128   0.0900    0.1028
   9    0.0001    0.0204   0.0900    0.1105
  10    0.0000    0.0329   0.0900    0.1229

lowest expected error at degree 5  (total 0.0984)

0.1 Seeing bias and variance directly¶

Plot the fitted function from 40 of the training sets, plus their average, against the truth. Bias is the gap between the average (dashed) and the truth (solid). Variance is the width of the grey band.

In [3]:
show_degrees = [1, 3, 5, 9]
fig, axs = plt.subplots(1, 4, figsize=(17, 4.2), sharey=True)
for ax, d in zip(axs, show_degrees):
    rec = [r for r in res100 if r[0] == d][0]
    preds = rec[5]
    for p in preds[:40]:
        ax.plot(xg_bv, p, color="0.55", lw=0.8, alpha=0.55)
    ax.plot(xg_bv, f_true(xg_bv), lw=3, color=C_DATA, label="truth $f(x)$")
    ax.plot(xg_bv, preds.mean(0), lw=3, ls="--", color=C_ALT, label="average fit")
    ax.set_title(f"degree {d}\nbias$^2$ = {rec[1]:.3f}\nvariance = {rec[2]:.3f}", fontsize=13)
    ax.set(xlabel="$x$", ylim=(-2.0, 2.0))
axs[0].set_ylabel("$y$"); axs[0].legend(fontsize=10, loc="lower left")
plt.tight_layout(); plt.show()
No description has been provided for this image
In [4]:
b2 = np.array([r[1] for r in res100])
vr = np.array([r[2] for r in res100])
nz = np.array([r[3] for r in res100])
tot = b2 + vr + nz

fig, ax = plt.subplots(figsize=(8.5, 5))
ax.plot(degrees, b2, "o-", lw=3, ms=8, color=C_DATA, label="bias$^2$")
ax.plot(degrees, vr, "s-", lw=3, ms=8, color=C_RESID, label="variance")
ax.plot(degrees, nz, ":", lw=2.5, color="0.4", label="noise $\\sigma^2$")
ax.plot(degrees, tot, "^-", lw=3, ms=8, color=C_SPAN, label="total")
ax.axvline(degrees[int(tot.argmin())], color=C_ALT, ls="--", lw=2.5,
           label=f"minimum at degree {degrees[int(tot.argmin())]}")
ax.set(xlabel="polynomial degree", ylabel="expected squared error", yscale="log",
       xticks=degrees)
ax.legend(fontsize=11); plt.tight_layout(); plt.show()
No description has been provided for this image

0.2 More data supports more complexity¶

Nothing about the tradeoff is a property of the model alone. Shrink the training set and the variance term grows at every degree, so the minimum moves left: with less data the best model is simpler than the truth.

In [5]:
res25 = decompose(25, degrees)
tot25 = np.array([r[4] for r in res25])

fig, ax = plt.subplots(figsize=(8.5, 5))
for res, tot_i, col, name in [(res25, tot25, C_RESID, "N = 25"),
                              (res100, tot, C_DATA, "N = 100")]:
    ax.plot(degrees, tot_i, "o-", lw=3, ms=8, color=col, label=f"{name}, total error")
    k = int(np.argmin(tot_i))
    ax.plot(degrees[k], tot_i[k], "*", ms=22, color=C_ALT, mec="k", zorder=5)
ax.set(xlabel="polynomial degree", ylabel="expected squared error", yscale="log",
       xticks=degrees)
ax.legend(fontsize=12); plt.tight_layout(); plt.show()

print(f"N = 25  -> minimum at degree {degrees[int(tot25.argmin())]} (total {tot25.min():.4f})")
print(f"N = 100 -> minimum at degree {degrees[int(tot.argmin())]} (total {tot.min():.4f})")
No description has been provided for this image 
N = 25  -> minimum at degree 3 (total 0.2151)
N = 100 -> minimum at degree 5 (total 0.0984)

0.3 Regularization moves along the same curve¶

Model complexity is not only the number of basis functions. Hold the degree fixed at 9 and turn $\lambda$ instead: large $\lambda$ buys bias and sells variance, exactly as reducing the degree does.

In [6]:
def decompose_ridge(n_train, deg, lambdas, seed=189):
    rng = np.random.default_rng(seed)
    out = []
    for lam in lambdas:
        preds = np.empty((NSETS, len(xg_bv)))
        for s in range(NSETS):
            xt = rng.uniform(-0.9, 0.9, n_train)
            tt = f_true(xt) + SIGMA_BV * rng.standard_normal(n_train)
            A = Phi_leg(xt, deg)
            P = np.eye(deg + 1); P[0, 0] = 0.0           # intercept left unpenalized
            w = np.linalg.solve(A.T @ A + lam * P, A.T @ tt)
            preds[s] = Phi_leg(xg_bv, deg) @ w
        bias2 = float(np.mean((preds.mean(0) - f_true(xg_bv)) ** 2))
        var = float(np.mean(preds.var(0)))
        out.append((lam, bias2, var, bias2 + var + SIGMA_BV ** 2))
    return out

lams = np.logspace(-4, 3, 22)
rr = decompose_ridge(100, 9, lams)
b2r = np.array([r[1] for r in rr]); vrr = np.array([r[2] for r in rr])
totr = np.array([r[3] for r in rr])

fig, ax = plt.subplots(figsize=(8.5, 5))
ax.plot(lams, b2r, "o-", lw=3, ms=7, color=C_DATA, label="bias$^2$")
ax.plot(lams, vrr, "s-", lw=3, ms=7, color=C_RESID, label="variance")
ax.plot(lams, totr, "^-", lw=3, ms=7, color=C_SPAN, label="total")
ax.axvline(lams[int(totr.argmin())], color=C_ALT, ls="--", lw=2.5,
           label=rf"best $\lambda$ = {lams[int(totr.argmin())]:.3g}")
ax.set(xscale="log", yscale="log", xlabel=r"$\lambda$  (degree fixed at 9)",
       ylabel="expected squared error")
ax.legend(fontsize=11); plt.tight_layout(); plt.show()

print(f"lambda -> 0     : bias^2 {b2r[0]:.4f}, variance {vrr[0]:.4f}")
print(f"best lambda     : bias^2 {b2r[totr.argmin()]:.4f}, variance {vrr[totr.argmin()]:.4f}"
      f"  (lambda = {lams[int(totr.argmin())]:.3g})")
print(f"lambda -> large : bias^2 {b2r[-1]:.4f}, variance {vrr[-1]:.4f}")
No description has been provided for this image 
lambda -> 0     : bias^2 0.0000, variance 0.0143
best lambda     : bias^2 0.0003, variance 0.0078  (lambda = 0.215)
lambda -> large : bias^2 0.4699, variance 0.0052

Everything above assumed squared error and a real-valued target. Today's second half changes the target to a label in $\{0, 1\}$. The first thing to check is whether squared error still makes sense, and the answer will be no: not for the model, and not for the error function.

The Data: Wisconsin Breast Cancer¶

569 tumor biopsies, 30 measurements each. In scikit-learn, target = 0 means malignant, so we flip it into a malignant column where 1 = malignant, 0 = benign.

In [7]:
data_dict = load_breast_cancer()
data = pd.DataFrame(data_dict["data"], columns=data_dict["feature_names"])
data["malignant"] = (data_dict["target"] == 0).astype(int)

print("Fraction malignant:", data["malignant"].mean().round(3))
data[["mean radius", "mean texture", "malignant"]].head()

Fraction malignant: 0.373
Out[7]:
mean radius mean texture malignant
0 17.99 10.38 1
1 20.57 17.77 1
2 19.69 21.25 1
3 11.42 20.38 1
4 20.29 14.34 1

One split for the whole unit¶

We make one train / validation / test split (60 / 20 / 20) and use it in both logistic regression lectures. Today we only touch train and test. The validation set is for choosing thresholds next lecture.

In [8]:
data_tr, data_te = train_test_split(data, test_size=0.2, random_state=42, stratify=data["malignant"])
data_tr, data_va = train_test_split(data_tr, test_size=0.25, random_state=42, stratify=data_tr["malignant"])

X_train, y_train = data_tr[["mean radius"]].to_numpy(), data_tr["malignant"].to_numpy()
X_test,  y_test  = data_te[["mean radius"]].to_numpy(), data_te["malignant"].to_numpy()
print(len(data_tr), "train /", len(data_va), "validation /", len(data_te), "test")

341 train / 114 validation / 114 test

Plotting helpers¶

Every plot below is "points colored by class plus some curves". We write that once. Points are jittered vertically so the 0s and 1s do not sit on top of each other.

In [9]:
rng = np.random.default_rng(0)

def jitter(y, amt=0.1):
    return y + amt * rng.random(len(y)) - amt / 2

def class_scatter(X, y, title=""):
    fig = go.Figure()
    for label, name, color in [(0, "Benign", "blue"), (1, "Malignant", "red")]:
        m = (y == label)
        fig.add_scatter(x=X[m].ravel(), y=jitter(y[m]), mode="markers", name=name,
                        marker=dict(color=color, opacity=0.5))
    fig.update_layout(title=title, template="plotly_white",
                      xaxis_title="Mean Radius", yaxis_title="Malignant")
    return fig

x_grid = np.linspace(5, 30, 200).reshape(-1, 1)
class_scatter(X_train, y_train, "Mean Radius vs Malignant (jittered)")

Question: Looking at this plot, what rule would you use to predict malignancy?

Step 1: Try Least Squares Anyway¶

"I suppose it is tempting, if the only tool you have is a hammer, to treat everything as if it were a nail." - Abraham Maslow

The labels are 0 and 1. Those are numbers, so nothing stops us from fitting ordinary least squares. No regularization here: one feature, 341 points.

In [10]:
lin_reg = LinearRegression().fit(X_train, y_train)

fig = class_scatter(X_train, y_train, "Least Squares Fit to a 0/1 Target")
fig.add_scatter(x=x_grid.ravel(), y=lin_reg.predict(x_grid), mode="lines",
                name="Least squares", line=dict(color="orange", width=3))
fig

To classify, we need a decision rule: predict malignant if $f(x) > 0.5$. For a line $f(x) = w_0 + w_1 x$ the boundary is where $f(x) = 0.5$, i.e. $x^* = (0.5 - w_0)/w_1$.

In [11]:
def ls_boundary(model):
    return (0.5 - model.intercept_) / model.coef_[0]

ls_x_star = ls_boundary(lin_reg)
ls_train_acc = np.mean((lin_reg.predict(X_train) > 0.5) == y_train)
print(f"Least-squares boundary: x* = {ls_x_star:.2f}")
print(f"Training accuracy:      {ls_train_acc:.3f}")

Least-squares boundary: x* = 15.45
Training accuracy:      0.886

Accuracy looks fine. But what does a prediction of 1.3, or -0.2, mean? Clipping to [0, 1] hides the problem; it does not make the output a probability.

In [12]:
fig = class_scatter(X_train, y_train, "Clipping the Line to [0, 1]")
fig.add_scatter(x=x_grid.ravel(), y=lin_reg.predict(x_grid), mode="lines",
                name="Least squares", line=dict(color="orange", dash="dot"))
fig.add_scatter(x=x_grid.ravel(), y=np.clip(lin_reg.predict(x_grid), 0, 1), mode="lines",
                name="Clipped", line=dict(color="green", width=5))
fig

The failure: one very large tumor¶

Suppose we observe a malignant tumor with mean radius 100. It is obviously malignant, and it is on the correct side of the boundary already. Watch what least squares does with it.

In [13]:
X_ex = np.vstack([X_train, [[100.0]]])
y_ex = np.append(y_train, 1)

lin_reg_ex = LinearRegression().fit(X_ex, y_ex)
ls_x_star_ex = ls_boundary(lin_reg_ex)
ls_acc_ex = np.mean((lin_reg_ex.predict(X_train) > 0.5) == y_train)

print(f"Boundary moved from {ls_x_star:.2f} to {ls_x_star_ex:.2f}")
print(f"Training accuracy (original points) fell from {ls_train_acc:.3f} to {ls_acc_ex:.3f}")

Boundary moved from 15.45 to 17.43
Training accuracy (original points) fell from 0.886 to 0.809
In [14]:
x_wide = np.linspace(5, 105, 300).reshape(-1, 1)
fig = class_scatter(X_ex, y_ex, "One Easy Point Moves the Boundary")
fig.add_scatter(x=x_wide.ravel(), y=lin_reg.predict(x_wide), mode="lines",
                name="Before", line=dict(color="gray"))
fig.add_scatter(x=x_wide.ravel(), y=lin_reg_ex.predict(x_wide), mode="lines",
                name="After adding x = 100", line=dict(color="orange"))
fig.add_vline(x=ls_x_star, line_dash="dot", line_color="gray")
fig.add_vline(x=ls_x_star_ex, line_dash="dash", line_color="black")
fig.update_yaxes(range=[-0.5, 1.5])
fig

Least squares penalizes the new point for being "too far above the line", so it tilts the line toward it. That is the Gaussian-noise assumption behind least squares at work: it treats a correct, confident point as a large residual. We need a different model and a different objective.

Step 2: The Logistic Sigmoid¶

$$\sigma(a) = \frac{1}{1 + e^{-a}}$$

We derive it in lecture from the assumption that the log-odds are linear in $x$.

In [15]:
def sigmoid(a):
    return 1 / (1 + np.exp(-a))

a = np.linspace(-8, 8, 200)
px.line(x=a, y=sigmoid(a), labels={"x": "a", "y": "σ(a)"},
        title="The Logistic Sigmoid", template="plotly_white")

Checking the properties numerically¶

Each property from the slide, verified on a grid of values instead of taken on faith.

In [16]:
a = np.linspace(-6, 6, 13)
h = 1e-5

print("Symmetry  1 - σ(a) = σ(-a):         ", np.allclose(1 - sigmoid(a), sigmoid(-a)))
numeric_deriv = (sigmoid(a + h) - sigmoid(a - h)) / (2 * h)
print("Derivative σ'(a) = σ(a)(1 - σ(a)):  ", np.allclose(numeric_deriv, sigmoid(a) * (1 - sigmoid(a))))
p = sigmoid(a)
print("Inverse   logit(σ(a)) = a:          ", np.allclose(np.log(p / (1 - p)), a))

Symmetry  1 - σ(a) = σ(-a):          True
Derivative σ'(a) = σ(a)(1 - σ(a)):   True
Inverse   logit(σ(a)) = a:           True

Shifting and scaling: $\sigma(Ax + B)$¶

Bigger $A$ makes the curve steeper. Changing $B$ shifts it; the midpoint sits at $x = -B/A$.

In [17]:
x = np.linspace(-6, 6, 200)
fig = go.Figure()
for A, B in [(1, 0), (3, 0), (1, 2)]:
    fig.add_scatter(x=x, y=sigmoid(A * x + B), mode="lines", name=f"A={A}, B={B}")
fig.update_layout(title="σ(Ax + B)", template="plotly_white")
fig

Step 3: Fit Logistic Regression¶

Heads up: scikit-learn's LogisticRegression adds L2 regularization by default (C=1.0). We have not introduced regularization yet, so we turn it off with C=np.inf. We come back to C next lecture.

In [18]:
log_reg = LogisticRegression(C=np.inf).fit(X_train, y_train)

w0, w1 = log_reg.intercept_[0], log_reg.coef_[0, 0]
print(f"w0 (intercept) = {w0:.3f}")
print(f"w1 (slope)     = {w1:.3f}")

w0 (intercept) = -14.809
w1 (slope)     = 0.997

predict_proba returns one column per class, in the order of log_reg.classes_. Column 1 is $P(\text{malignant} \mid x)$.

In [19]:
print(log_reg.classes_)
log_reg.predict_proba(X_train[:5]).round(3)

[0 1]
Out[19]:
array([[0.499, 0.501],
       [0.999, 0.001],
       [0.974, 0.026],
       [0.789, 0.211],
       [0.815, 0.185]])

The real decision boundary¶

$\sigma(w_0 + w_1 x) = 0.5$ exactly when $w_0 + w_1 x = 0$, so $x^* = -w_0 / w_1$. (The original notebook shaded the regions at the mean radius. That is not the boundary.)

In [20]:
lr_x_star = -w0 / w1
print(f"Logistic regression boundary: x* = {lr_x_star:.2f}")

fig = class_scatter(X_train, y_train, "Logistic Regression vs Least Squares")
fig.add_scatter(x=x_grid.ravel(), y=log_reg.predict_proba(x_grid)[:, 1], mode="lines",
                name="Logistic regression  P(y=1|x)", line=dict(color="black", width=3))
fig.add_scatter(x=x_grid.ravel(), y=lin_reg.predict(x_grid), mode="lines",
                name="Least squares", line=dict(color="orange", dash="dot"))
fig.add_vrect(x0=5, x1=lr_x_star, fillcolor="blue", opacity=0.07, line_width=0)
fig.add_vrect(x0=lr_x_star, x1=30, fillcolor="red", opacity=0.07, line_width=0)
fig.add_vline(x=lr_x_star, line_dash="dash", annotation_text="decision boundary")
fig

Logistic regression boundary: x* = 14.86

Does the large tumor still break things?¶

In [21]:
log_reg_ex = LogisticRegression(C=np.inf).fit(X_ex, y_ex)
lr_x_star_ex = -log_reg_ex.intercept_[0] / log_reg_ex.coef_[0, 0]

print(f"Least squares boundary:        {ls_x_star:.2f} -> {ls_x_star_ex:.2f}")
print(f"Logistic regression boundary:  {lr_x_star:.2f} -> {lr_x_star_ex:.2f}")

Least squares boundary:        15.45 -> 17.43
Logistic regression boundary:  14.86 -> 14.86

The point at radius 100 already has $\sigma(\cdot) \approx 1$, so it contributes almost nothing to the cross-entropy gradient. Logistic regression ignores points it already gets right with confidence.

In [22]:
print(f"Test accuracy, least squares:       {np.mean((lin_reg.predict(X_test) > 0.5) == y_test):.3f}")
print(f"Test accuracy, logistic regression: {log_reg.score(X_test, y_test):.3f}")

Test accuracy, least squares:       0.860
Test accuracy, logistic regression: 0.860

Same test accuracy on this split. The win is not accuracy on clean data: it is outputs we can read as probabilities and a boundary that does not move when we see an easy point.

Reading the weights: odds¶

The model says $\log \frac{p}{1-p} = w_0 + w_1 x$. Increasing $x$ by 1 adds $w_1$ to the log-odds, which multiplies the odds by $e^{w_1}$.

In [23]:
print(f"Each extra unit of mean radius multiplies the odds of malignancy by e^w1 = {np.exp(w1):.2f}")

for r in [12, 13]:
    p = log_reg.predict_proba([[r]])[0, 1]
    print(f"radius {r}: p = {p:.3f}, odds = {p / (1 - p):.3f}")

Each extra unit of mean radius multiplies the odds of malignancy by e^w1 = 2.71
radius 12: p = 0.055, odds = 0.058
radius 13: p = 0.136, odds = 0.157

Step 4: Two Features, Linear Boundary¶

With $D$ features the boundary $w^\top x + w_0 = 0$ is a $(D-1)$-dimensional hyperplane, and $w$ is orthogonal to it. Here $D = 2$: the boundary is a line. We standardize the features so the two axes are on the same scale and the orthogonality is visible.

In [24]:
feats2 = ["mean radius", "mean texture"]
scaler2 = StandardScaler().fit(data_tr[feats2])
X2 = scaler2.transform(data_tr[feats2])

log_reg2 = LogisticRegression(C=np.inf).fit(X2, y_train)
w, b = log_reg2.coef_[0], log_reg2.intercept_[0]
print("w =", w.round(3), " w0 =", round(b, 3))

w = [3.657 1.01 ]  w0 = -0.727
In [25]:
x1 = np.linspace(-2.5, 4, 50)
x2_boundary = -(b + w[0] * x1) / w[1]

fig = px.scatter(x=X2[:, 0], y=X2[:, 1], color=y_train.astype(str),
                 color_discrete_map={"0": "blue", "1": "red"}, opacity=0.5,
                 labels={"x": "mean radius (std.)", "y": "mean texture (std.)", "color": "malignant"},
                 template="plotly_white", title="2-D Decision Boundary and the Weight Vector")
fig.add_scatter(x=x1, y=x2_boundary, mode="lines", name="w·x + w0 = 0", line=dict(color="black"))
origin = np.array([x1[20], x2_boundary[20]])
fig.add_annotation(x=origin[0] + w[0] / 3, y=origin[1] + w[1] / 3, ax=origin[0], ay=origin[1],
                   xref="x", yref="y", axref="x", ayref="y", showarrow=True, arrowhead=3,
                   arrowwidth=3, text="w")
fig.update_yaxes(range=[-2.5, 4], scaleanchor="x", scaleratio=1)
fig

/opt/anaconda3/lib/python3.11/site-packages/plotly/express/_core.py:1979: FutureWarning:

When grouping with a length-1 list-like, you will need to pass a length-1 tuple to get_group in a future version of pandas. Pass `(name,)` instead of `name` to silence this warning.

Step 5: Which Objective? Squared Error vs Cross-Entropy¶

A small toy dataset, no intercept, so the model is $p = \sigma(w x)$ with a single scalar $w$.

In [26]:
toy_df = pd.DataFrame({"x": [-4, -2, -0.5, 1, 3, 5],
                       "y": [0, 0, 1, 0, 1, 1]})
toy_df
Out[26]:
x y
0 -4.0 0
1 -2.0 0
2 -0.5 1
3 1.0 0
4 3.0 1
5 5.0 1

Squared error on top of the sigmoid is non-convex¶

In [27]:
def ss_error_on_toy_data(w):
    p_hat = sigmoid(toy_df["x"] * w)
    return 0.5 * np.sum((toy_df["y"] - p_hat) ** 2)

w_grid = np.linspace(-10, 10, 400)
w_error = pd.DataFrame({"w": w_grid, "Squared Error": [ss_error_on_toy_data(w) for w in w_grid]})
px.line(w_error, x="w", y="Squared Error", title="Squared Error on Toy Data", template="plotly_white")

Run the optimizer from two starting points. Same data, same objective, different answers.

In [28]:
best_w = minimize(ss_error_on_toy_data, x0=0)["x"][0]
stuck_w = minimize(ss_error_on_toy_data, x0=-5)["x"][0]
print(f"start at w =  0  ->  w = {best_w:.3f}")
print(f"start at w = -5  ->  w = {stuck_w:.3f}   (stuck on the flat plateau)")

start at w =  0  ->  w = 0.545
start at w = -5  ->  w = -10.858   (stuck on the flat plateau)
In [29]:
fig = px.scatter(toy_df, x="x", y="y", template="plotly_white")
fig.update_traces(marker_size=14)
xs = np.linspace(-6, 6, 200)
fig.add_scatter(x=xs, y=sigmoid(xs * best_w), mode="lines", name=f"w = {best_w:.2f}")
fig.add_scatter(x=xs, y=sigmoid(xs * stuck_w), mode="lines", name=f"w = {stuck_w:.2f}")
fig

Squared error is also bounded¶

For one point with $y = 1$: the worst squared error (with the ½) is 0.5, no matter how confidently wrong the model is. Negative log-probability goes to infinity.

In [30]:
p = np.linspace(0.001, 0.999, 200)
loss_y1 = pd.DataFrame({"p": p, "Squared Error": 0.5 * (1 - p) ** 2, "Neg Log": -np.log(p)})
px.line(loss_y1, x="p", y=["Squared Error", "Neg Log"], template="plotly_white",
        title="Loss for One Point with y = 1", labels={"value": "loss", "p": "predicted P(y=1)"})
In [31]:
loss_y0 = pd.DataFrame({"p": p, "Squared Error": 0.5 * p ** 2, "Neg Log": -np.log(1 - p)})
px.line(loss_y0, x="p", y=["Squared Error", "Neg Log"], template="plotly_white",
        title="Loss for One Point with y = 0", labels={"value": "loss", "p": "predicted P(y=1)"})

Cross-entropy on the same toy data¶

Computing $\log \sigma(a)$ directly returns -inf once $\sigma(a)$ rounds to 0 or 1. We use the identity $-\log \sigma(a) = \log(1 + e^{-a})$, which np.logaddexp(0, -a) computes stably.

In [32]:
def cross_entropy_from_logits(t, a):
    # t*(-log σ(a)) + (1-t)*(-log(1-σ(a))), written stably
    return t * np.logaddexp(0, -a) + (1 - t) * np.logaddexp(0, a)

def mean_ce_on_toy_data(w):
    return np.mean(cross_entropy_from_logits(toy_df["y"], toy_df["x"] * w))

w_error["Cross-Entropy"] = [mean_ce_on_toy_data(w) for w in w_grid]
px.line(w_error, x="w", y="Cross-Entropy", title="Cross-Entropy on Toy Data (convex)",
        template="plotly_white")
In [33]:
ce_w_from_0 = minimize(mean_ce_on_toy_data, x0=0)["x"][0]
ce_w_from_neg5 = minimize(mean_ce_on_toy_data, x0=-5)["x"][0]
print(f"start at  0 -> w = {ce_w_from_0:.3f}")
print(f"start at -5 -> w = {ce_w_from_neg5:.3f}   (same answer)")

start at  0 -> w = 0.743
start at -5 -> w = 0.743   (same answer)