FYS-STK4155/doc/LectureNotes/_build/jupyter_execute/week37.ipynb

Week 37: Gradient descent methods

Morten Hjorth-Jensen, Department of Physics, University of Oslo, Norway

Date: September 8-12, 2025

Plans for week 37, lecture Monday

Plans and material for the lecture on Monday September 8.

The family of gradient descent methods

1. Plain gradient descent (constant learning rate), reminder from last week with examples using OLS and Ridge

2. Improving gradient descent with momentum

3. Introducing stochastic gradient descent

4. More advanced updates of the learning rate: ADAgrad, RMSprop and ADAM

Readings and Videos:

1. Recommended: Goodfellow et al, Deep Learning, introduction to gradient descent, see sections 4.3-4.5 at https://www.deeplearningbook.org/contents/numerical.html and chapter 8.3-8.5 at https://www.deeplearningbook.org/contents/optimization.html

2. Rashcka et al, pages 37-44 and pages 278-283 with focus on linear regression.

3. Video on gradient descent at https://www.youtube.com/watch?v=sDv4f4s2SB8

4. Video on Stochastic gradient descent at https://www.youtube.com/watch?v=vMh0zPT0tLI

Material for lecture Monday September 8

Gradient descent and revisiting Ordinary Least Squares from last week

Last week we started with linear regression as a case study for the gradient descent methods. Linear regression is a great test case for the gradient descent methods discussed in the lectures since it has several desirable properties such as:

1. An analytical solution (recall homework sets for week 35).

2. The gradient can be computed analytically.

3. The cost function is convex which guarantees that gradient descent converges for small enough learning rates

We revisit an example similar to what we had in the first homework set. We have a function of the type

In [1]:

import numpy as np
x = 2*np.random.rand(m,1)
y = 4+3*x+np.random.randn(m,1)

with x_i \in [0,1] is chosen randomly using a uniform distribution. Additionally we have a stochastic noise chosen according to a normal distribution \cal {N}(0,1). The linear regression model is given by

h_\theta(x) = \boldsymbol{y} = \theta_0 + \theta_1 x,

such that

\boldsymbol{y}_i = \theta_0 + \theta_1 x_i.

Gradient descent example

Let \mathbf{y} = (y_1,\cdots,y_n)^T, \mathbf{\boldsymbol{y}} = (\boldsymbol{y}_1,\cdots,\boldsymbol{y}_n)^T and \theta = (\theta_0, \theta_1)^T

It is convenient to write \mathbf{\boldsymbol{y}} = X\theta where X \in \mathbb{R}^{100 \times 2} is the design matrix given by (we keep the intercept here)

X \equiv \begin{bmatrix}
1 & x_1  \\
\vdots & \vdots  \\
1 & x_{100} &  \\
\end{bmatrix}.

The cost/loss/risk function is given by

C(\theta) = \frac{1}{n}||X\theta-\mathbf{y}||_{2}^{2} = \frac{1}{n}\sum_{i=1}^{100}\left[ (\theta_0 + \theta_1 x_i)^2 - 2 y_i (\theta_0 + \theta_1 x_i) + y_i^2\right]

and we want to find \theta such that C(\theta) is minimized.

The derivative of the cost/loss function

Computing \partial C(\theta) / \partial \theta_0 and \partial C(\theta) / \partial \theta_1 we can show that the gradient can be written as

\nabla_{\theta} C(\theta) = \frac{2}{n}\begin{bmatrix} \sum_{i=1}^{100} \left(\theta_0+\theta_1x_i-y_i\right) \\
\sum_{i=1}^{100}\left( x_i (\theta_0+\theta_1x_i)-y_ix_i\right) \\
\end{bmatrix} = \frac{2}{n}X^T(X\theta - \mathbf{y}),

where X is the design matrix defined above.

The Hessian matrix

The Hessian matrix of C(\theta) is given by

\boldsymbol{H} \equiv \begin{bmatrix}
\frac{\partial^2 C(\theta)}{\partial \theta_0^2} & \frac{\partial^2 C(\theta)}{\partial \theta_0 \partial \theta_1}  \\
\frac{\partial^2 C(\theta)}{\partial \theta_0 \partial \theta_1} & \frac{\partial^2 C(\theta)}{\partial \theta_1^2} &  \\
\end{bmatrix} = \frac{2}{n}X^T X.

This result implies that C(\theta) is a convex function since the matrix X^T X always is positive semi-definite.

Simple program

