We have many good discussions about perfect separation in logistic regression. Such as, Logistic regression in R resulted in perfect separation (Hauck-Donner phenomenon). Now what? and Logistic regression model does not converge .

I personally still feel it is not intuitive for why it will be a problem and why adding regularization will fix it. I made some animations and think it will be helpful. So post his question and answer it by myself to share with the community.


A 2D demo with toy data will be used to explain what was happening for perfect separation on logistic regression with and without regularization. The experiments started with an overlapping data set and we gradually move two classes apart. The objective function contour and optima (logistic loss) will be shown in the right sub figure. The data and the linear decision boundary are plotted in left sub figure.

First we try the logistic regression without regularization.

  • As we can see with the data moving apart, the objective function (logistic loss) is changing dramatically, and the optim is moving away to a larger value.
  • When we have completed the operation, the contour will not be a "closed shape". At this time, the objective function will always be smaller when the solution moves to upper right comer.

enter image description here

Next we try logistic regression with L2 regularization (L1 is is similar).

  • With the same setup, adding a very small L2 regularization will change the objective function changes respect to the separation of the data.

  • In this case, we will always have the "convex" objective. No matter how much separation the data has.

enter image description here

code (I also use same code for this answer: Regularization methods for logistic regression)

d=mlbench::mlbench.2dnormals(100, 2, r=1)

x = d$x
y = ifelse(d$classes==1, 1, 0)

logistic_loss <- function(w){
  p    = plogis(x %*% w)
  L    = -y*log(p) - (1-y)*log(1-p)
  LwR2 = sum(L) + lambda*t(w) %*% w

logistic_loss_gr <- function(w){
  p = plogis(x %*% w)
  v = t(x) %*% (p - y)
  return(c(v) + 2*lambda*w)

w_grid_v = seq(-10, 10, 0.1)
w_grid   = expand.grid(w_grid_v, w_grid_v)

lambda = 0
opt1   = optimx::optimx(c(1,1), fn=logistic_loss, gr=logistic_loss_gr, method="BFGS")
z1     = matrix(apply(w_grid,1,logistic_loss), ncol=length(w_grid_v))

lambda = 5
opt2   = optimx::optimx(c(1,1), fn=logistic_loss, method="BFGS")
z2     = matrix(apply(w_grid,1,logistic_loss), ncol=length(w_grid_v))

plot(d, xlim=c(-3,3), ylim=c(-3,3))
abline(0, -opt1$p2/opt1$p1, col='blue',  lwd=2)
abline(0, -opt2$p2/opt2$p1, col='black', lwd=2)
contour(w_grid_v, w_grid_v, z1, col='blue',  lwd=2, nlevels=8)
contour(w_grid_v, w_grid_v, z2, col='black', lwd=2, nlevels=8, add=T)
points(opt1$p1, opt1$p2, col='blue',  pch=19)
points(opt2$p1, opt2$p2, col='black', pch=19)
  • 2
    $\begingroup$ This is very reminiscent of viscosity solutions for hyperbolic PDEs. I wonder if it is possible to directly get the limiting $\lambda\to 0^+$ solution without needing to approximate $\lambda=\epsilon$, as it is for some PDEs? $\endgroup$ – GeoMatt22 Oct 13 '16 at 16:19
  • 3
    $\begingroup$ These visualizations are fantastic. $\endgroup$ – Matthew Drury May 8 '17 at 4:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.