Gradient Descent in R
This R script illustrates the difference between stochastic gradient descent (SGD) and gradient descent (GD) in a simple linear regression problem.
-
Synthetic data generation
The dataset is simulated according to the model
y = X θ* + ε
where X is a 100×2 matrix of Gaussian random variables, θ* = (2, −1), and ε is Gaussian noise with standard deviation 0.5.
-
Loss surface
A grid of values for θ₁ and θ₂ is created around the true parameters.
For each pair (θ₁, θ₂), the script computes the sum of squared residuals
L(θ₁, θ₂) = Σ (xᵢᵀ θ − yᵢ)².
A contour plot displays the loss surface in the (θ₁, θ₂) plane.
-
Stochastic Gradient Descent (SGD)
Starting from θ = (0,0), the algorithm updates the parameters after each observation:
θ ← θ − δ × 2xᵢ(xᵢᵀθ − yᵢ),
with learning rate δ = 0.05.
This is repeated for 10 epochs, randomly shuffling the data at each pass.
The parameter trajectories (θ₁, θ₂) are stored and plotted.
SGD produces a noisy trajectory but usually converges quickly toward the minimum.
-
Gradient Descent (GD)
Here, the gradient is computed over the whole dataset at each iteration:
θ ← θ − δ × 2Xᵀ(Xθ − y),
with δ = 0.001 and 50 iterations.
This yields a much smoother and more regular trajectory than SGD.
-
Visual comparison
The final contour plot shows both optimization paths:
– red line: stochastic gradient descent
– black line: batch gradient descent
Both methods move toward the same minimum (close to θ* = (2, −1)), but SGD fluctuates around the valley of the loss function, whereas GD follows a smooth, deterministic path.
This example demonstrates how SGD and GD explore the parameter space differently while aiming to minimize the same quadratic loss function.
set.seed(1974)
N <- 100
X <- matrix(rnorm(N * 2), ncol = 2)
theta_true <- c(2, -1)
y <- as.vector(X %*% theta_true + rnorm(N, sd = 0.5))
t1_seq <- seq(theta_true[1] - 2, theta_true[1] + 2, length.out = 120)
t2_seq <- seq(theta_true[2] - 2, theta_true[2] + 2, length.out = 120)
L <- matrix(0, nrow = length(t1_seq), ncol = length(t2_seq))
for (i in 1:length(t1_seq)) for (j in 1:length(t2_seq)) L[i, j] <- sum((X %*% c(t1_seq[i],t2_seq[j]) - y)^2)
contour(t1_seq, t2_seq, L,
nlevels = 20,
xlab = expression(theta[1]),
ylab = expression(theta[2]))
epoch <- 10
theta.sgd <- matrix(0,N*epoch+1,2)
delta <- 0.05
for (k in 1:epoch)
{
indi <- sample(1:N,N)
for (i in 1:N)
{
grad <- 2*X[indi[i],]*as.vector(X[indi[i],]%*%theta.sgd[i+(k-1)*N,]-y[indi[i]])
theta.sgd[i+(k-1)*N+1,] <- theta.sgd[i+(k-1)*N,]-delta*grad
}
}
plot(1:(N*epoch+1),theta.sgd[,1],type="l")
plot(1:(N*epoch+1),theta.sgd[,2],type="l")
epoch <- 50
theta.gd <- matrix(0,epoch+1,2)
delta <- 0.001
for (k in 1:epoch)
{
grad <- 2*t(X)%*%(X%*%theta.gd[k,]-y)
theta.gd[k+1,] <- theta.gd[k,]-delta*grad
}
plot(1:(epoch+1),theta.gd[,1],type="l")
plot(1:(epoch+1),theta.gd[,2],type="l")
contour(t1_seq, t2_seq, L,
nlevels = 20,
xlab = expression(theta[1]),
ylab = expression(theta[2]))
lines(theta.sgd[,1],theta.sgd[,2],lwd=2,col="red")
lines(theta.gd[,1],theta.gd[,2],lwd=2,col="black")