Gradient descent doesn't find solution to ordinary least squares on this dataset?

Question

I have been studying linear regression and tried it on below set {(x,y)}, where x specified the area of house in square-feet, and y specified the price in dollars. This is the first example in Andrew Ng Notes.

2104,400
1600,330
2400,369
1416,232
3000,540

I developed a sample code but when I run it, the cost is increasing with each step whereas it should be decreasing with each step. Code and output given below. bias is W₀X₀, where X₀=1. featureWeights is an array of [X₁,X₂,...,X_N]

I also tried an online python solution available here, and explained here. But this example is also giving the same output.

Where is the gap in understanding the concept?

Code:

package com.practice.cnn;

import java.util.Arrays;

public class LinearRegressionExample {

    private float ALPHA = 0.0001f;
    private int featureCount = 0;
    private int rowCount = 0;

    private float bias = 1.0f;
    private float[] featureWeights = null;

    private float optimumCost = Float.MAX_VALUE;

    private boolean status = true;

    private float trainingInput[][] = null;
    private float trainingOutput[] = null;

    public void train(float[][] input, float[] output) {
        if (input == null || output == null) {
            return;
        }

        if (input.length != output.length) {
            return;
        }

        if (input.length == 0) {
            return;
        }

        rowCount = input.length;
        featureCount = input[0].length;

        for (int i = 1; i < rowCount; i++) {
            if (input[i] == null) {
                return;
            }

            if (featureCount != input[i].length) {
                return;
            }
        }

        featureWeights = new float[featureCount];
        Arrays.fill(featureWeights, 1.0f);

        bias = 0;   //temp-update-1
        featureWeights[0] = 0;  //temp-update-1

        this.trainingInput = input;
        this.trainingOutput = output;

        int count = 0;
        while (true) {
            float cost = getCost();

            System.out.print("Iteration[" + (count++) + "] ==> ");
            System.out.print("bias -> " + bias);
            for (int i = 0; i < featureCount; i++) {
                System.out.print(", featureWeights[" + i + "] -> " + featureWeights[i]);
            }
            System.out.print(", cost -> " + cost);
            System.out.println();

//          if (cost > optimumCost) {
//              status = false;
//              break;
//          } else {
//              optimumCost = cost;
//          }

            optimumCost = cost;

            float newBias = bias + (ALPHA * getGradientDescent(-1));

            float[] newFeaturesWeights = new float[featureCount];
            for (int i = 0; i < featureCount; i++) {
                newFeaturesWeights[i] = featureWeights[i] + (ALPHA * getGradientDescent(i));
            }

            bias = newBias;

            for (int i = 0; i < featureCount; i++) {
                featureWeights[i] = newFeaturesWeights[i];
            }
        }
    }

    private float getCost() {
        float sum = 0;
        for (int i = 0; i < rowCount; i++) {
            float temp = bias;
            for (int j = 0; j < featureCount; j++) {
                temp += featureWeights[j] * trainingInput[i][j];
            }

            float x = (temp - trainingOutput[i]) * (temp - trainingOutput[i]);
            sum += x;
        }
        return (sum / rowCount);
    }

    private float getGradientDescent(final int index) {
        float sum = 0;
        for (int i = 0; i < rowCount; i++) {
            float temp = bias;
            for (int j = 0; j < featureCount; j++) {
                temp += featureWeights[j] * trainingInput[i][j];
            }

            float x = trainingOutput[i] - (temp);
            sum += (index == -1) ? x : (x * trainingInput[i][index]);
        }
        return ((sum * 2) / rowCount);
    }

    public static void main(String[] args) {
        float[][] input = new float[][] { { 2104 }, { 1600 }, { 2400 }, { 1416 }, { 3000 } };

        float[] output = new float[] { 400, 330, 369, 232, 540 };

        LinearRegressionExample example = new LinearRegressionExample();
        example.train(input, output);
    }
}

Output:

