Data: 17,898 human-labeled pulsar candidates from a radio-telescope survey

• Pulsars (\(Y = 1\)): 1,639, 9.2%
• Noise (\(Y = 0\)): 16,259

Measurement \(x\): the mean of the candidate’s pulse profile

• Pulsars average 56.7, noise 116.6

Decision: pulsar or noise, from \(x\) alone

Overlap: the same \(x\) occurs in both classes

• \(x\) from 70 to 80: 145 pulsars and 250 noise
• Any rule based on \(x\) alone misclassifies some candidates

Data: HTRU2, R. J. Lyon et al., MNRAS 2016; UCI Machine Learning Repository (CC BY 4.0)

Rule: decide pulsar if \(x < T\), noise otherwise

Two errors:

• False alarm: noise decided as pulsar, \(P(Y{=}0,\ X < T)\)
• Miss: pulsar decided as noise, \(P(Y{=}1,\ X \ge T)\)

Example, \(T = 80\), counted over the 17,898 candidates:

• False alarms: 309 candidates (1.7% of all, 1.9% of noise)
• Misses: 404 candidates (2.3% of all, 25% of pulsars)

Raising \(T\): fewer misses, more false alarms

• The classes overlap: both errors are nonzero for every \(T\)

Decision rule \(\delta(\mathbf{x})\): an action for each \(\mathbf{x}\)

• Binary: decide 1 (pulsar) or 0 (noise)

Loss \(c(y, a)\): cost of action \(a\) when the class is \(y\)

• Row: class, column: action
• \(c_{01} = 1\): one wasted telescope observation per false alarm
• \(c_{10} = 10\): a missed pulsar, valued at 10 such observations

Risk:

Example, counted over the candidates:

• \(T = 80\): \(R = 0.017 + 10 \times 0.023 = 0.24\)
• \(T = 86.4\): \(R = 0.21\), the smallest

Bayes rule \(\delta^*\): the rule with the smallest risk

Risk by iterated expectations:

• The inner term depends only on the action at that \(\mathbf{x}\)
• Choosing the smallest inner term at every \(\mathbf{x}\) minimizes the risk

Bayes rule:

Binary, with \(\pi(\mathbf{x}) = P(Y=1 \mid \mathbf{x})\):

• Expected loss of deciding noise: \(c_{10}\, \pi(\mathbf{x})\)
• Expected loss of deciding pulsar: \(c_{01}\, (1 - \pi(\mathbf{x}))\)
• Bayes rule: the smaller of the two

\(\pi(\mathbf{x}) = \mathbb{E}[Y \mid X = \mathbf{x}]\): for \(Y \in \{0,1\}\), the posterior is the MMSE predictor

Example: a Gaussian fitted to each class, \(\mathcal{N}(116.6, 17.5^2)\) for noise and \(\mathcal{N}(56.7, 30.0^2)\) for pulsars, with \(P(Y{=}1) = 0.092\)

• Expected losses equal at \(x = 90.1\)
• Minimum-risk threshold from counting: 86.4

0-1 loss: every error costs 1, so the risk is the probability of error

Expected loss of deciding class \(a\): \(1 - P(Y=a \mid \mathbf{x})\)

MAP rule: the class with the largest posterior

• Binary: decide 1 when \(\pi(\mathbf{x}) > 1/2\)

Bayes error: the error rate of the MAP rule

• No rule has a lower error rate
• Counterpart of \(\mathbb{E}[\text{Var}(Y \mid X)]\), the minimum MSE

Example: decide pulsar when \(x < 73.8\), where the fitted \(\pi(x) = 1/2\)

• 3.4% of candidates misclassified
• Deciding noise for every candidate: 9.2%
• 30% of pulsars decided as noise

Decide pulsar when its expected loss is smaller:

Threshold on the posterior:

• Equal costs: \(\tau^* = 1/2\), the MAP rule
• Costlier misses: \(\tau^* < 1/2\)

Example, \(c_{01} = 1\), \(c_{10} = 10\): \(\tau^* = 1/11 = 0.091\)

Odds of an event: \(\text{odds}(A) = \dfrac{P(A)}{P(A^c)} = \dfrac{P(A)}{1 - P(A)}\)

• \(P = 0.092\): odds \(0.10\), written 1:9.9

Likelihood ratio:

• How many times more probable \(\mathbf{x}\) is under class 1 than under class 0
• Depends on the class densities only: no prior, no costs

Bayes theorem, odds form (\(p(\mathbf{x})\) cancels):

Log form: \(\log\) posterior odds \(= \log \Lambda(\mathbf{x}) + \log\) prior odds

• \(\log \Lambda(\mathbf{x})\) is the change from prior to posterior log-odds
• Independent observations: their \(\log \Lambda\) values add

Example at \(x = 80\), fitted densities:

• \(\Lambda = 0.00983 / 0.00256 = 3.84\)
• Posterior odds \(= 3.84 \times 1/9.9 = 0.39\), \(\pi = 0.28\)
• \(\Lambda > 1\): \(x\) is more probable under the pulsar class
• \(\pi < 1/2\): noise remains the more probable class

Decide pulsar when posterior odds \(> c_{01}/c_{10}\)

Substituting posterior odds \(= \Lambda(\mathbf{x}) \times\) prior odds:

• \(\Lambda(\mathbf{x})\) depends on the data only
• \(\eta\) depends on the costs and the prior only

Example: \(\eta = \dfrac{1}{10} \cdot \dfrac{16{,}259}{1{,}639} = 0.99\)

