# Why using Newton's method for logistic regression optimization is called iterative re-weighted least squares?

Why using Newton's method for logistic regression optimization is called iterative re-weighted least squares?

It seems not clear to me because logistic loss and least squares loss are completely different things.

• I don't think they are the same. IRLS is Newton-Raphson with the expected Hessian rather than the observed Hessian. Commented May 3, 2018 at 20:44
• @DimitriyV.Masterov thanks, could you tell me more on Expected Hessian vs Observed? Also, what do you think about this explanation Commented May 3, 2018 at 20:46
• Commented Jun 14, 2019 at 12:19
• I would better say, not exactly: IR in general is a broader term, can reweight iteratively by any optimization algorithm till either tolerance_condition (for flat function) or minimization of convex_function. Newton's methods use both gradient & hessian and thus are appropriate exactly for convex optimization. And IRLS yes, really have above-mentioned implementation Commented Mar 15 at 17:12

Summary: GLMs are fit via Fisher scoring which, as Dimitriy V. Masterov notes, is Newton-Raphson with the expected Hessian instead (i.e. we use an estimate of the Fisher information instead of the observed information). If we are using the canonical link function it turns out that the observed Hessian equals the expected Hessian so NR and Fisher scoring are the same in that case. Either way, we'll see that Fisher scoring is actually fitting a weighted least squares linear model, and the coefficient estimates from this converge* on a maximum of the logistic regression likelihood. Aside from reducing fitting a logistic regression to an already solved problem, we also get the benefit of being able to use linear regression diagnostics on the final WLS fit to learn about our logistic regression.

I'm going to keep this focused on logistic regression, but for a more general perspective on maximum likelihood in GLMs I recommend section 15.3 of this chapter which goes through this and derives IRLS in a more general setting (I think it's from John Fox's Applied Regression Analysis and Generalized Linear Models).

$^*$ see comments at the end

## The likelihood and score function

We will be fitting our GLM by iterating something of the form $$b^{(m+1)} = b^{(m)} - J^{-1}_{(m)}\nabla \ell(b^{(m)})$$ where $\ell$ is the log likelihood and $J_{m}$ will be either the observed or expected Hessian of the log likelihood.

Our link function is a function $g$ that maps the conditional mean $\mu_i = E(y_i | x_i)$ to our linear predictor, so our model for the mean is $g(\mu_i) = x_i^T\beta$. Let $h$ be the inverse link function mapping the linear predictor to the mean.

For a logistic regression we have a Bernoulli likelihood with independent observations so $$\ell(b; y) = \sum_{i=1}^n y_i\log h(x_i^T b) + (1 - y_i) \log(1 - h(x_i^Tb)).$$ Taking derivatives, $$\frac{\partial \ell}{\partial b_j} = \sum_{i=1}^n \frac{y_i}{h(x_i^T b)} h'(x_i^T b) x_{ij} - \frac{1 - y_i}{1 - h(x_i^T b)} h'(x_i^T b) x_{ij}$$ $$= \sum_{i=1}^n x_{ij} h'(x_i^T b) \left(\frac{y_i}{h(x_i^T b)} - \frac{1 - y_i}{1 - h(x_i^T b)} \right)$$ $$= \sum_i x_{ij} \frac{h'(x_i^T b)}{h(x_i^T b)(1 - h(x_i^T b))}(y_i - h(x_i^T b)).$$

Now let's suppose we're using the canonical link function $g_c = \text{logit}$. Then $g^{-1}_c(x) := h_c(x) = \frac{1}{1+e^{-x}}$ so $h_c' = h_c \cdot (1-h_c)$ which means this simplifies to $$\frac{\partial \ell}{\partial b_j} = \sum_i x_{ij} (y_i - h_c(x_i^T b))$$ so $$\nabla \ell (b; y) = X^T (y - \hat y).$$ Furthermore, still using $h_c$, $$\frac{\partial^2 \ell}{\partial b_k \partial b_j} = - \sum_i x_{ij} \frac{\partial}{\partial b_k} h_c(x_i^T b) = - \sum_i x_{ij}x_{ik} \left[h_c(x_i^T b) (1 - h_c(x_i^T b))\right].$$