We can now write a program that minimizes C(\theta) using the gradient descent method with a constant learning rate \eta according to

\theta_{k+1} = \theta_k - \eta \nabla_\theta C(\theta_k), \ k=0,1,\cdots

We can use the expression we computed for the gradient and let use a \theta_0 be chosen randomly and let \eta = 0.001. Stop iterating when ||\nabla_\theta C(\theta_k) || \leq \epsilon = 10^{-8}. Note that the code below does not include the latter stop criterion.

And finally we can compare our solution for \theta with the analytic result given by \theta= (X^TX)^{-1} X^T \mathbf{y}.

Gradient Descent Example

Here our simple example

In [2]:

%matplotlib inline


# Importing various packages
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import sys

# the number of datapoints
n = 100
x = 2*np.random.rand(n,1)
y = 4+3*x+np.random.randn(n,1)

X = np.c_[np.ones((n,1)), x]
# Hessian matrix
H = (2.0/n)* X.T @ X
# Get the eigenvalues
EigValues, EigVectors = np.linalg.eig(H)
print(f"Eigenvalues of Hessian Matrix:{EigValues}")

theta_linreg = np.linalg.inv(X.T @ X) @ X.T @ y
print(theta_linreg)
theta = np.random.randn(2,1)

eta = 1.0/np.max(EigValues)
Niterations = 1000

for iter in range(Niterations):
    gradient = (2.0/n)*X.T @ (X @ theta-y)
    theta -= eta*gradient

print(theta)
xnew = np.array([[0],[2]])
xbnew = np.c_[np.ones((2,1)), xnew]
ypredict = xbnew.dot(theta)
ypredict2 = xbnew.dot(theta_linreg)
plt.plot(xnew, ypredict, "r-")
plt.plot(xnew, ypredict2, "b-")
plt.plot(x, y ,'ro')
plt.axis([0,2.0,0, 15.0])
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'Gradient descent example')
plt.show()

Gradient descent and Ridge

We have also discussed Ridge regression where the loss function contains a regularized term given by the L_2 norm of \theta,

C_{\text{ridge}}(\theta) = \frac{1}{n}||X\theta -\mathbf{y}||^2 + \lambda ||\theta||^2, \ \lambda \geq 0.

In order to minimize C_{\text{ridge}}(\theta) using GD we adjust the gradient as follows

\nabla_\theta C_{\text{ridge}}(\theta)  = \frac{2}{n}\begin{bmatrix} \sum_{i=1}^{100} \left(\theta_0+\theta_1x_i-y_i\right) \\
\sum_{i=1}^{100}\left( x_i (\theta_0+\theta_1x_i)-y_ix_i\right) \\
\end{bmatrix} + 2\lambda\begin{bmatrix} \theta_0 \\ \theta_1\end{bmatrix} = 2 (\frac{1}{n}X^T(X\theta - \mathbf{y})+\lambda \theta).

We can easily extend our program to minimize C_{\text{ridge}}(\theta) using gradient descent and compare with the analytical solution given by

\theta_{\text{ridge}} = \left(X^T X + n\lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y}.

The Hessian matrix for Ridge Regression

The Hessian matrix of Ridge Regression for our simple example is given by

\boldsymbol{H} \equiv \begin{bmatrix}
\frac{\partial^2 C(\theta)}{\partial \theta_0^2} & \frac{\partial^2 C(\theta)}{\partial \theta_0 \partial \theta_1}  \\
\frac{\partial^2 C(\theta)}{\partial \theta_0 \partial \theta_1} & \frac{\partial^2 C(\theta)}{\partial \theta_1^2} &  \\
\end{bmatrix} = \frac{2}{n}X^T X+2\lambda\boldsymbol{I}.

This implies that the Hessian matrix is positive definite, hence the stationary point is a minimum. Note that the Ridge cost function is convex being a sum of two convex functions. Therefore, the stationary point is a global minimum of this function.

Program example for gradient descent with Ridge Regression

In [3]:

from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import sys

# the number of datapoints
n = 100
x = 2*np.random.rand(n,1)
y = 4+3*x+np.random.randn(n,1)

X = np.c_[np.ones((n,1)), x]
XT_X = X.T @ X

#Ridge parameter lambda
lmbda  = 0.001
Id = n*lmbda* np.eye(XT_X.shape[0])