• Prior factor \(P(Y{=}0)/P(Y{=}1) = 9.9\)
• Cost factor \(c_{01}/c_{10} = 1/10\)
• For the pulsar data the two nearly cancel

Reject option: a third action, defer the decision to a human reviewer, at cost \(c_r\)

Expected loss of each action:

• Decide noise: \(c_{10}\,\pi\)
• Decide pulsar: \(c_{01}(1-\pi)\)
• Defer: \(c_r\), for every \(\pi\)

Bayes rule, the smallest of the three:

• Noise: \(\pi < c_r / c_{10}\)
• Pulsar: \(\pi > 1 - c_r / c_{01}\)
• Defer: in between

Two thresholds only if \(c_r < \dfrac{c_{01}\,c_{10}}{c_{01} + c_{10}}\), otherwise a single threshold, \(\tau^*\)

Example, \(c_{01} = 1\), \(c_{10} = 10\), \(c_r = 0.3 < 10/11 = 0.91\), two thresholds:

• Noise: \(\pi < 0.03\)
• Pulsar: \(\pi > 0.70\)
• Defer: \(0.03 \le \pi \le 0.70\)

Gaussian fits to the pulsar data (spreads 17.5 and 30.0): log-odds \(= \log \Lambda(x) +\) log prior odds

• Quadratic in \(x\), not linear

General: classes \(\mathcal{N}(\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)\)

Shared \(\boldsymbol{\Sigma}\): the quadratic terms cancel

Different variances, one feature: a quadratic term remains

• Linear in the features \((x, x^2)\)

Data: labeled pairs \((\mathbf{x}_i, y_i)\), not densities and not \(\pi(\mathbf{x})\)

Estimating \(\pi(\mathbf{x})\) from labeled pairs:

1. Estimate \(p(\mathbf{x} \mid Y=k)\) and \(P(Y=k)\), then apply Bayes theorem
2. Fit \(\mathbf{w}\) and \(b\) in \(\pi(\mathbf{x}) = 1/(1 + e^{-(\mathbf{w}^T\mathbf{x} + b)})\): logistic regression

Least squares on 0/1 labels estimates \(\mathbb{E}[Y \mid X = x] = P(Y = 1 \mid x)\)

• The right quantity: the MMSE predictor

Linear model \(\hat\pi(x) = w\,x + b\):

• No bounds: below 0 at one end, above 1 at the other
• One slope: it cannot be flat in the tails and steep in between

Example: least-squares line through the 17,898 pulsar labels

• \(\hat\pi(x) = 0.932 - 0.0076\,x\)
• Negative for \(x > 123\), a third of the candidates
• 0.78 at \(x = 20\), where the observed pulsar fraction is 1.00

Odds \(\dfrac{\pi}{1 - \pi}\): from 0 to \(\infty\)

• \(\pi = 0.5\): odds 1 (1:1)
• \(\pi = 0.75\): odds 3 (3:1)
• \(\pi = 0.25\): odds 1/3 (1:3)

Log-odds \(\log \dfrac{\pi}{1 - \pi}\): any real number

• \(\pi = 0.5\): 0
• \(\pi = 0.75\): \(+1.10\), \(\pi = 0.25\): \(-1.10\)
• Each real value matches exactly one \(\pi\) in \((0, 1)\)

Example: observed pulsar fraction per bin of \(x\) (figure)

• Bounded on the probability scale
• Close to a line on the log-odds scale from 55 to 140

Linear score \(z = \mathbf{w}^T\mathbf{x} + b\), which can be any real number

Model: the score is the log-odds

Inverting: odds \(= e^{z}\), so \(\pi = e^{z}/(1 + e^{z})\)

Sigmoid:

• Between 0 and 1 for every \(z\)
• \(\pi = 1/2\) at \(z = 0\)
• Flat in both tails, steepest at \(z = 0\)

Not least squares on log-odds: a label’s log-odds is \(\log(1/0)\) or \(\log(0/1)\)

• \(\mathbf{w}\) and \(b\) come from maximum likelihood

Example, fitted to the pulsar data: \(w = -0.111\), \(b = 7.97\)

• \(x = 80\): \(z = -0.87\), \(\pi = 0.29\)
• Steepest at \(x = 72.1\), slope \(-0.028\) per unit
• Least-squares line: slope \(-0.0076\) everywhere

Log-odds \(= \log \Lambda(\mathbf{x}) + \log\) prior odds

• Linear in \(\mathbf{x}\) whenever \(\log \Lambda(\mathbf{x})\) is linear

\(\log \Lambda\) linear in \(\mathbf{x}\):

• Gaussian classes, shared covariance
• Poisson counts, rates \(\lambda_0, \lambda_1\): \(\log \Lambda = x \log \frac{\lambda_1}{\lambda_0} - (\lambda_1 - \lambda_0)\)
• Independent binary features, \(P(x_j = 1 \mid Y = k) = \theta_{kj}\):

Word counts, as in spam filtering: independent Poisson counts per word give \(\log \Lambda\) linear in the count vector

• Logistic regression on counts has the form these class models imply

Two alternatives to the sigmoid:

• Probit: \(\pi = \Phi(\mathbf{w}^T\mathbf{x} + b)\), the Gaussian CDF
• Complementary log-log: \(\pi = 1 - \exp(-e^{\mathbf{w}^T\mathbf{x} + b})\)

Fit: nearly the same

• Probit, rescaled, stays within 0.0095 of the sigmoid
• Complementary log-log: asymmetric, approaching 1 faster than 0

Log-odds only:

• Evidence and prior add: log-odds \(= \log \Lambda(\mathbf{x}) + \log\) prior odds
• \(e^{w_j}\): the factor on the odds per unit of \(x_j\)
• Probit and log-log scores have no such reading