Iteration[0] ==> bias -> 0.0, featureWeights[0] -> 0.0, cost -> 150097.0
Iteration[1] ==> bias -> 0.07484, featureWeights[0] -> 168.14847, cost -> 1.34029099E11
Iteration[2] ==> bias -> -70.60721, featureWeights[0] -> -159417.34, cost -> 1.20725801E17
Iteration[3] ==> bias -> 67012.305, featureWeights[0] -> 1.51299168E8, cost -> 1.0874295E23
Iteration[4] ==> bias -> -6.3599688E7, featureWeights[0] -> -1.43594258E11, cost -> 9.794949E28
Iteration[5] ==> bias -> 6.036088E10, featureWeights[0] -> 1.36281745E14, cost -> 8.822738E34
Iteration[6] ==> bias -> -5.7287012E13, featureWeights[0] -> -1.29341617E17, cost -> Infinity
Iteration[7] ==> bias -> 5.4369677E16, featureWeights[0] -> 1.2275491E20, cost -> Infinity
Iteration[8] ==> bias -> -5.1600908E19, featureWeights[0] -> -1.1650362E23, cost -> Infinity
Iteration[9] ==> bias -> 4.897313E22, featureWeights[0] -> 1.1057068E26, cost -> Infinity
Iteration[10] ==> bias -> -4.6479177E25, featureWeights[0] -> -1.0493987E29, cost -> Infinity
Iteration[11] ==> bias -> 4.411223E28, featureWeights[0] -> 9.959581E31, cost -> Infinity
Iteration[12] ==> bias -> -4.186581E31, featureWeights[0] -> -Infinity, cost -> Infinity
Iteration[13] ==> bias -> Infinity, featureWeights[0] -> NaN, cost -> NaN
Iteration[14] ==> bias -> NaN, featureWeights[0] -> NaN, cost -> NaN

Answer 1

The short answer is that your step size is too big. Instead of descending the canyon wall, your step is so big that you're jumping across from one side to higher up on the other!

Cost function below:

The long answer is that it's difficult for a naive gradient descent to solve this problem because the level sets of your cost function are highly elongated ellipses rather than circles. To robustly solve this problem, note that there are more sophisticated ways to choose:

• a step size (than hardcoding a constant).
• a step direction (than gradient descent).

Underlying problem

The underlying problem is that level sets of your cost function are highly elongated ellipses, and this causes problems for gradient descent. The below figure shows level sets for the cost function.

• With highly elliptical level sets, the direction of steepest descent may barely align with the direction of the solution. For example in this problem, the intercept term (what you call "bias") needs to travel a great distance (from $0$ to $\approx 26.789$ along the canyon floor) but it is for the other feature where the partial derivative has a much larger slope.
• If step size is too big, you will literally jump over the lower blue region and ascend instead of descend.
• BUT if you if you reduce your step size, your progress in getting $\theta_0$ to the proper value becomes painfully slow.

I suggest reading this answer on Quora.

Quick fix 1:

Change your code to private float ALPHA = 0.0000002f; and you'll stop overshooting.

Quick fix 2:

If you rescale your X data to 2.104, 1.600, etc... your level sets become spherical and gradient descent quickly converges with a higher learning rate. This lowers the condition number of your design matrix $X'X$.

More advanced fixes

If the goal were to efficiently solve ordinary least squares rather than simply learn gradient descent for a class, observe that:

• There are more sophisticated ways to calculate step size, such as line search and the Armijo rule.
• Near an answer where local conditions prevail, Newton's method obtains quadratic convergence and is a great way to choose a step direction and size.
• Solving least squares is equivalent to solving a linear system. Modern algorithms don't use naive gradient descent. Instead:
  • For small systems ($k$ on the order of several thousand or less), they use something like QR decomposition with partial pivoting.
  • For large systems, they do formulate it is an optimization problem and use iterative methods such as the Krylov subspace methods.