# Hessian matrix
H = (2.0/n)* XT_X+2*lmbda* np.eye(XT_X.shape[0])
# Get the eigenvalues
EigValues, EigVectors = np.linalg.eig(H)
print(f"Eigenvalues of Hessian Matrix:{EigValues}")


theta_linreg = np.linalg.inv(XT_X+Id) @ X.T @ y
print(theta_linreg)
# Start plain gradient descent
theta = np.random.randn(2,1)

eta = 1.0/np.max(EigValues)
Niterations = 100

for iter in range(Niterations):
    gradients = 2.0/n*X.T @ (X @ (theta)-y)+2*lmbda*theta
    theta -= eta*gradients

print(theta)
ypredict = X @ theta
ypredict2 = X @ theta_linreg
plt.plot(x, ypredict, "r-")
plt.plot(x, ypredict2, "b-")
plt.plot(x, y ,'ro')
plt.axis([0,2.0,0, 15.0])
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'Gradient descent example for Ridge')
plt.show()

Using gradient descent methods, limitations

• Gradient descent (GD) finds local minima of our function. Since the GD algorithm is deterministic, if it converges, it will converge to a local minimum of our cost/loss/risk function. Because in ML we are often dealing with extremely rugged landscapes with many local minima, this can lead to poor performance.

• GD is sensitive to initial conditions. One consequence of the local nature of GD is that initial conditions matter. Depending on where one starts, one will end up at a different local minima. Therefore, it is very important to think about how one initializes the training process. This is true for GD as well as more complicated variants of GD.

• Gradients are computationally expensive to calculate for large datasets. In many cases in statistics and ML, the cost/loss/risk function is a sum of terms, with one term for each data point. For example, in linear regression, E \propto \sum_{i=1}^n (y_i - \mathbf{w}^T\cdot\mathbf{x}_i)^2; for logistic regression, the square error is replaced by the cross entropy. To calculate the gradient we have to sum over all n data points. Doing this at every GD step becomes extremely computationally expensive. An ingenious solution to this, is to calculate the gradients using small subsets of the data called "mini batches". This has the added benefit of introducing stochasticity into our algorithm.

• GD is very sensitive to choices of learning rates. GD is extremely sensitive to the choice of learning rates. If the learning rate is very small, the training process take an extremely long time. For larger learning rates, GD can diverge and give poor results. Furthermore, depending on what the local landscape looks like, we have to modify the learning rates to ensure convergence. Ideally, we would adaptively choose the learning rates to match the landscape.

• GD treats all directions in parameter space uniformly. Another major drawback of GD is that unlike Newton's method, the learning rate for GD is the same in all directions in parameter space. For this reason, the maximum learning rate is set by the behavior of the steepest direction and this can significantly slow down training. Ideally, we would like to take large steps in flat directions and small steps in steep directions. Since we are exploring rugged landscapes where curvatures change, this requires us to keep track of not only the gradient but second derivatives. The ideal scenario would be to calculate the Hessian but this proves to be too computationally expensive.

• GD can take exponential time to escape saddle points, even with random initialization. As we mentioned, GD is extremely sensitive to initial condition since it determines the particular local minimum GD would eventually reach. However, even with a good initialization scheme, through the introduction of randomness, GD can still take exponential time to escape saddle points.

Momentum based GD

We discuss here some simple examples where we introduce what is called 'memory'about previous steps, or what is normally called momentum gradient descent. For the mathematical details, see whiteboad notes from lecture on September 8, 2025.

Improving gradient descent with momentum

In [4]:

from numpy import asarray
from numpy import arange
from numpy.random import rand
from numpy.random import seed
from matplotlib import pyplot
 
# objective function
def objective(x):
	return x**2.0
 
# derivative of objective function
def derivative(x):
	return x * 2.0
 
# gradient descent algorithm
def gradient_descent(objective, derivative, bounds, n_iter, step_size):
	# track all solutions
	solutions, scores = list(), list()
	# generate an initial point
	solution = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])
	# run the gradient descent
	for i in range(n_iter):
		# calculate gradient
		gradient = derivative(solution)
		# take a step
		solution = solution - step_size * gradient
		# evaluate candidate point
		solution_eval = objective(solution)
		# store solution
		solutions.append(solution)
		scores.append(solution_eval)
		# report progress
		print('>%d f(%s) = %.5f' % (i, solution, solution_eval))
	return [solutions, scores]
 