Reason for the sigmoid: this reading, not its fit

Estimating \(\pi(\mathbf{x})\), with \(d\) input features:

1. Class-density model: estimate the class densities and prior, then apply Bayes theorem
2. Logistic regression: fit \(\mathbf{w}\) and \(b\) in \(\pi(\mathbf{x}) = \sigma(\mathbf{w}^T\mathbf{x} + b)\)

Parameters, Gaussian classes with a shared covariance:

• Class densities: means \(2d\), covariance \(d(d+1)/2\), prior 1
• Logistic regression: \(d + 1\)

Assumptions:

• Class densities: Gaussian classes with one covariance
• Logistic regression: log-odds linear in \(\mathbf{x}\)

Limitations:

• Class densities: wrong densities give a wrong posterior, even when the log-odds is linear
• Logistic regression: no model of \(p(\mathbf{x})\), and it needs more samples when the Gaussian model is correct

Log-odds \(= \sum_j w_j x_j + b\)

• \(w_j\): change in log-odds when \(x_j\) rises by 1, others fixed
• \(e^{w_j}\): the factor on the odds for that increase
• \(b\): log-odds at \(\mathbf{x} = \mathbf{0}\)

Example, pulsar fit \(w = -0.111\):

• Each unit of \(x\) multiplies the odds of a pulsar by \(e^{-0.111} = 0.90\)
• Each 10 units: by \(e^{-1.11} = 0.33\)

Decision boundary \(\pi(\mathbf{x}) = \tau\): \(\;\mathbf{w}^T\mathbf{x} + b = \log \frac{\tau}{1 - \tau}\)

• A hyperplane in \(\mathbf{x}\)
• A different \(\tau\) gives a parallel hyperplane: same \(\mathbf{w}\), shifted
• Example: \(\tau = 1/2\) at \(x = 72.1\), \(\tau = 1/11\) at \(x = 92.9\)

Hard threshold on the score: decide 1 if \(z > 0\)

• The Bayes rule at \(\tau = 1/2\), written on the log-odds
• Output 0 or 1, slope 0 everywhere except the jump at \(z = 0\)

Clamp \(0.5 + z/4\), cut to \([0, 1]\)

• Slope 0 outside \(|z| < 2\)

Sigmoid \(\sigma(z)\): rises from 0 to 1 around \(z = 0\)

• Output a probability, not a label
• Thresholding \(\sigma(z)\) at \(\tau\) recovers the decision
• Slope positive at every \(z\): any change in \(\mathbf{w}\) or \(b\) changes \(\pi\)

Fitting by gradient needs that slope: with the threshold or the clamp, most changes in \(\mathbf{w}\) leave the output unchanged

Chain rule on \(\sigma(z) = (1 + e^{-z})^{-1}\):

