Algorithm explained: Linear regression using gradient descent with PHP
Tags: machinelearning, php, computerscience, tutorial
Published: 2020-12-11 11:02:47

> Part 4 of **Algorithms explained**! Every few weeks I write about an algorithm and explain and implement it!
**Is there an algorithm you always wanted to know about? Leave a comment!**

I once took an online class on machine learning by the amazing Andrew Ng from Stanford for my bachelors studies. One of the first ML algorithms I ever had the chance to look at, was linear regression.

This algorithm really struck me. It allows machines to do what humans are incredibly good at: Pattern recognition. It allows machines to create a linear function that _more or less_ matches seemingly unrelated data by applying mathematics to it!

**In this post I'll explain how linear regression works and implement it in PHP!**

# So what does linear regression actually _do_?

In one sentence: It creates a linear function that represents the trend of something. The go-to example is usually house prices compared to their size. The larger a house, the more expensive it is. Sounds reasonable, right? You would put a known house size, for example in square meters, into the function and you would get an estimated price out of it. 

Let's have a look at this data set:

![House prices and sizes compared to each other in a two dimensional graph. There's some correlation visible.](https://dev-to-uploads.s3.amazonaws.com/i/vfrmsopqohvpkfa4dyx3.png)

On first glance, the data seems to correlate. The larger the house, the more expensive it is. The data is a bit scattered, but there's an overall trend visible:

![The same graph as above, now with an added trend line](https://dev-to-uploads.s3.amazonaws.com/i/ps4ub4lebew8ekatcbl4.png)

_(This is an approximation)_

Of course there's usually more factors to the price of a house in real life: Is there a grocery store nearby, is it modern, was it once considered haunted, what's the size of its backyard? These factors make up the variation in the data - it's not linear.

**The algorithm itself works like this:** 
 1. Define any odd linear function as a trend line (usually random)
 2. Measure how far off it is by calculating the average distance between predicted Y and actual Y of every data point (the so called "error")
 3. Adjust the trend line based on the measurement (the "gradient descent")
 4. Repeat Step 2 and 3 until the average distance has reached a minimum (not necessarily 0)

Let me illustrate with the above example:

![The different steps illustrated](https://dev-to-uploads.s3.amazonaws.com/i/rz8y85oqi457al4lnz6q.png)

_(Again, this is an approximation)_

So, to implement the algorithm, we need three different pieces:

 * The definition of a linear function
 * The measurement
 * The adjustment
 * Repeating until it's at a minimum

Let's start with...

# The definition of a linear function

A linear function {% katex inline %}f(x){% endkatex %} usually comes in the form of:

{% katex %}
f(x) = ax + b
{% endkatex %}

Whereas {%katex inline%}a{%endkatex%} basically defines the *slope* of the function and {%katex inline%}b{%endkatex%} moves it along the Y axis. 

We can make a more general definition of a linear function by introducing a variable called {%katex inline%}x_0{%endkatex%} into the function and renaming a and b to {%katex inline%}c_0{%endkatex%}  and {%katex inline%}c_0{%endkatex%} ("c" stands for "coefficient"):

{%katex%}
f(x_0, x_1) = c_0x_0 + c_1x_1
{%endkatex%}

In PHP, this is pretty straight forward:

```php
<?php

declare(strict_types=1);

function linearFunction (float $c0, float $x0, float $c1, float $x1) : float {
    return $c0 * $x0 + $c1 * $x1;
}
```

We can randomly define `$c0` and `$c1` to get a first starting point:

```php
$c0 = mt_rand(-100, 100); // -100 and 100 are almost parallel to the X and Y axis
$c1 = mt_rand(-100, 100);
```

So far, so good. Let's generate some data next!

# Getting data

I'm going to generate the data randomly with a correlation. To get the linear function {%katex inline%}f(x_0, x_1) = c_0x_0 + c_1x_1{%endkatex%} back into it's original form {%katex inline%}f(x) = c_0 + c_1x{%endkatex%}, we set {%katex inline%}x_0{%endkatex%} to always be 1:

```php
$data = [
    [1, mt_rand(0, 100), mt_rand(0, 100)],
    [1, mt_rand(100, 200), mt_rand(100, 200)],
    [1, mt_rand(200, 300), mt_rand(200, 300)],
    [1, mt_rand(300, 400), mt_rand(300, 400)],
    [1, mt_rand(400, 500), mt_rand(400, 500)],
    [1, mt_rand(500, 600), mt_rand(500, 600)],
    [1, mt_rand(600, 700), mt_rand(600, 700)],
    [1, mt_rand(700, 800), mt_rand(700, 800)],
    [1, mt_rand(800, 900), mt_rand(800, 900)],
    [1, mt_rand(900, 1000), mt_rand(900, 1000)],
];
```

For illustration purposes, I'll be going to work with this data set and predefined values for `$c0` and `$c1`:

```php
$data = [
  [1, 48, 36],
  [1, 121, 104],
  [1, 285, 208],
  [1, 345, 373],
  [1, 402, 442],
  [1, 576, 582],
  [1, 617, 629],
  [1, 723, 787],
  [1, 852, 857],
  [1, 965, 950],
];
$c0 = 16;
$c1 = 15;
```

To graph a few things, I'll use Desmos, which is an amazing tool to fiddle around with maths. On a graph, the data looks like this:

![The above data, graphed as a scatter plot](https://dev-to-uploads.s3.amazonaws.com/i/x3qi3yaw2x1fnbu9ak73.png)

We can see a clear trend there. The trend line will likely be going from bottom left to top right.

With the data at hand, we can start to...

# Measure the error

So, to get a good **estimate** of how far off the linear function is, we calculate the average squared distance of the function's output Y and the actual Y. Squared, because otherwise negative numbers might emerge, which is something we don't want.

For this, we'll define yet another function that looks like this:

{%katex%}
error(x_0^{(1..n)}, x_1^{(1..n)}, y^{(1..n)}) = \frac{1}{n}\sum_{i=1}^{n}{(y^{(i)} - f(x_0^{(i)}, x_1^{(i)}))^2}
{%endkatex%}

**(Notice: {%katex inline%}x^{(i)}{%endkatex%} is not "x to the power of (i)", but rather "the x at position (i) in the list)**

Or, in PHP:

```php
function squaredError(float $c0, float $c1, array $data): float {
  return array_sum(
    array_map(
      function ($point) use ($c0, $c1) {
        return ($point[2] - linearFunction($c0, $point[0], $c1, $point[1])) ** 2;
      },
      $data
    )
  ) / count($data);
}

var_dump(squaredError($c0, $c1, $data)); // 73679488.2
```

This formula sums up the squares (because otherwise there's negative numbers) of the differences between predicted Y of the trend line and actual Y for every data point. It then divides it by the total number of data points, resulting in an average.

Ok, but `73679488.2` is quite the number. Let's see what this would look like by adding the linear function to the graph:

![Desmos showing the data and the graph](https://dev-to-uploads.s3.amazonaws.com/i/6lba991qfjoywm5b1ee5.png)

Yup, that's pretty far off. The Y-value for each point is far off to the top, accounting for this very large error. But that was the random definition. When can start improving from here.

# Calculating the adjustment

To improve the trend line's fit to the data, we'll use "gradient descent". According to Wikipedia, gradient descent is "a first-order iterative optimization algorithm for finding a local minimum of a differentiable function." And that pretty much sums it up.

The "differentiable function" is, in the case of linear regression, the error function. So gradient descent basically tries to find the local minimum of the error function (thus maximizing the fit of the original linear function to the data). To achieve this, it can adjust two things.

Since X is always given (it's the data, we can't change that), gradient descent can adjust {%katex inline%}c_{0}{%endkatex%} and {%katex inline%}c_{1}{%endkatex%}.

Gradient descent decides on a level of adjustment by looking at the error function's steepness. Again, Wikipedia has a good analogy on this: Imagine being stuck in the mountains. It's very foggy and you need to get to a valley to escape. All you can rely on is local information (stuff that is visible around you). Using gradient descent, you would go in the steepest direction downwards, eventually reaching a minimum - a valley. The larger the steep, the larger the step you would take.

Let me illustrate in 2D:

![Gradient descent in two dimensions. A parabola showing the relation between a single coefficient C and the error. Red dots indicate which C is chosen and in which direction it is altered.](https://dev-to-uploads.s3.amazonaws.com/i/jlxzh2pu49nhkevpohlh.png)

The red dots on the graph show the value for C for each iteration of gradient descent. The arrows indicate in which direction the algorithm would alter C in order to find a smaller error.

(Since we have two C's ({%katex inline%}c_{0}{%endkatex%} and {%katex inline%}c_{1}{%endkatex%}) and one extra dimension for the error, the graph would be 3D, though.)

To get to the _amount_ of change for each value of {%katex inline%}c_{0}{%endkatex%} and {%katex inline%}c_{1}{%endkatex%}, we can do partial derivatives with respect to {%katex inline%}c_{0}{%endkatex%} and {%katex inline%}c_{1}{%endkatex%} (i.e. what's the steepness of the function for one variable, while the other is held constant?). I won't go too deep into the calculus here. We would end up with this formula:

{%katex%}
D(c_{m}) = \frac{-2}{n} \sum_{i=1}^{n}{((y^{(i)} - f(x_0^{(i)}, x_1^{(i)})) \times x_m^{(i)})}
{%endkatex%}

{%katex inline%}D(c_{m}){%endkatex%} meaning the descent of {%katex inline%}c_{m}{%endkatex%} and {%katex inline%}c_{m}{%endkatex%} being {%katex inline%}c_{0}{%endkatex%} or {%katex inline%}c_{1}{%endkatex%}. {%katex inline%}x^{(i)}_m{%endkatex%} basically means "the x at position (i) in the data for the coefficient m". So for example, {%katex inline%}x^{(i)}_0{%endkatex%} would always be 1, whereas {%katex inline%}x^{(i)}_1{%endkatex%} is the data.

We can implement this more or less the same way we implemented the error function:

```php
function descent(int $m, float $c0, float $c1, array $data): float {
  return (-2 / count($data)) * array_sum(
    array_map(
      function ($point) use ($c0, $c1, $m) {
        return ($point[2] - linearFunction($c0, $point[0], $c1, $point[1])) * $point[$m];
      },
      $data
    )
  );
}
```

This gives us the *descent* that needs to be applied to {%katex inline%}c_m{%endkatex%} in order to bring the linear function closer to the data. Now we can...

# Apply the adjustment

This is rather straight forward: We calculate the descents for {%katex inline%}c_0{%endkatex%} and {%katex inline%}c_1{%endkatex%}, apply some *learning rate* (to adjust the step size) and set the new values.

Wait, now what's a learning rate?

The gradient descent might try to take steps that are way too large to yield good results. Remember the analogy of being stuck on a mountain? Now imagine being 50 meters tall. The smallest step you could take would still be way too large, especially if the valley you want to get in is rather narrow. You wouldn't go into the valley, but end up on the opposite side of it. The next step you would take would be too large, too, so you end up, yet again, on the other side of the valley, eventually going even more uphill than downhill:

![Gradient descent going in the wrong direction, an approximation.](https://dev-to-uploads.s3.amazonaws.com/i/6zwfxjq2g6japfp3jnl5.png)

The learning rate acts as a counter measurement to this: It adjusts your size so your steps are small enough for you to step into the valley.

The learning rate depends on the kind of data you have. A learning rate that's too small will drastically increase the number of steps to take to reach a good fit, a learning rate that's too large will worsen the problem of taking steps that are too large.

There's a sweet spot to every data set, but there's no "one fits all" learning rate. One could apply other ML algorithms to find the perfect learning rate, though, but for small examples, it's best to try a few values and settle for the one that gives the best results.

With this information, we can define two functions to calculate the adjusted values of {%katex inline%}c_0{%endkatex%} and {%katex inline%}c_1{%endkatex%}:

```php
function adaptC0(float $c0, float $c1, array $data, float $learningRate): float {
    return $c0 - $learningRate * descent(0, $c0, $c1, $data);
}

function adaptC1(float $c0, float $c1, array $data, float $learningRate): float {
    return $c1 - $learningRate * descent(1, $c0, $c1, $data);
}
```

Almost there. Now we need a loop to adapt the values until a certain goal is reached. Let's just do 50 iterations and see where the values end up:

```php
// Rather small learning rate
$learningRate = 0.000001;

$errors = [];
for ($i = 0; $i < 50; $i++) {
    // Keep the errors so we can graph them later.
    $errors[] = squaredError($c0, $c1, $data);

    // Do not assign immediately, because otherwise 
    // the result of $c0 would influence the descent 
    // of $c1!
    $newC0 = adaptC0($c0, $c1, $data, $learningRate);
    $newC1 = adaptC1($c0, $c1, $data, $learningRate);

    $c0 = $newC0;
    $c1 = $newC1;
}

echo $c0 . ', ' . $c1; // 14.976594533192, 0.99210801880553
```

Ok, that seems to have worked. Let's see how good those values fit our data:

![Desmos plotting the linear function above the data: It fits well.](https://dev-to-uploads.s3.amazonaws.com/i/vwwc85vjnak58isek260.png)

That's as good as it gets and actually pretty close! This function can now do rather accurate predictions for unknown values of x of the same data set. How did the error do, though?

![Demos showing the errors for each iteration: They get small really fast and then settle near the X axis.](https://dev-to-uploads.s3.amazonaws.com/i/1v98kvqbp52iwr4apper.png)

We see that the error is large at first, but gets smaller really fast: Within around 5 iterations the error stabilizes at around `1429.28` - our implementation works!

# Takeaway thoughts

Linear regression is applied mathematics. Just like most ML algorithms. Machine learning and AI may seem like utter magic at first, but once I started reading into the topic, I realized that the people applying it are basically working with the same tools everyone else does.

It's never the less amazing, that a mathematical model is able to recognize patterns. The more complex algorithms even find patterns a human being probably wouldn't, which serves as a great tool.

---

*I hope you enjoyed reading this article as much as I enjoyed writing it! If so, leave a* ❤️ *or a* 🦄*! I write tech articles in my free time and like to drink coffee every once in a while.*

*If you want to support my efforts, please consider* [*buying me a coffee* ☕](http://buymeacoffee.com/pthormeier) *or* [*follow me on Twitter* 🐦](https://twitter.com/pthormeier)*!*

[![Buy me a coffee button](https://dev-to-uploads.s3.amazonaws.com/i/x55hp6jopwyy161d8e2u.png)](http://buymeacoffee.com/pthormeier)