# seed the pseudo random number generator
seed(4)
# define range for input
bounds = asarray([[-1.0, 1.0]])
# define the total iterations
n_iter = 30
# define the step size
step_size = 0.1
# perform the gradient descent search
solutions, scores = gradient_descent(objective, derivative, bounds, n_iter, step_size)
# sample input range uniformly at 0.1 increments
inputs = arange(bounds[0,0], bounds[0,1]+0.1, 0.1)
# compute targets
results = objective(inputs)
# create a line plot of input vs result
pyplot.plot(inputs, results)
# plot the solutions found
pyplot.plot(solutions, scores, '.-', color='red')
# show the plot
pyplot.show()

Same code but now with momentum gradient descent

In [5]:

from numpy import asarray
from numpy import arange
from numpy.random import rand
from numpy.random import seed
from matplotlib import pyplot
 
# objective function
def objective(x):
	return x**2.0
 
# derivative of objective function
def derivative(x):
	return x * 2.0
 
# gradient descent algorithm
def gradient_descent(objective, derivative, bounds, n_iter, step_size, momentum):
	# track all solutions
	solutions, scores = list(), list()
	# generate an initial point
	solution = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])
	# keep track of the change
	change = 0.0
	# run the gradient descent
	for i in range(n_iter):
		# calculate gradient
		gradient = derivative(solution)
		# calculate update
		new_change = step_size * gradient + momentum * change
		# take a step
		solution = solution - new_change
		# save the change
		change = new_change
		# evaluate candidate point
		solution_eval = objective(solution)
		# store solution
		solutions.append(solution)
		scores.append(solution_eval)
		# report progress
		print('>%d f(%s) = %.5f' % (i, solution, solution_eval))
	return [solutions, scores]
 
# seed the pseudo random number generator
seed(4)
# define range for input
bounds = asarray([[-1.0, 1.0]])
# define the total iterations
n_iter = 30
# define the step size
step_size = 0.1
# define momentum
momentum = 0.3
# perform the gradient descent search with momentum
solutions, scores = gradient_descent(objective, derivative, bounds, n_iter, step_size, momentum)
# sample input range uniformly at 0.1 increments
inputs = arange(bounds[0,0], bounds[0,1]+0.1, 0.1)
# compute targets
results = objective(inputs)
# create a line plot of input vs result
pyplot.plot(inputs, results)
# plot the solutions found
pyplot.plot(solutions, scores, '.-', color='red')
# show the plot
pyplot.show()

Overview video on Stochastic Gradient Descent (SGD)

What is Stochastic Gradient Descent There are several reasons for using stochastic gradient descent. Some of these are:

1. Efficiency: Updates weights more frequently using a single or a small batch of samples, which speeds up convergence.

2. Hopefully avoid Local Minima

3. Memory Usage: Requires less memory compared to computing gradients for the entire dataset.

Batches and mini-batches

In gradient descent we compute the cost function and its gradient for all data points we have.

In large-scale applications such as the ILSVRC challenge, the training data can have on order of millions of examples. Hence, it seems wasteful to compute the full cost function over the entire training set in order to perform only a single parameter update. A very common approach to addressing this challenge is to compute the gradient over batches of the training data. For example, a typical batch could contain some thousand examples from an entire training set of several millions. This batch is then used to perform a parameter update.

Pros and cons

1. Speed: SGD is faster than gradient descent because it uses only one training example per iteration, whereas gradient descent requires the entire dataset. This speed advantage becomes more significant as the size of the dataset increases.

2. Convergence: Gradient descent has a more predictable convergence behaviour because it uses the average gradient of the entire dataset. In contrast, SGD’s convergence behaviour can be more erratic due to its random sampling of individual training examples.

3. Memory: Gradient descent requires more memory than SGD because it must store the entire dataset for each iteration. SGD only needs to store the current training example, making it more memory-efficient.

Convergence rates

1. Stochastic Gradient Descent has a faster convergence rate due to the use of single training examples in each iteration.

2. Gradient Descent as a slower convergence rate, as it uses the entire dataset for each iteration.

Accuracy

In general, stochastic Gradient Descent is Less accurate than gradient descent, as it calculates the gradient on single examples, which may not accurately represent the overall dataset. Gradient Descent is more accurate because it uses the average gradient calculated over the entire dataset.

There are other disadvantages to using SGD. The main drawback is that its convergence behaviour can be more erratic due to the random sampling of individual training examples. This can lead to less accurate results, as the algorithm may not converge to the true minimum of the cost function. Additionally, the learning rate, which determines the step size of each update to the model’s parameters, must be carefully chosen to ensure convergence.

