Gradient descent in Octave

Gradient descent is one of the simplest method to fit a model of a given form from a bunch of data. For example, you may want to know which is the best (in terms of mean squared error) line fitting a series of points. In that case, you're trying to fit:

y = m * x + q

where y and x are known vectors, collected from some source:

y = [1 2 3.1 4 5.5]
x = [0 1 2 3 4]

and m and q are the parameters you want to discover.

In reality, the function may be the predicted rating of a book or movie given the preferences of a user, or the number of views of a page as a function of length and tags.

The (high-school) mathematics behind it

Given a vector of parameters, in the one-dimension linear case [m, q], one way to quantify the error on the training set is to calculate the (mean) squared error:

SE = sum{[y(i) - m * x(i) - q]^2}

where i varies from 1 to the number of samples N. The mean error should be divided by N, but for the purpose of our discussion nothing changes as minimizing this error or its quotient by N is the same.

When you wanted to minimize a function in high-school, what you do is to take the derivatives with respect to the variables you have to find:

dSE/dm = sum{2*[y(i) - m * x(i) - q]*(-x(i))}
dSE/dq = sum{2*[y(i) - m * x(i) - q]*(-1)}

These expressions are usually too difficult to minimize globally in a closed form (coming up with an algebric solution), but we can use iterative methods to get close to a point where the gradient G = [dSE/dm dSE/dq] is close as possible to [0 0], which means the SE is at a local minimum. There's no way to be sure it is the unique minimum or that it is the global one without exploring the whole domain.

Our iterative solution, gradient descent, is to:

1. pick starting points at random for m and q.
2. Move a bit into the opposite direction of the gradient G (which is the fastest direction to decrease SE).
3. If we have not reached a fixed number of iterations or an acceptable value for SE, we restart from the new values of m and q.

We moderate this movement by multiplying the gradient for a small number called learning rate, essentially to avoid big jumps that miss a minimum.

A complete example

With this Octave code, we try to fit a parabola of the form y = a*x^2 + c. We could do the same with y = a*x^2 + b*x + c, and hopefully we would find b would be close to 0 for these data.

This is the training set (it is usually bigger):

% training set: samples from the real world
x = [-4, -3, -2, -1, 0, 1, 2, 3, 4];
y = [26.1, 19, 14.1, 10.9, 10.1, 11.1, 13.9, 19.1, 25.9];
% our true model would produce: 26, 19, 14, 11, 10, 11, 14, 19, 26
% it is in fact y = x^2 + 10

It is handy to be able to calculate the MSE:

% looking for a model y = f(x) = a*x^2 + c
% MSE = sum{[y-f(x)]^2}
function result = mse(y, x, a, c)
    result = 0;
    for i=1:size(x, 2)  % could be done with matrices, but I prefer to be clear
        result = result + (y(i) - (a * x(i)^2 + c))^2;
    end
    result = result / size(x, 2);
end

Since the number of samples is constant, we can minimize the SE instead of the MSE.

These are the derivatives with respect to the parameters to be found, a and c:

% MSE derivative with respect to a: sum{2[y-f(x)]*(-2x)}
function d = a_direction_component(y, x, a, c)
    d = sum(2*(y - a * (x .^ 2) - c) .* (-2 * x .^ 2));
end

% MSE derivative with respect to c: sum{2[y-f(x)]*(-1)}
function d = c_direction_component(y, x, a, c)
    d = sum(2*(y - a * x .^ 2 - c) .* (-1));
end

We start from [0, 0], but we there is no preferred point and we could pick a 2D point at random. The learning rate step is a small number:

% initialization
a = 0;
c = 0;
step = 0.0001;
mse(y, x, a, c)

We perform some thosands iteration, as gradient descent is one of the slowest methods to converge:

for i=1:10000
    sprintf("Iteration %d", i)
    gradient = [a_direction_component(y, x, a, c) c_direction_component(y, x, a, c)]
    a = a - step * gradient(1);
    c = c - step * gradient(2);
    sprintf("Model learned: y = %f.2 * x^2 + %f.2", a, c)
    current_mse = mse(y, x, a, c)
end

The execution is the following:

ans = Iteration 1
gradient =

  -5235.60   -300.40

ans = Model learned: y = 0.523560 * x^2 + 0.030040
current_mse =  181.14
ans = Iteration 2
gradient =

  -3745.67   -237.03

ans = Model learned: y = 0.898127 * x^2 + 0.053743
current_mse =  113.73
ans = Iteration 3
gradient =

  -2679.21   -191.66

ans = Model learned: y = 1.166047 * x^2 + 0.072909
current_mse =  79.155
...
ans = Iteration 10000
gradient =

   0.0026914  -0.0316712

ans = Model learned: y = 0.998938 * x^2 + 10.027546
current_mse =  0.0083345

Conclusions

Gradient descent is a very simple method to implement, and can in theory learn any model with any number of parameters (linear, quadratic, exponential...). However, it is prone to slow convergence and to fall into local minima of the MSE (think of holes in a garden) without being able to escape, since it's like a ball rolling down a slope.

In any case, the basic flow of iterative algorithms that try to improve a metric on a training set is very popular in machine learning and gradient descent is the first example of them.