Stochastic Gradient Descent Code Check for Least Squares I have an R-based implementation of the gradient descent and am trying to also get it to work as SGD. The function matches R's lm function when using the entire data set. But, when I sample from the full data I am nowhere near matching even with different values for learning rates. I fear I have a bug or a misunderstanding.
R Code for implementation is:
sgd <- function(x, y, alpha, n = nrow(x), R = 500) {
    samp <- sample(nrow(x), n)
    y <- y[samp]
    y <- as.matrix(y)
    x <- x[samp,]
    M <- length(y)
    theta <- t(as.matrix(numeric(ncol(X))))
    result <- matrix(0, nrow=R, ncol = ncol(x))
    for (i in 1:R) {
        shuffle <- sample(1:nrow(x), nrow(x), replace = FALSE)
        x <- x[shuffle,]
        y <- y[shuffle]
        xtx <- crossprod(x)
        xty <- crossprod(x,y)
        gradient <- (1/M) * (xtx %*% t(theta)  - xty)
        theta <- theta - alpha * t(gradient)  
        result[i,] <- theta
    }
    list(coef = result[R,])
}

Here is some toy code to run. The run using all 10,000 cases matches lm. The second version with a subsample of 1000 and a learning rate of .01 (or any other value) does not match.
Is there an obvious error in my code or perhaps in my implementation of the concept such that it is not working on samples of data?
N <- 10000
y = rnorm(N)
x1 = rnorm(N)
x2 = rnorm(N)
x3 = rnorm(N)
res <- lm(y ~ x1+x2+x3)
X <- model.matrix(res)

sgd(X, y, 1, 10000) ## Run on full data sample (regular gradient descent)
coef(res)

sgd(X, y, .01, 1000) ## Run on data sample

 A: Here's a quick and dirty implementation showing coefficients estimated using the normal equation, using a manual implementation of gradient descent, and using a manual implementation of stochastic gradient descent:
set.seed(12345)

n_obs <- 1000
df <- data.frame(x1=rnorm(n_obs), x2=runif(n_obs), epsilon=rnorm(n_obs))
df$y <- 5 + 2 * df$x1 + 5 * df$x2 + df$epsilon

X <- as.matrix(cbind(rep(1, n_obs), df[, c("x1", "x2")]))
colnames(X) <- c("constant", "x1", "x2")

beta_hat_norm_eq <- solve(t(X) %*% X) %*% t(X) %*% df$y

## This should be (nearly) identical to coefficients(lm(y ~ x1 + x2, data=df))
beta_hat_norm_eq

## Manual gradient descent (very quick and dirty)
learning_rate <- 0.1
n_iter <- 1000
beta_hat_gd <- matrix(0, nrow=n_iter, ncol=ncol(X))
for(iter in seq_len(n_iter - 1)) {
    residuals <- df$y - X %*% beta_hat_gd[iter, ]
    gradient <- -(2 / n_obs) * t(X) %*% residuals
    beta_hat_gd[iter + 1, ] <- beta_hat_gd[iter, ] - learning_rate * gradient
}

## This should be close to beta_hat_norm_eq
beta_hat_gd[n_iter, ]

## Manual stochastic gradient descent (also _very_ quick and dirty)
learning_rate <- 0.1
batch_size <- 32
n_iter <- 1000
beta_hat_sgd <- matrix(0, nrow=n_iter, ncol=ncol(X))
for(iter in seq_len(n_iter - 1)) {
    ## Keep it simple: sample a random subset of batch_size rows on every iteration
    row_idx <- sample(seq_len(n_obs), size=batch_size)
    residuals <- df$y[row_idx] - X[row_idx, ] %*% beta_hat_sgd[iter, ]
    gradient <- -(2 / batch_size) * t(X[row_idx, ]) %*% residuals
    beta_hat_sgd[iter + 1, ] <- beta_hat_sgd[iter, ] - learning_rate * gradient
}

beta_hat_sgd[n_iter, ]

plot(beta_hat_sgd[, c(2, 3)], type="l", col="red", xlab="beta_1", ylab="beta_2")
lines(beta_hat_gd[, c(2, 3)], type="l", col="blue", lty=2)
points(x=beta_hat_norm_eq[2], y=beta_hat_norm_eq[3], col="black", pch=4, cex=3)


The normal equation coefficients are the black X, and gradient descent and stochastic gradient descent are shown in blue and red, respectively.  They all end up at roughly the same place.