It is however the method of choice in deep learning algorithms where SGD is often used in combination with other optimization techniques, such as momentum or adaptive learning rates

Stochastic Gradient Descent (SGD)

In stochastic gradient descent, the extreme case is the case where we have only one batch, that is we include the whole data set.

This process is called Stochastic Gradient Descent (SGD) (or also sometimes on-line gradient descent). This is relatively less common to see because in practice due to vectorized code optimizations it can be computationally much more efficient to evaluate the gradient for 100 examples, than the gradient for one example 100 times. Even though SGD technically refers to using a single example at a time to evaluate the gradient, you will hear people use the term SGD even when referring to mini-batch gradient descent (i.e. mentions of MGD for “Minibatch Gradient Descent”, or BGD for “Batch gradient descent” are rare to see), where it is usually assumed that mini-batches are used. The size of the mini-batch is a hyperparameter but it is not very common to cross-validate or bootstrap it. It is usually based on memory constraints (if any), or set to some value, e.g. 32, 64 or 128. We use powers of 2 in practice because many vectorized operation implementations work faster when their inputs are sized in powers of 2.

In our notes with SGD we mean stochastic gradient descent with mini-batches.

Stochastic Gradient Descent

Stochastic gradient descent (SGD) and variants thereof address some of the shortcomings of the Gradient descent method discussed above.

The underlying idea of SGD comes from the observation that the cost function, which we want to minimize, can almost always be written as a sum over n data points \{\mathbf{x}_i\}_{i=1}^n,

C(\mathbf{\theta}) = \sum_{i=1}^n c_i(\mathbf{x}_i,
\mathbf{\theta}).

Computation of gradients

This in turn means that the gradient can be computed as a sum over $i$-gradients

\nabla_\theta C(\mathbf{\theta}) = \sum_i^n \nabla_\theta c_i(\mathbf{x}_i,
\mathbf{\theta}).

Stochasticity/randomness is introduced by only taking the gradient on a subset of the data called minibatches. If there are n data points and the size of each minibatch is M, there will be n/M minibatches. We denote these minibatches by B_k where k=1,\cdots,n/M.

SGD example

As an example, suppose we have 10 data points (\mathbf{x}_1,\cdots, \mathbf{x}_{10}) and we choose to have M=5 minibathces, then each minibatch contains two data points. In particular we have $B_1 = (\mathbf{x}_1,\mathbf{x}_2), \cdots, B_5 = (\mathbf{x}9,\mathbf{x}{10})$. Note that if you choose M=1 you have only a single batch with all data points and on the other extreme, you may choose M=n resulting in a minibatch for each datapoint, i.e B_k = \mathbf{x}_k.

The idea is now to approximate the gradient by replacing the sum over all data points with a sum over the data points in one the minibatches picked at random in each gradient descent step

\nabla_{\theta}
C(\mathbf{\theta}) = \sum_{i=1}^n \nabla_\theta c_i(\mathbf{x}_i,
\mathbf{\theta}) \rightarrow \sum_{i \in B_k}^n \nabla_\theta
c_i(\mathbf{x}_i, \mathbf{\theta}).

The gradient step

Thus a gradient descent step now looks like

\theta_{j+1} = \theta_j - \eta_j \sum_{i \in B_k}^n \nabla_\theta c_i(\mathbf{x}_i,
\mathbf{\theta})

where k is picked at random with equal probability from [1,n/M]. An iteration over the number of minibathces (n/M) is commonly referred to as an epoch. Thus it is typical to choose a number of epochs and for each epoch iterate over the number of minibatches, as exemplified in the code below.

Simple example code

In [6]:

import numpy as np 

n = 100 #100 datapoints 
M = 5   #size of each minibatch
m = int(n/M) #number of minibatches
n_epochs = 10 #number of epochs

j = 0
for epoch in range(1,n_epochs+1):
    for i in range(m):
        k = np.random.randint(m) #Pick the k-th minibatch at random
        #Compute the gradient using the data in minibatch Bk
        #Compute new suggestion for 
        j += 1

Taking the gradient only on a subset of the data has two important benefits. First, it introduces randomness which decreases the chance that our opmization scheme gets stuck in a local minima. Second, if the size of the minibatches are small relative to the number of datapoints (M < n), the computation of the gradient is much cheaper since we sum over the datapoints in the k-th minibatch and not all n datapoints.

When do we stop?