Let $$W = \text{diag}\left(h_c(x_1^T b)(1 - h_c(x_1^T b)), \dots, h_c(x_n^T b)(1 - h_c(x_n^T b))\right) = \text{diag}\left(\hat y_1(1 - \hat y_1), \dots, \hat y_n (1 - \hat y_n)\right).$$ Then we have $$H = -X^TWX$$ and note how this doesn't have any $y_i$ in it anymore, so $E(H) = H$ (we're viewing this as a function of $b$ so the only random thing is $y$ itself). Thus we've shown that Fisher scoring is equivalent to Newton-Raphson when we use the canonical link in logistic regression. Also by virtue of $\hat y_i \in (0,1)$ $-X^TWX$ will always be strictly negative definite, although numerically if $\hat y_i$ gets too close to $0$ or $1$ then we may have weights round to $0$ which can make $H$ negative semidefinite and therefore computationally singular.

Now create the working response $z = W^{-1}(y - \hat y)$ and note that $$\nabla \ell = X^T(y - \hat y) = X^T W z.$$

All together this means that we can optimize the log likelihood by iterating $$b^{(m+1)} = b^{(m)} + (X^T W_{(m)} X)^{-1}X^T W_{(m)} z_{(m)}$$ and $(X^T W_{(m)} X)^{-1}X^T W_{(m)} z_{(m)}$ is exactly $\hat \beta$ for a weighted least squares regression of $z_{(m)}$ on $X$.

Checking this in R:

set.seed(123)
p <- 5
n <- 500
x <- matrix(rnorm(n * p), n, p)
betas <- runif(p, -2, 2)
hc <- function(x) 1 /(1 + exp(-x)) # inverse canonical link
p.true <- hc(x %*% betas)
y <- rbinom(n, 1, p.true)