• Computed from \(\sigma(z)\) alone
• Largest at \(z = 0\): \(\sigma'(0) = 1/4\)
• Near 0 for large \(|z|\)

Example: the slope of \(\pi(x)\) is \(w\,\pi(1 - \pi)\), largest in size at \(x = 72.1\), where it is \(-0.111 \times 1/4 = -0.028\)

Each label is a Bernoulli draw given its input:

• \(y_i = 1\): \(\pi_i\)
• \(y_i = 0\): \(1 - \pi_i\)

Independent labels: the likelihood is a product

Negative log: the product becomes a sum

• Each term: the binary cross-entropy of one label
• Maximizing \(L\) is the same as minimizing \(\mathcal{L}\)

In floating point, the product underflows to 0 and the sum does not

• Pulsar labels under the fitted model: the product is exactly 0 from \(n = 4{,}232\)

Example: four labeled points \((x_i, y_i)\)

• \((0, 0)\), \((1, 1)\), \((1.5, 0)\), \((2.5, 1)\)
• Overlap: \(y = 1\) at \(x = 1\), \(y = 0\) at \(x = 1.5\)

At \(w = 2\), \(b = -2\):

\(\mathcal{L} = 2.18\)

• The point \((1.5, 0)\), predicted at \(\pi = 0.73\), contributes 60% of the loss
• \(-\log(1 - \pi) \to \infty\) as \(\pi \to 1\): no limit on the cost of a confident mistake

Pulsar fit: the worst-predicted 1% of candidates carry 37% of \(\mathcal{L}\)

• Pulsars at \(x\) from 130 to 139, given \(\pi \approx 0.001\): loss up to 7.4 each

One label: \(z = \mathbf{w}^T\mathbf{x} + b\), \(\pi = \sigma(z)\)

Chain rule:

Weights and bias:

• Prediction error times input
• Confident and right: near 0
• Confident and wrong: near \(\pm\mathbf{x}\)

Squared error on \(\pi\) keeps the factor \(\pi(1 - \pi)\), which is near 0 when the prediction is confidently wrong (figure)

Stack the data: \(\mathbf{X} \in \mathbb{R}^{n \times d}\) with rows \(\mathbf{x}_i^T\), \(\boldsymbol{\pi} = \sigma(\mathbf{X}\mathbf{w} + b\mathbf{1})\)

Sum over labels:

Least squares, same form:

• Both: (prediction − label), weighted by the inputs
• Least squares predicts \(\mathbf{X}\mathbf{w}\); logistic predicts \(\sigma(\mathbf{X}\mathbf{w} + b\mathbf{1})\)
• Cost: \(O(nd)\) per gradient in both

Example, four points at \(w = 2\), \(b = -2\) (figure):

• \(\partial \mathcal{L} / \partial w = 0.48 > 0\): lowering \(w\) lowers the loss
• \(\partial \mathcal{L} / \partial b = 0.30 > 0\): lowering \(b\) lowers the loss
• The point \((1.5, 0)\) contributes most of both

Stationary point:

• \(\mathbf{w}\) appears inside \(\sigma\), so the equations are nonlinear in \(\mathbf{w}\)
• Least squares: \(\mathbf{X}^T\mathbf{X}\mathbf{w} = \mathbf{X}^T\mathbf{y}\) is linear and is solved in one step
• Logistic: solved by iteration

The \(b\) equation: \(\sum_i \pi_i = \sum_i y_i\)

• At the fit, the predicted number of cases equals the observed number

Example, four points: \(w = 1.60\), \(b = -2.00\)

• \(\sum_i \pi_i = 2.00 = \sum_i y_i\)
• Pulsar fit: \(\sum_i \pi_i = 1{,}639 = \sum_i y_i\)

Zero start: \(\mathbf{w} = \mathbf{0}\), \(b = 0\)

• \(\pi = 1/2\) for every input
• \(\mathcal{L} = n \log 2\): 2.77 for the four points

Class-fraction start: \(\mathbf{w} = \mathbf{0}\), \(b = \log \dfrac{\bar y}{1 - \bar y}\)

• \(\pi = \bar y\), the fraction of class-1 labels, for every input
• Pulsar data, 9.2% pulsars: \(b = -2.29\)
• Starts closer to the answer on imbalanced data

One minimum: the loss is convex, so every start ends at the same \((\mathbf{w}, b)\)

• The start sets the number of steps, not the answer

Random starts: needed for networks, where units that start equal stay equal

One iteration:

1. Predict: \(\pi_i = \sigma(\mathbf{w}^T\mathbf{x}_i + b)\)
2. Differentiate: \(\mathbf{g}_w = \mathbf{X}^T(\boldsymbol{\pi} - \mathbf{y})\), \(\;g_b = \sum_i (\pi_i - y_i)\)
3. Update: \(\mathbf{w} \leftarrow \mathbf{w} - \eta\,\mathbf{g}_w\), \(\;b \leftarrow b - \eta\,g_b\)
4. Stop when \(\|\mathbf{g}\|\) is small; otherwise repeat

Example, four points, zero start, \(\eta = 0.5\):

• Iteration 1: \(\mathbf{g} = (-1.0, 0)\) gives \(w = 0.5\), \(b = 0\), and \(\mathcal{L}\) falls from 2.77 to 2.56
• 190 iterations to \(\|\mathbf{g}\| < 10^{-6}\): \(w = 1.60\), \(b = -2.00\)

Step size: \(\eta < 2/\lambda_{\max}\) of the Hessian

• At the zero start: \(2 / 3.11 = 0.64\)
• \(\eta = 1\): the loss oscillates and stays above the minimum

Hessian over \(\boldsymbol{\theta} = (\mathbf{w}, b)\), with \(\tilde{\mathbf{X}} = [\mathbf{X} \;\; \mathbf{1}]\):

Positive semidefinite: for any \(\mathbf{v}\)

• Positive definite when \(\tilde{\mathbf{X}}\) has full column rank

Example at the zero start: \(\mathbf{H} = \begin{bmatrix} 2.375 & 1.25 \\ 1.25 & 1 \end{bmatrix}\), eigenvalues 3.11 and 0.26

A minimum exists only if no hyperplane separates the classes

• Separable: \(\mathcal{L} \to 0\) as \(\|\mathbf{w}\| \to \infty\), a value never reached
• Pulsar subsamples of \(n = 50\): 33% separable or without a pulsar, so no minimum

Newton: \(\boldsymbol{\theta} \leftarrow \boldsymbol{\theta} - \mathbf{H}^{-1}\mathbf{g}\), with \(\mathbf{g} = \tilde{\mathbf{X}}^T(\boldsymbol{\pi} - \mathbf{y})\), \(\mathbf{H} = \tilde{\mathbf{X}}^T\mathbf{S}\tilde{\mathbf{X}}\)

Rearranged:

• The normal equations, each point weighted by \(\pi_i(1 - \pi_i)\)
• New weights every step: iteratively reweighted least squares

Cost per step:

• Newton: \(O(nd^2)\) to form \(\mathbf{H}\), \(O(d^3)\) to solve
• Gradient descent: \(O(nd)\)

Example, pulsar data, zero start, steps to \(\hat w\) within a relative \(10^{-6}\):

• Newton, \(x\) unscaled: 7 steps
• Gradient descent, \(x\) standardized: 477 steps
• Gradient descent, \(x\) unscaled: condition number of \(\mathbf{H}\) is \(2.6 \times 10^5\)
• After 200,000 steps at the largest stable \(\eta\): \(w = -0.053\) against \(-0.111\)

Reference: the fit to all 17,898 pulsar candidates, \(\hat w = -0.111\)

Subsamples: 500 random subsets at each \(n\), each fit by Newton

• Spread shrinks roughly as \(1/\sqrt{n}\)
• \(n = 50\): about 5 pulsars per subset
• 33% of those subsets are separable or have no pulsar, so no finite estimate

Subsets of one data set stand in for new samples from the same survey

Data: 5,574 human-labeled English text messages

• 747 spam (13.4%), 4,827 legitimate
• 4,000 training messages, 1,574 held-out test messages

Legitimate:

• “When are you guys leaving?”
• “Sorry, I’ll call later in meeting”

Spam:

• “Santa Calling! Would your little ones like a call from Santa Xmas eve? Call 09058094583 to book your time.”
• “Free Msg: get Gnarls Barkleys “Crazy” ringtone TOTALLY FREE just reply GO to this message right now!”

Decision: deliver the message, or send it to the spam folder

• Class 1: spam
• False alarm: a legitimate message sent to the spam folder
• Miss: a spam message delivered

Costs: a lost legitimate message is worse than a delivered spam

• \(c_{01} = 10\), \(c_{10} = 1\)
• The reverse of the pulsar costs, where a miss was costlier

Model: logistic regression on the words of each message

Data: T. A. Almeida, J. M. Gómez Hidalgo, A. Yamakami, SMS Spam Collection, UCI Machine Learning Repository (2011), CC BY 4.0.

Vocabulary: every word that appears in at least 5 of the 4,000 training messages

• \(d = 1{,}471\) words

Features: \(x_j\) = number of times word \(j\) appears in the message

• No word order
• The median message has 10 nonzero entries out of 1,471

def tokenize(message):
    message = re.sub(r"&\w+;", " ", message)       # HTML escapes in the raw text
    return re.findall(r"[a-z0-9£]+", message.lower())

def count_vector(message):
    x = np.zeros(len(vocab))
    for word in tokenize(message):
        if word in column:           # word -> column index
            x[column[word]] += 1
    return x

Example: “Free Msg: get Gnarls Barkleys “Crazy” ringtone TOTALLY FREE just reply GO to this message right now!”

x = count_vector('Free Msg: get Gnarls Barkleys "Crazy" ringtone '
                 'TOTALLY FREE just reply GO to this message right now!')
print({vocab[j]: int(x[j]) for j in np.flatnonzero(x)})

{'crazy': 1, 'free': 2, 'get': 1, 'go': 1, 'just': 1, 'message': 1, 'msg': 1, 'now': 1, 'reply': 1, 'right': 1, 'ringtone': 1, 'this': 1, 'to': 1}

• 13 nonzero counts
• “gnarls”, “barkleys”, “totally”: not in the vocabulary

The loop from the four-point example, with the gradient averaged over the \(n = 4{,}000\) messages:

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

n = len(y_train)
w = np.zeros(len(vocab))
b = np.log(y_train.mean() / (1 - y_train.mean()))   # class-fraction start
history = []
for step in range(1, 3001):
    p = sigmoid(X_train @ w + b)
    w -= 0.5 * X_train.T @ (p - y_train) / n
    b -= 0.5 * (p - y_train).mean()
    if step % 25 == 0:
        p_test = sigmoid(X_test @ w + b)
        history.append((step, ((p > 0.5) == y_train).mean(),
                        ((p_test > 0.5) == y_test).mean(), np.linalg.norm(w)))

Newton: a \(1{,}472 \times 1{,}472\) Hessian to form and solve at every step

• Test accuracy: 98.3% after 3,000 steps
• \(\|\mathbf{w}\|\) still growing: the training messages are close to separable

Odds factor \(e^{w_j}\): each occurrence of word \(j\) multiplies the odds of spam by it

Largest factors:

• “txt” 8.8, “call” 7.0, “uk” 6.2, “stop” 6.0, “text” 5.7
• Words from sign-up and reply instructions

Smallest factors:

• “i” 0.26, “me” 0.36, “my” 0.38
• First-person words, common in messages between people

Bias \(b = -4.55\): odds of spam 0.011 for a message with no vocabulary words

Score: \(z = b + \sum_j w_j x_j\), one term per word occurrence

Three test messages:

• Legitimate message: “call” and “free” raise \(\pi\) to 0.21, still delivered
• Spam without spam vocabulary: \(\pi = 0.04\), missed
• Word counts only: no feature for meaning

Threshold from the costs, \(c_{01} = 10\), \(c_{10} = 1\):

Test messages (1,574, of which 200 spam):

• \(\tau^*\): 22 more spam delivered, 5 fewer legitimate messages lost

Lowest test cost: 37, at \(\tau = 0.69\), not at \(\tau^*\)

• 18 test messages have \(0.69 < \pi \le 0.91\), mean \(\pi = 0.80\)
• 17 of them are spam (94%): the fitted \(\pi\) is too low in this range
• Choosing \(\tau = 0.69\) from this curve tunes on the test set
• Threshold selection uses a validation split, not the test set

Spam test set: 1,574 messages, 200 spam (12.7%)

Accuracy, fraction correct:

• Deliver every message: 87.3%, and no spam caught
• Fitted model at \(\tau = 1/2\): 98.3%
• On the 4,000 training messages: 99.3%, higher than on held-out messages

Confusion matrix, \(\tau = 1/2\):

• 27 errors of two kinds: 22 spam delivered, 5 legitimate messages lost
• Accuracy: each error counts 1
• Cost: 1 per spam delivered, 10 per legitimate message lost

Four counts:

• TP: spam sent to spam
• FN: spam delivered (miss)
• FP: legitimate message sent to spam (false alarm)
• TN: legitimate message delivered

At \(\tau^* = 0.91\): TP 156, FN 44, FP 0, TN 1,374

• Accuracy falls to 97.2%
• Cost falls from 72 to 44

Accuracy and cost disagree whenever the two errors cost different amounts

Two rates, one per true class:

• TPR: fraction of spam caught (recall)
• FPR: fraction of legitimate messages lost
• Each divides by one class’s size: unchanged by the spam fraction

At \(\tau = 1/2\): TPR \(= 178/200 = 0.89\), FPR \(= 5/1{,}374 = 0.0036\)

ROC curve: (FPR, TPR) as \(\tau\) runs from 1 to 0

Area under the curve (AUC): 0.979

• Probability that a random spam scores above a random legitimate message
• Depends only on the ordering of the \(\pi\) values

Precision: fraction of messages sent to spam that are spam

• At \(\tau = 1/2\): \(178/183 = 0.97\)

Same TPR and FPR, rarer spam, by Bayes theorem:

• At 0.1% spam, legitimate messages sent to spam outnumber the spam caught 4 to 1
• The ROC curve does not change with the spam fraction

Calibrated: among messages given \(\pi \approx q\), a fraction \(q\) is spam

Test messages, grouped by \(\pi\):

• Close in the two outer groups, 1,505 of the 1,574 messages
• Too low from 0.5 to 0.9: mean \(\pi\) 0.72, observed 0.81
• The messages that put the lowest test cost below \(\tau^*\) are in this range
• Total: \(\sum_i \pi_i = 196\) against 200 spam

Separate from ranking: rescaling \(\pi\) without reordering changes calibration and leaves the ROC curve unchanged

Log-loss: the negative log-likelihood per message

• Measures ranking and calibration together
• On the training messages: the quantity gradient descent minimizes

Spam model, by gradient step:

• Training loss keeps falling as \(\|\mathbf{w}\|\) grows
• Test loss is lowest near 0.063, around step 7,000, then rises
• After step 7,000, further steps fit the training messages, not new ones

Posterior: a vector \(\boldsymbol{\pi}(\mathbf{x}) = [\pi_1(\mathbf{x}), \ldots, \pi_K(\mathbf{x})]^T\), with \(\pi_k(\mathbf{x}) = P(Y = k \mid \mathbf{x})\)

• Nonnegative entries that sum to 1: a point on the probability simplex

Cost matrix \(\mathbf{C} \in \mathbb{R}^{K \times K}\): \(c_{ka}\) = cost of deciding \(a\) when the class is \(k\)

Expected loss of deciding \(a\):

Bayes rule: \(\delta^*(\mathbf{x}) = \arg\min_a R(a \mid \mathbf{x})\)

• 0-1 loss: \(\arg\max_k \pi_k(\mathbf{x})\), the MAP rule
• No single threshold: straight lines on the simplex separate the decisions

Example, \(K = 3\), class 3 costly to miss:

• Class 3 is decided on 76% of the simplex, against 33% under 0-1 loss

One-vs-rest: \(K\) binary logistic models, class \(k\) against the other \(K - 1\)

• Reuses the two-class model
• Each fit separately: nothing ties the \(K\) outputs together

Example data: 342 penguins, three species

• 151 Adelie, 68 Chinstrap, 123 Gentoo
• \(\mathbf{x}\) = bill length and flipper length (mm)

One-vs-rest on these birds:

• Sum of the three outputs: 0.28 to 1.97 across the birds
• 24% of the birds: sum outside 0.9 to 1.1

Argmax of the \(K\) outputs still gives a label

Expected losses \(\mathbf{C}^T\boldsymbol{\pi}\) need a probability vector

Data: Palmer Penguins, A. M. Horst, A. P. Hill, K. B. Gorman, 2020 (CC0); collected by K. Gorman and Palmer Station LTER

Softmax of scores \(\mathbf{z} \in \mathbb{R}^K\):

• Exponentiating: every entry positive
• Normalizing: entries sum to 1
• Order kept: a larger score gives a larger probability
• One score far above the rest: output near one-hot

Example, \(\mathbf{z} = [2, 1, -1]\):

• \(e^{\mathbf{z}} = [7.39, 2.72, 0.37]\), sum 10.48
• \(\operatorname{softmax}(\mathbf{z}) = [0.705, 0.259, 0.035]\)

Scores tripled, \(\mathbf{z} = [6, 3, -3]\): \([0.953, 0.047, 0.0001]\)

• Larger gaps, closer to one-hot

Log of each numerator of Bayes theorem, classes \(\mathcal{N}(\boldsymbol{\mu}_k, \boldsymbol{\Sigma})\):

Normalizing: the shared term cancels

• \(K = 2\): \(\pi_2 = \sigma(z_2 - z_1)\), the two-class logistic posterior
• Boundaries between classes: straight lines, \(z_k = z_j\)

Example: one Gaussian per species with a shared covariance, fit to the 240 training penguins

• 99 of the 102 test birds decided correctly

Model: one weight vector and one bias per class

• \(\mathbf{W} \in \mathbb{R}^{K \times d}\), row \(k\) is \(\mathbf{w}_k^T\); \(\mathbf{b} \in \mathbb{R}^K\)
• \(K(d + 1)\) parameters, fit by maximum likelihood; penguins: 9

Pairwise log-odds, linear:

• Each pair of classes: a two-class logistic model
• Boundary between \(k\) and \(j\): the hyperplane \(z_k = z_j\)
• Region of class \(k\): the \(K - 1\) half-spaces \(z_k > z_j\), intersected, a convex set

Assumption: every pairwise log-odds is linear in \(\mathbf{x}\); no model of \(p(\mathbf{x} \mid Y = k)\)

Example, flipper 205 mm: Adelie below a bill length of 44.1 mm, Gentoo to 47.6 mm, Chinstrap above

One label \(y \in \{1, \ldots, K\}\), one-hot \(\mathbf{y}\): the categorical cross-entropy, the negative log-likelihood of a discrete label

Softmax substituted:

• Second term, the log-sum-exp: at least \(\max_k z_k\)
• One wrong class far ahead: \(\ell \approx z_{\max} - z_y\), linear in the gap

Data set: \(\mathcal{L}(\mathbf{W}, \mathbf{b}) = \sum_{i=1}^{n} \ell_i\)

\(K = 2\), \(z_1 = 0\), labels \(y \in \{0, 1\}\): \(\ell = \log(1 + e^{z}) - y\,z\), the binary cross-entropy

Example: a test penguin, Chinstrap, bill 50.8 mm, flipper 210 mm

• Decided Gentoo: the Gentoo score is 0.58 higher
• \(\ell = -1.35 + \log(e^{-3.28} + e^{1.35} + e^{1.93}) = 1.03 = -\log 0.357\)

Overflow: \(e^{z}\) exceeds the largest float for

• \(z > 709.78\) in float64
• \(z > 88.72\) in float32

Shift by \(m = \max_k z_k\), which changes nothing:

• Largest term becomes \(e^0 = 1\), so the sum lies between 1 and \(K\)
• Terms that underflow to 0 are below \(10^{-308}\) next to that 1

Loss from the scores: \(\ell = -z_y + m + \log \sum_k e^{z_k - m}\)

• Computing \(\pi_y\) first and then \(-\log \pi_y\) fails when \(\pi_y\) underflows to 0

\(K = 2\): \(\ell = \log(1 + e^{z}) - y\,z\) is np.logaddexp(0, z) - y * z

Differentiating \(\ell = -z_y + \log \sum_k e^{z_k}\):

• Through the softmax Jacobian \(\mathbf{J} = \operatorname{diag}(\boldsymbol{\pi}) - \boldsymbol{\pi}\boldsymbol{\pi}^T\): \(\mathbf{J}^T(-\mathbf{y} \oslash \boldsymbol{\pi}) = \boldsymbol{\pi} - \mathbf{y}\), the \(1/\pi\) cancels as in the two-class case

Weights, \(\mathbf{z} = \mathbf{W}\mathbf{x} + \mathbf{b}\):

Data set, one row per sample, \(\boldsymbol{\Pi}, \mathbf{Y} \in \mathbb{R}^{n \times K}\):

• \(K = 2\) with a reference class: the two-class \(\mathbf{X}^T(\boldsymbol{\pi} - \mathbf{y})\)
• \(O(ndK)\) per gradient

Example, the test Chinstrap, \(\mathbf{y} = [0, 1, 0]\), standardized \(\mathbf{x} = [1.26, 0.59]\):

• \(\boldsymbol{\pi} - \mathbf{y} = [0.003,\ -0.643,\ 0.640]\)
• Row \(k\) of \(\nabla_{\mathbf{W}}\,\ell\) is \((\pi_k - y_k)\,\mathbf{x}^T\):

Descent step on this bird:

• Chinstrap row moves toward \(\mathbf{x}\): its score here rises
• Gentoo row moves away from \(\mathbf{x}\): its score here falls
• Each class’s error \(\pi_k - y_k\) sets the size of its row’s change

Standardized features, averaged gradient:

W, b = np.zeros((3, 2)), np.zeros(3)
for step in range(20000):
    P = softmax(X_train @ W.T + b)       # n x K
    W -= 0.5 * (P - Y_train).T @ X_train / n
    b -= 0.5 * (P - Y_train).mean(axis=0)

Penguins, 240 training and 102 test birds:

• Converged: gradient norm below \(10^{-9}\), training loss 0.126 per bird
• 99 of 102 test birds correct

Test confusion matrix, rows the true species:

• \(K(K - 1) = 6\) kinds of error; all 3 errors are Chinstraps

Separable: a hyperplane (in one dimension, a threshold) with every training point on the side of its class

Example: 50 pulsar candidates, 6 of them pulsars

• Every pulsar: profile mean at most 79.1
• Every noise candidate: at least 83.1

Sharper threshold, lower loss: multiplying \(\mathbf{w}\) and \(b\) by \(c > 1\) keeps the boundary and steepens \(\pi(x)\)

• Every point’s loss decreases
• \(\mathcal{L} \to 0\) as \(\|\mathbf{w}\| \to \infty\): no finite minimizer
• In the limit, \(\pi(x)\) is a step: 0 or 1 at every training point
• Gradient descent: \(|w|\) = 0.8, 5.4, 28.3 after 10, 1,000, 100,000 steps

Spam training messages: close to separable, \(\|\mathbf{w}\| = 22.0\) at step 20,000 and still growing

Penalized loss, bias not penalized:

• Same penalty as ridge regression

Unique: Hessian \(\frac{1}{n}\tilde{\mathbf{X}}^T\mathbf{S}\tilde{\mathbf{X}} + \lambda\,\operatorname{diag}(1, \ldots, 1, 0)\) is positive definite

Finite: \(R_\lambda \to \infty\) as \(\|\mathbf{w}\| \to \infty\), also for separable data

Gradient step:

• Factor \(1 - \eta\lambda\) on \(\mathbf{w}\) every step: weight decay

Example, the separable subset, \(\lambda = 0.01\): \(w = -2.60\) at step 1,000 and at step 100,000

MAP estimate:

Gaussian prior, independent \(w_j \sim \mathcal{N}(0, \tau^2)\):

• L2 penalty, \(\lambda = 1/(n\tau^2)\)
• Ridge regression: \(\lambda = \sigma^2/(n\tau^2)\), \(\sigma^2\) the noise variance

Laplace prior, independent \(p(w_j) = \frac{1}{2s}e^{-|w_j|/s}\):

• L1 penalty \(\lambda_1\|\mathbf{w}\|_1\), \(\lambda_1 = 1/(ns)\)

Narrower prior: larger penalty

One weight, loss \(\tfrac{1}{2}(w - a)^2\) around the unpenalized value \(a\):

Penalty slope at \(w = 0\):

• L2: \(\lambda w = 0\), so any nonzero loss slope moves \(\hat w\) off 0
• L1: \(\pm\lambda\), so a loss slope smaller than \(\lambda\) leaves \(\hat w\) at 0

Many weights: \(\hat w_j = 0\) whenever \(\lvert \partial(\text{loss}) / \partial w_j \rvert \le \lambda\) at the solution

Proximal gradient: a gradient step, then the soft threshold on every weight

def soft(v, t):
    return np.sign(v) * np.maximum(np.abs(v) - t, 0)

p = sigmoid(X_train @ w + b)
w = soft(w - eta * X_train.T @ (p - y_train) / n, eta * lam)

Spam model, 5,000 steps:

• Positive at \(\lambda = 3 \times 10^{-3}\): “txt”, “call”, “uk”, “stop”, “text”
• Negative: “i”, “me”, “my”

Subgradient steps, \(\lambda\operatorname{sign}(\mathbf{w})\) added to the gradient: no weight exactly 0 after 5,000 steps

5-fold cross-validation on the 4,000 training messages:

• For each \(\lambda\): fit on 4 of 5 parts, log-loss on the fifth
• Average over the 5 choices of held-out part
• Cost: 5 fits per \(\lambda\), 45 fits for 9 values

Spam model, L2: lowest average log-loss 0.058, at \(\lambda = 10^{-4}\)

• Refit on all 4,000, then one test evaluation
• Test log-loss 0.063, accuracy 98.2%
• Equal to the best test loss of the unpenalized run, without using the test set

Validation split: 800 of the 4,000 training messages; gradient descent on the other 3,200

Stopping rule: the step with the lowest validation loss

• Spam model: step 6,200, validation log-loss 0.048
• One test evaluation: log-loss 0.065, accuracy 98.2%

Relation to L2: from \(\mathbf{w} = \mathbf{0}\), \(\|\mathbf{w}\|\) grows with the number of steps

• Stopping early keeps \(\|\mathbf{w}\|\) smaller: 15.2 at step 6,200
• The step count is tuned in place of \(\lambda\)
• One run, no refits, 800 fewer training messages

XOR: class 1 where \(x_1\) and \(x_2\) have opposite signs

• The two class-1 quadrants lie on a diagonal
• A hyperplane splits the plane into two half-planes
• No line puts both class-1 quadrants on one side: any line misclassifies at least one quadrant

Logistic fit on 400 points:

• Weights near 0, \(\pi\) between 0.49 and 0.62 everywhere
• 54% correct

Every model in this lecture has this limit

• Logistic regression: one hyperplane
• Softmax regression: \(K\) linear scores, hyperplane boundaries

Added feature \(x_1 x_2\): logistic regression on \((x_1, x_2, x_1 x_2)\)

• A hyperplane in the three features separates the classes: 99% correct
• Still linear in the parameters, as for polynomial regression

Choosing the feature: \(x_1 x_2\) comes from knowing the pattern

• Without that knowledge: every product of inputs up to degree \(k\)

Number of features, all monomials of degree 1 to \(k\) in \(d\) inputs: \(\binom{d + k}{k} - 1\)

Neuron: a weighted sum plus a bias, then a nonlinearity

• Logistic regression: one sigmoid neuron, \(a = \pi(\mathbf{x})\)
• Softmax regression: a layer of \(K\) neurons on the same input

Fitting the neuron, as in this lecture:

• Loss: cross-entropy, the negative log-likelihood
• Gradient: \(\partial \ell / \partial z = \pi - y\), \(\;\partial \ell / \partial \mathbf{w} = (\pi - y)\,\mathbf{x}\)
• Update: a gradient-descent step on \(\mathbf{w}\) and \(b\)

Label at the output: \(y\) enters the gradient through \(\pi - y\)

Hidden layer: \(m\) logistic units on the input

Output unit: logistic regression on \(\mathbf{h}\)

• The features \(\mathbf{h}\) replace the hand-chosen products; \(\mathbf{W}_1\) is fit with \(\mathbf{w}_2\)
• \((d + 1)m + m + 1\) parameters: \(d = 784\), \(m = 100\) gives 78,601, against 308,504 degree-2 features

XOR with \(m = 2\), weights set by hand (figure):

• \(h_1 = \sigma(10(x_1 + x_2 + 1))\), \(\;h_2 = \sigma(10(1 - x_1 - x_2))\)
• \(\pi = \sigma(10(h_1 + h_2 - 1.5))\)
• 95.5% of the 400 points correct
• Weights not scaled by 10: \(\pi\) between 0.43 and 0.49 at every corner

Output weights \(\mathbf{w}_2\): logistic regression on \(\mathbf{h}\)

• Gradient \((\pi - y)\,\mathbf{h}\), with the label at the output

Hidden weights \(\mathbf{W}_1\): no label at \(\mathbf{h}\), so no \(\pi - y\) there

• A change in \(\mathbf{W}_1\) changes \(\mathbf{h}\), then the score \(\mathbf{w}_2^T\mathbf{h} + b_2\), then \(\pi\), then \(\ell\)
• The factors along that path: \(\pi - y\), the weights \(\mathbf{w}_2\), and \(\sigma' = \sigma(1 - \sigma)\)

Cost of the hidden layer:

• Loss not convex in \(\mathbf{W}_1\): several minima, and the start matters
• All-zero start: every hidden unit has the same output and gradient, and they stay identical