A natural question is when do we stop the search for a new minimum? One possibility is to compute the full gradient after a given number of epochs and check if the norm of the gradient is smaller than some threshold and stop if true. However, the condition that the gradient is zero is valid also for local minima, so this would only tell us that we are close to a local/global minimum. However, we could also evaluate the cost function at this point, store the result and continue the search. If the test kicks in at a later stage we can compare the values of the cost function and keep the \theta that gave the lowest value.

Slightly different approach

Another approach is to let the step length \eta_j depend on the number of epochs in such a way that it becomes very small after a reasonable time such that we do not move at all. Such approaches are also called scaling. There are many such ways to scale the learning rate and discussions here. See also https://towardsdatascience.com/learning-rate-schedules-and-adaptive-learning-rate-methods-for-deep-learning-2c8f433990d1 for a discussion of different scaling functions for the learning rate.

Time decay rate

As an example, let e = 0,1,2,3,\cdots denote the current epoch and let t_0, t_1 > 0 be two fixed numbers. Furthermore, let t = e \cdot m + i where m is the number of minibatches and i=0,\cdots,m-1. Then the function \eta_j(t; t_0, t_1) = \frac{t_0}{t+t_1} goes to zero as the number of epochs gets large. I.e. we start with a step length \eta_j (0; t_0, t_1) = t_0/t_1 which decays in time t.

In this way we can fix the number of epochs, compute \theta and evaluate the cost function at the end. Repeating the computation will give a different result since the scheme is random by design. Then we pick the final \theta that gives the lowest value of the cost function.

In [7]:

import numpy as np 

def step_length(t,t0,t1):
    return t0/(t+t1)

n = 100 #100 datapoints 
M = 5   #size of each minibatch
m = int(n/M) #number of minibatches
n_epochs = 500 #number of epochs
t0 = 1.0
t1 = 10

eta_j = t0/t1
j = 0
for epoch in range(1,n_epochs+1):
    for i in range(m):
        k = np.random.randint(m) #Pick the k-th minibatch at random
        #Compute the gradient using the data in minibatch Bk
        #Compute new suggestion for theta
        t = epoch*m+i
        eta_j = step_length(t,t0,t1)
        j += 1

print("eta_j after %d epochs: %g" % (n_epochs,eta_j))

Code with a Number of Minibatches which varies

In the code here we vary the number of mini-batches.

In [8]:

# Importing various packages
from math import exp, sqrt
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt

n = 100
x = 2*np.random.rand(n,1)
y = 4+3*x+np.random.randn(n,1)

X = np.c_[np.ones((n,1)), x]
XT_X = X.T @ X
theta_linreg = np.linalg.inv(X.T @ X) @ (X.T @ y)
print("Own inversion")
print(theta_linreg)
# Hessian matrix
H = (2.0/n)* XT_X
EigValues, EigVectors = np.linalg.eig(H)
print(f"Eigenvalues of Hessian Matrix:{EigValues}")

theta = np.random.randn(2,1)
eta = 1.0/np.max(EigValues)
Niterations = 1000


for iter in range(Niterations):
    gradients = 2.0/n*X.T @ ((X @ theta)-y)
    theta -= eta*gradients
print("theta from own gd")
print(theta)

xnew = np.array([[0],[2]])
Xnew = np.c_[np.ones((2,1)), xnew]
ypredict = Xnew.dot(theta)
ypredict2 = Xnew.dot(theta_linreg)

n_epochs = 50
M = 5   #size of each minibatch
m = int(n/M) #number of minibatches
t0, t1 = 5, 50

def learning_schedule(t):
    return t0/(t+t1)

theta = np.random.randn(2,1)

for epoch in range(n_epochs):
# Can you figure out a better way of setting up the contributions to each batch?
    for i in range(m):
        random_index = M*np.random.randint(m)
        xi = X[random_index:random_index+M]
        yi = y[random_index:random_index+M]
        gradients = (2.0/M)* xi.T @ ((xi @ theta)-yi)
        eta = learning_schedule(epoch*m+i)
        theta = theta - eta*gradients
print("theta from own sdg")
print(theta)

plt.plot(xnew, ypredict, "r-")
plt.plot(xnew, ypredict2, "b-")
plt.plot(x, y ,'ro')
plt.axis([0,2.0,0, 15.0])
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'Random numbers ')
plt.show()

Replace or not

In the above code, we have use replacement in setting up the mini-batches. The discussion here may be useful.

Second moment of the gradient