Note that there are many packages which will solve the linear system $(X'X) b = X'y$ for $b$ and you can check the results of your gradient descent algorithm against that.

The actual solution is

  26.789880528523071
   0.165118878075797

You will find that those achieve the minimum value for the cost function.

Answer 2

As Matthew (Gunn) has already indicated, the contours of the 3-dimensional cost or performance function are highly elliptical in this case. Since your Java code uses a single step-size value for the gradient descent calculations, the updates to the weights (i.e. the y-axis intercept and slope of the linear function) are both governed by this single step-size.

As a result, the very small step-size that is required to control the update of the weight associated with the larger gradient (in this case, the slope of the linear function) drastically limits how fast the other weight with the smaller gradient (the y-axis intercept of the linear function) is updated. Under the current conditions, the latter weight does not converge to its true value of approximately 26.7.

Given the time and effort that you have invested in writing your Java code, I would suggest modifying it to use two discrete step-size values, an appropriate step-size for each weight. Andrew Ng suggests in his notes that it is better to use feature scaling to ensure that the contours of the cost function are more regular (i.e. circular) in form. However, modifying your Java code to use a different step-size for each weight might be a good exercise in addition to looking at feature scaling.

Another idea to consider is how the initial weight values are picked. In your Java code you initialized both values to zero. It is also quite common to initialize the weights to small, fractional values. In this particular case, however, both of these approaches would not work in light of the highly elliptical (i.e. non-circular) contours of the three-dimensional cost function. Given the weights for this problem can be found using other methods, like the solution for the linear system suggested by Matthew at the end of his post, you could try to initialize the weights to values closer to the correct weights and see how your original code using a single step-size converges.

The Python code you found approaches the solution in the same way as your Java code - both use a single step-size parameter. I modified this Python code to use a different step-size for each weight. I have included it below.

from numpy import *

def compute_error_for_line_given_points(b, m, points):
    totalError = 0
    for i in range(0, len(points)):
        x = points[i, 0]
        y = points[i, 1]
        totalError += (y - (m * x + b)) ** 2
    return totalError / float(len(points))

def step_gradient(b_current, m_current, points, learningRate_1, learningRate_2):
    b_gradient = 0
    m_gradient = 0
    N = float(len(points))
    for i in range(0, len(points)):
        x = points[i, 0]
        y = points[i, 1]
        b_gradient += -(2/N) * (y - ((m_current * x) + b_current))
        m_gradient += -(2/N) * x * (y - ((m_current * x) + b_current))
    new_b = b_current - (learningRate_1 * b_gradient)
    new_m = m_current - (learningRate_2 * m_gradient)
    return [new_b, new_m]

def gradient_descent_runner(points, starting_b, starting_m, learning_rate_1, learning_rate_2, num_iterations):
    b = starting_b
    m = starting_m
    for i in range(num_iterations):
        b, m = step_gradient(b, m, array(points), learning_rate_1, learning_rate_2)
    return [b, m]

def run():
    #points = genfromtxt("data.csv", delimiter=",")
    #learning_rate = 0.0001
    #num_iterations = 200

    points = genfromtxt("test_set.csv", delimiter=",")
    learning_rate_1 = 0.5
    learning_rate_2 = 0.0000001
    num_iterations = 1000

    initial_b = 0 # initial y-intercept guess
    initial_m = 0 # initial slope guess


    print("Starting gradient descent at b = {0}, m = {1}, error = {2}".format(initial_b, initial_m, compute_error_for_line_given_points(initial_b, initial_m, points)))
    print("Running...")

    [b, m] = gradient_descent_runner(points, initial_b, initial_m, learning_rate_1, learning_rate_2, num_iterations)

    print("After {0} iterations b = {1}, m = {2}, error = {3}".format(num_iterations, b, m, compute_error_for_line_given_points(b, m, points)))

if __name__ == '__main__':
    run()

It runs under Python 3, which requires the parentheses around the argument for the "print" statements. Otherwise it will run under Python 2 by removing the parentheses. You'll need to create a CSV-file with the data from Andrew Ng's example.

Use can cross-reference the Python code to check your Java code.