# What are some reasons iteratively reweighted least squares would not converge when used for logistic regression?

I've been using the glm.fit function in R to fit parameters to a logistic regression model. By default, glm.fit uses iteratively reweighted least squares to fit the parameters. What are some reasons this algorithm would fail to converge, when used for logistic regression?

• I don't know how to judge "most common", since what's common for people working in one area might be uncommon for another. [One reason (of many possible) could be complete separation - where along some linear combination of predictors all the 0's are either above or below all the 1's. You can sometimes see when it happens because at least one parameter will tend to head off toward infinity.] – Glen_b Jan 13 '15 at 19:50
• @Glen_b: Thanks for your comment, I'll change it to "some reasons." – Jessica Jan 13 '15 at 19:57

In case the two classes are separable, iteratively reweighted least squares (IRLS) would break. In such a scenario, any hyperplane that separates the two classes is a solution and there are infinitely many of them. IRLS is meant to find a maximum likelihood solution. Maximum likelihood does not have a mechanism to favor any of these solutions over the other (e.g. no concept of maximum margin). Depending on the initialization, IRLS should go toward one of these solutions and would break due to numerical problems (don't know the details of IRLS; an educated guess).

Another problem arises in case of linear-separability of the training data. Any of the hyperplane solutions corresponds to a heaviside function. Therefore, all the probabilities are either 0 or 1. The linear regression solution would be a hard classifier rather than a probabilistic classifier.

To clarify using mathematical symbols, heaviside function is $\lim_{|\mathbf{w}| \rightarrow \infty}\sigma(\mathbf{w}^T x + b)$, the limit of sigmoid function, where $\sigma$ is the sigmoid function and $(\mathbf{w}, b)$ determines the hyperplane solution. So IRLS theroretically does not stop and goes toward a $\mathbf{w}$ with increasing magnitude but would break in practice due to numerical problems.

On top of linear separation (in which the MLE is at the boundary of the parameter space), the Fisher Scoring procedure in R is not completely numerically stable. It takes steps of fixed size, which in certain pathological cases can lead to non-convergence (when the true MLE is indeed an interior point).

For example,

y <- c(1,1,1,0)
x <- rep(1,4)
fit1 <- glm.fit(x,y, family=binomial(link="logit"),start=-1.81)


yields a coefficient of $2 \times 10^{15}$ rather than the expected logit$(3/4) \approx 1.0986$.

The CRAN package glm2 provides a drop-in replacement for glm.fit that adjusts step size to ensure monotone convergence.