In stochastic gradient descent, with and without momentum, we still have to specify a schedule for tuning the learning rates \eta_t as a function of time. As discussed in the context of Newton's method, this presents a number of dilemmas. The learning rate is limited by the steepest direction which can change depending on the current position in the landscape. To circumvent this problem, ideally our algorithm would keep track of curvature and take large steps in shallow, flat directions and small steps in steep, narrow directions. Second-order methods accomplish this by calculating or approximating the Hessian and normalizing the learning rate by the curvature. However, this is very computationally expensive for extremely large models. Ideally, we would like to be able to adaptively change the step size to match the landscape without paying the steep computational price of calculating or approximating Hessians.

During the last decade a number of methods have been introduced that accomplish this by tracking not only the gradient, but also the second moment of the gradient. These methods include AdaGrad, AdaDelta, Root Mean Squared Propagation (RMS-Prop), and ADAM.

Challenge: Choosing a Fixed Learning Rate

A fixed \eta is hard to get right:

1. If \eta is too large, the updates can overshoot the minimum, causing oscillations or divergence

2. If \eta is too small, convergence is very slow (many iterations to make progress)

In practice, one often uses trial-and-error or schedules (decaying \eta over time) to find a workable balance. For a function with steep directions and flat directions, a single global \eta may be inappropriate:

1. Steep coordinates require a smaller step size to avoid oscillation.

2. Flat/shallow coordinates could use a larger step to speed up progress.

3. This issue is pronounced in high-dimensional problems with sparse or varying-scale features – we need a method to adjust step sizesper feature.

Motivation for Adaptive Step Sizes

1. Instead of a fixed global \eta, use an adaptive learning rate for each parameter that depends on the history of gradients.

2. Parameters that have large accumulated gradient magnitude should get smaller steps (they've been changing a lot), whereas parameters with small or infrequent gradients can have larger relative steps.

3. This is especially useful for sparse features: Rarely active features accumulate little gradient, so their learning rate remains comparatively high, ensuring they are not neglected

4. Conversely, frequently active features accumulate large gradient sums, and their learning rate automatically decreases, preventing too-large updates

5. Several algorithms implement this idea (AdaGrad, RMSProp, AdaDelta, Adam, etc.). We will derive AdaGrad, one of the first adaptive methods.

AdaGrad algorithm, taken from Goodfellow et al

Figure 1:

Derivation of the AdaGrad Algorithm

Accumulating Gradient History.

1. AdaGrad maintains a running sum of squared gradients for each parameter (coordinate)

2. Let g_t = \nabla C_{i_t}(x_t) be the gradient at step t (or a subgradient for nondifferentiable cases).

3. Initialize r_0 = 0 (an all-zero vector in \mathbb{R}^d).

4. At each iteration t, update the accumulation:

r_t = r_{t-1} + g_t \circ g_t,

1. Here g_t \circ g_t denotes element-wise square of the gradient vector. g_t^{(j)} = g_{t-1}^{(j)} + (g_{t,j})^2 for each parameter j.

2. We can view H_t = \mathrm{diag}(r_t) as a diagonal matrix of past squared gradients. Initially H_0 = 0.

AdaGrad Update Rule Derivation

We scale the gradient by the inverse square root of the accumulated matrix H_t. The AdaGrad update at step t is:

\theta_{t+1} =\theta_t - \eta H_t^{-1/2} g_t,

where H_t^{-1/2} is the diagonal matrix with entries (r_{t}^{(1)})^{-1/2}, \dots, (r_{t}^{(d)})^{-1/2} In coordinates, this means each parameter j has an individual step size:

\theta_{t+1,j} =\theta_{t,j} -\frac{\eta}{\sqrt{r_{t,j}}}g_{t,j}.

In practice we add a small constant \epsilon in the denominator for numerical stability to avoid division by zero:

\theta_{t+1,j}= \theta_{t,j}-\frac{\eta}{\sqrt{\epsilon + r_{t,j}}}g_{t,j}.

Equivalently, the effective learning rate for parameter j at time t is \displaystyle \alpha_{t,j} = \frac{\eta}{\sqrt{\epsilon + r_{t,j}}}. This decreases over time as r_{t,j} grows.

AdaGrad Properties

1. AdaGrad automatically tunes the step size for each parameter. Parameters with more volatile or large gradients get smaller steps, and those with small or infrequent gradients get relatively larger steps

2. No manual schedule needed: The accumulation r_t keeps increasing (or stays the same if gradient is zero), so step sizes \eta/\sqrt{r_t} are non-increasing. This has a similar effect to a learning rate schedule, but individualized per coordinate.

3. Sparse data benefit: For very sparse features, r_{t,j} grows slowly, so that feature’s parameter retains a higher learning rate for longer, allowing it to make significant updates when it does get a gradient signal

4. Convergence: In convex optimization, AdaGrad can be shown to achieve a sub-linear convergence rate comparable to the best fixed learning rate tuned for the problem

It effectively reduces the need to tune \eta by hand.

1. Limitations: Because r_t accumulates without bound, AdaGrad’s learning rates can become extremely small over long training, potentially slowing progress. (Later variants like RMSProp, AdaDelta, Adam address this by modifying the accumulation rule.)

RMSProp: Adaptive Learning Rates

Addresses AdaGrad’s diminishing learning rate issue. Uses a decaying average of squared gradients (instead of a cumulative sum):

v_t = \rho v_{t-1} + (1-\rho)(\nabla C(\theta_t))^2,

with \rho typically 0.9 (or 0.99).

1. Update: \theta_{t+1} = \theta_t - \frac{\eta}{\sqrt{v_t + \epsilon}} \nabla C(\theta_t).

2. Recent gradients have more weight, so v_t adapts to the current landscape.

3. Avoids AdaGrad’s “infinite memory” problem – learning rate does not continuously decay to zero.

RMSProp was first proposed in lecture notes by Geoff Hinton, 2012 - unpublished.)