# fitting with our procedure
my_IRLS_canonical <- function(x, y, b.init, hc, tol=1e-8) {
change <- Inf
b.old <- b.init
while(change > tol) {
eta <- x %*% b.old  # linear predictor
y.hat <- hc(eta)
h.prime_eta <- y.hat * (1 - y.hat)
z <- (y - y.hat) / h.prime_eta

b.new <- b.old + lm(z ~ x - 1, weights = h.prime_eta)$coef # WLS regression change <- sqrt(sum((b.new - b.old)^2)) b.old <- b.new } b.new } my_IRLS_canonical(x, y, rep(1,p), hc) # x1 x2 x3 x4 x5 # -1.1149687 2.1897992 1.0271298 0.8702975 -1.2074851 glm(y ~ x - 1, family=binomial())$coef
# x1         x2         x3         x4         x5
# -1.1149687  2.1897992  1.0271298  0.8702975 -1.2074851


and they agree.

Now if we're not using the canonical link we don't get the simplification of $\frac{h'}{h(1-h)} = 1$ in $\nabla \ell$ so $H$ becomes much more complicated, and we therefore see a noticeable difference by using $E(H)$ in our Fisher scoring.

Here's how this will go: we already worked out the general $\nabla \ell$ so the Hessian will be the main difficulty. We need $$\frac{\partial^2 \ell}{\partial b_k \partial b_j} = \sum_i x_{ij} \frac{\partial}{\partial b_k}h'(x_i^T b) \left(\frac{y_i}{h(x_i^T b)} - \frac{1 - y_i}{1 - h(x_i^T b)} \right)$$ $$= \sum_i x_{ij}x_{ik} \left[h''(x_i^T b) \left(\frac{y_i}{h(x_i^T b)} - \frac{1 - y_i}{1 - h(x_i^T b)} \right) - h'(x_i^T b)^2\left(\frac{y_i}{h(x_i^T b)^2} + \frac{1-y_i}{(1-h(x_i^T b))^2} \right)\right]$$

Via the linearity of expectation all we need to do to get $E(H)$ is replace each occurrence of $y_i$ with its mean under our model which is $\mu_i=h(x_i^T\beta)$. Each term in the summand will therefore contain a factor of the form $$h''(x_i^T b) \left(\frac{h(x_i^T \beta)}{h(x_i^T b)} - \frac{1 - h(x_i^T \beta)}{1 - h(x_i^T b)} \right) - h'(x_i^T b)^2\left(\frac{h(x_i^T \beta)}{h(x_i^T b)^2} + \frac{1-h(x_i^T \beta)}{(1-h(x_i^T b))^2} \right).$$ But to actually do our optimization we'll need to estimate each $\beta$, and at step $m$ $b^{(m)}$ is the best guess we have. This means that this will reduce to $$h''(x_i^T b) \left(\frac{h(x_i^T b)}{h(x_i^T b)} - \frac{1 - h(x_i^T b)}{1 - h(x_i^T b)} \right) - h'(x_i^T b)^2\left(\frac{h(x_i^T b)}{h(x_i^T b)^2} + \frac{1-h(x_i^T b)}{(1-h(x_i^T b))^2} \right)$$ $$= - h'(x_i^T b)^2\left(\frac{1}{h(x_i^T b)} + \frac{1}{1-h(x_i^T b)} \right)$$ $$= -\frac{h'(x_i^T b)^2}{h(x_i^T b)(1-h(x_i^T b))}.$$ This means we will use $J$ with $$J_{jk} = -\sum_i x_{ij}x_{ik} \frac{h'(x_i^T b)^2}{h(x_i^T b)(1-h(x_i^T b))}.$$

Now let $$W^* = \text{diag}\left(\frac{h'(x_1^T b)^2}{h(x_1^T b)(1-h(x_1^T b))} ,\dots, \frac{h'(x_n^T b)^2}{h(x_n^T b)(1-h(x_n^T b))}\right)$$ and note how under the canonical link $h_c' = h_c \cdot (1-h_c)$ reduces $W^*$ to $W$ from the previous section. This lets us write $$J = -X^TW^*X$$ except this is now $\hat E(H)$ rather than necessarily being $H$ itself, so this can differ from Newton-Raphson. For all $i$ $W_{ii}^* > 0$ so aside from numerical issues $J$ will be negative definite.

We have $$\frac{\partial \ell}{\partial b_j} = \sum_i x_{ij} \frac{h'(x_i^T b)}{h(x_i^T b)(1 - h(x_i^T b))}(y_i - h(x_i^T b))$$ so letting our new working response be $z^* = D^{-1}(y-\hat y)$ with $D=\text{diag}\left(h'(x_1^T b), \dots, h'(x_n^T b)\right)$, we have $\nabla \ell = X^TW^*z^*$.

All together we are iterating $$b^{(m+1)} = b^{(m)} + (X^T W_{(m)}^* X)^{-1}X^T W_{(m)}^* z_{(m)}^*$$ so this is still a sequence of WLS regressions except now it's not necessarily Newton-Raphson.

I've written it out this way to emphasize the connection to Newton-Raphson, but frequently people will factor the updates so that each new point $b^{(m+1)}$ is itself the WLS solution, rather than a WLS solution added to the current point $b^{(m)}$. If we wanted to do this, we can do the following: $$b^{(m+1)} = b^{(m)} + (X^T W_{(m)}^* X)^{-1}X^T W_{(m)}^* z_{(m)}^*$$ $$= (X^T W_{(m)}^* X)^{-1}\left(X^T W_{(m)}^* Xb^{(m)}+ X^TW^*_{(m)}z_{(m)}^* \right)$$ $$= (X^T W_{(m)}^* X)^{-1}X^TW_{(m)}^*\left(Xb^{(m)}+ z_{(m)}^* \right)$$ so if we're going this way you'll see the working response take the form $\eta^{(m)} + D^{-1}_{(m)}(y - \hat y^{(m)})$, but it's the same thing.

Let's confirm that this works by using it to perform a probit regression on the same simulated data as before (and this is not the canonical link, so we need this more general form of IRLS).

my_IRLS_general <- function(x, y, b.init, h, h.prime, tol=1e-8) {
change <- Inf
b.old <- b.init
while(change > tol) {
eta <- x %*% b.old  # linear predictor
y.hat <- h(eta)
h.prime_eta <- h.prime(eta)
w_star <- h.prime_eta^2 / (y.hat * (1 - y.hat))
z_star <- (y - y.hat) / h.prime_eta

b.new <- b.old + lm(z_star ~ x - 1, weights = w_star)$coef # WLS change <- sqrt(sum((b.new - b.old)^2)) b.old <- b.new } b.new } # probit inverse link and derivative h_probit <- function(x) pnorm(x, 0, 1) h.prime_probit <- function(x) dnorm(x, 0, 1) my_IRLS_general(x, y, rep(0,p), h_probit, h.prime_probit) # x1 x2 x3 x4 x5 # -0.6456508 1.2520266 0.5820856 0.4982678 -0.6768585 glm(y~x-1, family=binomial(link="probit"))$coef
# x1         x2         x3         x4         x5
# -0.6456490  1.2520241  0.5820835  0.4982663 -0.6768581


and again the two agree.

Finally, a few quick comments on convergence (I'll keep this brief as this is getting really long and I'm no expert at optimization). Even though theoretically each $J_{(m)}$ is negative definite, bad initial conditions can still prevent this algorithm from converging. In the probit example above, changing the initial conditions to b.init=rep(1,p) results in this, and that doesn't even look like a suspicious initial condition. If you step through the IRLS procedure with that initialization and these simulated data, by the second time through the loop there are some $\hat y_i$ that round to exactly $1$ and so the weights become undefined. If we're using the canonical link in the algorithm I gave we won't ever be dividing by $\hat y_i (1 - \hat y_i)$ to get undefined weights, but if we've got a situation where some $\hat y_i$ are approaching $0$ or $1$, such as in the case of perfect separation, then we'll still get non-convergence as the gradient dies without us reaching anything.

• +1. I love how detailed your answers often are. Commented May 5, 2018 at 20:01
• You stated "the coefficient estimates from this converge on a maximum of the logistic regression likelihood." Is that necessarily so, from any initial values? Commented May 5, 2018 at 20:05
• @MarkL.Stone ah I was being too casual there, didn't mean to offend the optimization people :) I'll add some more details (and would appreciate your thoughts on them when I do)
– jld
Commented May 5, 2018 at 20:12
• any chance you watched the link I posted? Seems that video is talking from machine learning perspective, just optimize logistic loss, without talking about Hessain expectation? Commented May 7, 2018 at 15:06
• @hxd1011 in that pdf i linked to (link again: sagepub.com/sites/default/files/upm-binaries/…) on page 24 of it the author goes into the theory and explains what exactly makes a link function canonical. I found that pdf extremely helpful when I first came across this (although it took me a while to get through).
– jld
Commented May 7, 2018 at 15:36

Because error is being calculated according Mean-Squared-Error being minimized - as is Least-Squared technique/approach, used iteratively in Newton-method.

in Python Newtons-method looks like (extracted from here):

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import optimize
from scipy.special import expit

X = df[['gre', 'gpa', 'rank']].values

# intercept column (for standartized data no intercept needed)
X= np.column_stack((X, np.ones(X.shape[0])))

## sigmoid(x):
def logistic_model(weights, data):
"""Logistic model."""
return expit(data.dot(weights))

def log_likelihood(weights, data, labels):
"""Normalized - negative log likelihood."""
# OR CROSS-ENTROPY:
# !!!!! calc. acc. mse(w, x, y)= np.mean((f(w, x) - y) ** 2), where f(w, x) is Hypothesis
scores = np.dot(data, weights)
return -np.mean(labels*scores - np.log(1 + np.exp(scores)))

predictions = logistic_model(weights, data)
return -data.T.dot(labels - predictions)/len(data)

def log_likelihood_hessian(weights, data, labels):
"""Hessian of negative log likelihood."""
predictions = logistic_model(weights, data)
diag = np.diag(predictions * (1 - predictions))
return data.T.dot(diag.dot(data))/len(data)

def newton_method(init, steps, grad, hessian, num_to_keep=None):
xs = [init]
# Iteratively !!!!!
for step in steps:
Hinv = np.linalg.pinv(hessian(xs[-1]))
if num_to_keep:
xs = xs[-num_to_keep:]
return xs

objective = lambda w: log_likelihood(w, X, y)
hessian = lambda w: log_likelihood_hessian(w, X, y)
init = np.zeros(X.shape[1])

ws_newton = newton_method(init, np.ones(10), gradient, hessian)

print(ws_newton[-1])
plt.plot([log_likelihood(w, X, y) for w in ws_newton], label="cross-entropy error")
_ = plt.legend()
plt.title("Newton-method stepwise cost reduction")
plt.show()


final thetas: [ 2.29395950e-03 7.77013574e-01 -5.60031387e-01 -3.44954840e+00]

IRLS can not use 2nd derivative, only LS for Logistic-Regression, like here

P.S. GLM with IRLS - Newton's implementations with QR & SVD decomposition