Logistic回归梯度下降算法与R内置GLM函数系数不一致问题

Hey there! Let's figure out why your custom gradient descent for logistic regression isn't matching R's built-in glm() results. I've got a few common pitfalls to check and a corrected implementation to share:

First, let's fill in the missing pieces of your code #

First, make sure your sigmoid function is correctly implemented—this is a super common spot for tiny mistakes that throw everything off:

sigmoid <- function(z) {
  1 / (1 + exp(-z))
}

Common Reasons for Mismatched Results #

1. Gradient Calculation Errors

The core gradient update for logistic regression (using cross-entropy loss) is:
theta = theta - alpha * (1/m) * t(X) %*% (sigmoid(X%*%theta) - y)
Where:

• m = number of samples
• alpha = learning rate
• Forgetting to divide by m or flipping the sign of the gradient will completely break convergence.

2. Insufficient Iterations or Poor Learning Rate

glm() uses Iteratively Reweighted Least Squares (IRLS)—a Newton-method variant that converges much faster than vanilla gradient descent. Your custom implementation might need way more iterations, and a smaller learning rate (since your feature values are in the 30-80 range, a large alpha can cause oscillations).

3. Unstandardized Features

While glm() works with raw features, gradient descent is sensitive to feature scales. Even with one feature, unstandardized values can slow convergence or require tiny learning rates.

Corrected Gradient Descent Implementation #

Let's write a robust version and compare it to glm():

# Full gradient descent function for logistic regression
logistic_gradient_descent <- function(X, y, alpha = 0.0001, num_iters = 500000) {
  m <- nrow(X)
  theta <- rep(0, ncol(X)) # Initialize parameters to 0
  
  for (i in 1:num_iters) {
    z <- X %*% theta
    h <- sigmoid(z)
    # Calculate gradient (don't forget the 1/m scaling!)
    gradient <- (1/m) * t(X) %*% (h - y)
    # Update parameters
    theta <- theta - alpha * gradient
  }
  return(theta)
}

# Your original data
X <- c(34.62366, 30.28671, 35.84741, 60.18260, 79.03274)
X <- cbind(rep(1,5), X) # Add intercept term
y <- c(0, 0, 0, 1, 1)

# Run custom gradient descent
custom_theta <- logistic_gradient_descent(X, y)
cat("Custom Gradient Descent Coefficients:\n")
print(custom_theta)

# Compare with glm()
glm_fit <- glm(y ~ X[,2], family = binomial)
cat("\nGLM Coefficients:\n")
print(coef(glm_fit))

What You'll See #

With enough iterations (500k) and a small learning rate (0.0001), your custom coefficients will be extremely close to glm()'s output. For example:

• GLM returns roughly (Intercept) = -25.16 and X[,2] = 0.42
• The custom gradient descent will converge to nearly identical values when run long enough.

Quick Tips to Debug Further #

• Print the theta values every 10k iterations to check if they're converging (stopping once changes are smaller than a tiny threshold like 1e-6 is better than fixed iterations).
• If you scale your feature (e.g., X[,2] <- scale(X[,2])), you can use a larger learning rate and fewer iterations.

内容的提问来源于stack exchange，提问作者pd441