RMSProp algorithm, taken from Goodfellow et al

Figure 1:

Adam Optimizer

Why combine Momentum and RMSProp? Motivation for Adam: Adaptive Moment Estimation (Adam) was introduced by Kingma an Ba (2014) to combine the benefits of momentum and RMSProp.

1. Fast convergence by smoothing gradients (accelerates in long-term gradient direction).

2. Adaptive rates (RMSProp): Per-dimension learning rate scaling for stability (handles different feature scales, sparse gradients).

3. Adam uses both: maintains moving averages of both first moment (gradients) and second moment (squared gradients)

4. Additionally, includes a mechanism to correct the bias in these moving averages (crucial in early iterations)

Result: Adam is robust, achieves faster convergence with less tuning, and often outperforms SGD (with momentum) in practice.

ADAM optimizer

In ADAM, we keep a running average of both the first and second moment of the gradient and use this information to adaptively change the learning rate for different parameters. The method is efficient when working with large problems involving lots data and/or parameters. It is a combination of the gradient descent with momentum algorithm and the RMSprop algorithm discussed above.

Why Combine Momentum and RMSProp?

1. Momentum: Fast convergence by smoothing gradients (accelerates in long-term gradient direction).

2. Adaptive rates (RMSProp): Per-dimension learning rate scaling for stability (handles different feature scales, sparse gradients).

3. Adam uses both: maintains moving averages of both first moment (gradients) and second moment (squared gradients)

4. Additionally, includes a mechanism to correct the bias in these moving averages (crucial in early iterations)

Result: Adam is robust, achieves faster convergence with less tuning, and often outperforms SGD (with momentum) in practice

Adam: Exponential Moving Averages (Moments)

Adam maintains two moving averages at each time step t for each parameter w: First moment (mean) m_t.

The Momentum term

m_t = \beta_1m_{t-1} + (1-\beta_1)\, \nabla C(\theta_t),

Second moment (uncentered variance) v_t.

The RMS term

v_t = \beta_2v_{t-1} + (1-\beta_2)(\nabla C(\theta_t))^2,

with typical \beta_1 = 0.9, \beta_2 = 0.999. Initialize m_0 = 0, v_0 = 0.

These are biased estimators of the true first and second moment of the gradients, especially at the start (since m_0,v_0 are zero)

Adam: Bias Correction

To counteract initialization bias in m_t, v_t, Adam computes bias-corrected estimates

\hat{m}_t = \frac{m_t}{1 - \beta_1^t}, \qquad \hat{v}_t = \frac{v_t}{1 - \beta_2^t}.

• When t is small, 1-\beta_i^t \approx 0, so \hat{m}_t, \hat{v}_t significantly larger than raw m_t, v_t, compensating for the initial zero bias.

• As t increases, 1-\beta_i^t \to 1, and \hat{m}_t, \hat{v}_t converge to m_t, v_t.

• Bias correction is important for Adam’s stability in early iterations

Warning:

Output truncated. This notebook contains too many cells to display efficiently.