When people say:

“I used machine learning to predict the price using several features.”

it can sound like something far more complicated than the simple line we started with.

But if we strip away the terminology, multiple linear regression is the same mathematical problem as simple linear regression, with more inputs.

At its core, we are still just trying to find a handful of numbers:

The intercept — $\beta_0$
One coefficient per input — $\beta_1, \beta_2, \dots, \beta_p$

Once we have those numbers, making a prediction is simply a matter of putting values into an equation.

This article builds on simple linear regression. If you have not read that one yet, start there. Here we reuse the same ideas: the error, the SSE, partial derivatives, and setting them to zero.

1. The basic idea

Suppose we want to predict someone’s salary. Years of experience alone is a good start, but it is not the whole story. Let’s add a second input: the number of professional certifications they hold.

We have some historical data (salary in dollars):

Years of Experience	Certifications	Salary
1	0	24,500
2	1	30,300
3	0	32,300
4	2	41,600
5	1	43,600
6	3	52,800
7	2	53,500
8	3	61,400

With one input, we fitted a line. With two inputs, we fit a plane. With three or more, we fit a hyperplane, which is hard to draw but works exactly the same way.

The equation of our model is:

\hat{y} = \beta_0 + \beta_1 x_1 + \beta_2 x_2

Where:

$\hat{y}$ = predicted salary
$x_1$ = years of experience
$x_2$ = number of certifications
$\beta_0$ = intercept
$\beta_1$ = coefficient of years of experience
$\beta_2$ = coefficient of certifications

For any number of inputs $p$ , the general form is:

\hat{y} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p

That’s essentially our entire prediction model.

2. What does each coefficient mean?

Each coefficient tells us how much the predicted output changes when that one input increases by one unit, while all the other inputs stay the same.

That last part is the important difference from simple regression.

For example, suppose we calculate:

\beta_1 = 4{,}161

This means that for every additional year of experience, with the number of certifications held fixed, the predicted salary increases by approximately 4,161.

Similarly:

\beta_2 = 2{,}596

means that each additional certification adds approximately 2,596 to the predicted salary, with years of experience held fixed.

So if our model is:

\hat{y} = 19{,}881 + 4{,}161x_1 + 2{,}596x_2

then:

3 years, 0 certifications → 32,364
4 years, 0 certifications → 36,525
4 years, 1 certification → 39,121

Going from the first row to the second changes only experience, and the prediction moves by $\beta_1$ . Going from the second row to the third changes only certifications, and the prediction moves by $\beta_2$ .

Each coefficient controls how steep the plane is in the direction of its own input.

3. What does the intercept mean?

The intercept is the predicted value of $y$ when all inputs are zero.

Suppose:

\beta_0 = 19{,}881

Then our model says:

\hat{y} = 19{,}881 + 4{,}161x_1 + 2{,}596x_2

When experience and certifications are both zero:

\hat{y} = 19{,}881 + 4{,}161(0) + 2{,}596(0)

Therefore:

\hat{y} = 19{,}881

So the intercept is:

The point where the regression plane crosses the $y$ -axis.

Depending on the problem, the intercept may or may not have a meaningful real-world interpretation. Mathematically, however, it is one of the parameters needed to define the plane.

4. So where do the coefficients come from?

This is where the interesting mathematics happens.

We don’t simply guess:

\beta_0 = 19{,}881, \quad \beta_1 = 4{,}161, \quad \beta_2 = 2{,}596

Instead, we calculate values that produce the best-fitting plane for the observed data.

We use the same approach as before: Ordinary Least Squares (OLS).

The idea behind OLS is unchanged:

Find the plane that minimizes the total squared difference between the actual values and the predicted values.

For each observation, we have an error:

e_i = y_i - \hat{y}_i

Since some errors are positive and others are negative, simply adding them together could cause them to cancel each other out.

Instead, we square each error.

The objective becomes:

\text{SSE} = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2

Since:

\hat{y}_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}

we can write:

\text{SSE} = \sum_{i=1}^{n} \left( y_i - \beta_0 - \beta_1 x_{i1} - \beta_2 x_{i2} \right)^2

Here $x_{i1}$ means the value of input 1 for observation $i$ , and $x_{i2}$ is the value of input 2 for the same observation.

OLS finds the values of $\beta_0$ , $\beta_1$ and $\beta_2$ that minimize this quantity.

5. Deriving the coefficients using partial derivatives

We have already defined the Sum of Squared Errors:

\text{SSE} = \sum_{i=1}^{n} \left(y_i - \beta_0 - \beta_1 x_{i1} - \beta_2 x_{i2}\right)^2

Our goal is to find the values of $\beta_0$ , $\beta_1$ and $\beta_2$ that make this SSE as small as possible.

Because we now have three unknowns, we take the partial derivative with respect to each one. That gives us three equations, one for each unknown.

Step 1: Find $\beta_0$

We differentiate with respect to $\beta_0$ . The summation can stay outside, and we use the chain rule. If:

u = y_i - \beta_0 - \beta_1 x_{i1} - \beta_2 x_{i2}

then:

\frac{\partial}{\partial \beta_0}\left(u^2\right) = 2u \, \frac{\partial u}{\partial \beta_0}

Since:

\frac{\partial u}{\partial \beta_0} = -1

we get:

\frac{\partial \text{SSE}}{\partial \beta_0} = -2 \sum_{i=1}^{n} \left(y_i - \beta_0 - \beta_1 x_{i1} - \beta_2 x_{i2}\right)

At the minimum, the derivative must be zero:

2 $\sum_{i=1}^{n} \left(y_i - \beta_0 - \beta_1 x_{i1} - \beta_2 x_{i2}\right) = 0$

Divide by $-2$ and expand the summation:

\sum y_i - n\beta_0 - \beta_1 \sum x_{i1} - \beta_2 \sum x_{i2} = 0

Divide by $n$ :

\bar{y} - \beta_0 - \beta_1 \bar{x}_1 - \beta_2 \bar{x}_2 = 0

Therefore:

\boxed{\beta_0 = \bar{y} - \beta_1 \bar{x}_1 - \beta_2 \bar{x}_2}

Compare this with simple regression, where we had $\beta_0 = \bar{y} - \beta_1 \bar{x}$ . The idea is identical: the fitted plane always passes through the point of means.

Step 2: Find $\beta_1$

Now we differentiate the same SSE with respect to $\beta_1$ . Using the chain rule again, the derivative of the inside is:

\frac{\partial u}{\partial \beta_1} = -x_{i1}

Therefore:

\frac{\partial \text{SSE}}{\partial \beta_1} = -2 \sum_{i=1}^{n} x_{i1} \left(y_i - \beta_0 - \beta_1 x_{i1} - \beta_2 x_{i2}\right)

Set it to zero, divide by $-2$ , and expand:

\sum x_{i1} y_i - \beta_0 \sum x_{i1} - \beta_1 \sum x_{i1}^2 - \beta_2 \sum x_{i1} x_{i2} = 0

Step 3: Find $\beta_2$

The same steps for $\beta_2$ , where the derivative of the inside is $-x_{i2}$ , give:

\sum x_{i2} y_i - \beta_0 \sum x_{i2} - \beta_1 \sum x_{i1} x_{i2} - \beta_2 \sum x_{i2}^2 = 0

Step 4: The normal equations

Now we collect all three equations, moving the terms with $\beta$ to the left:

\begin{aligned} n\beta_0 + \beta_1 \sum x_{i1} + \beta_2 \sum x_{i2} &= \sum y_i \\ \beta_0 \sum x_{i1} + \beta_1 \sum x_{i1}^2 + \beta_2 \sum x_{i1} x_{i2} &= \sum x_{i1} y_i \\ \beta_0 \sum x_{i2} + \beta_1 \sum x_{i1} x_{i2} + \beta_2 \sum x_{i2}^2 &= \sum x_{i2} y_i \end{aligned}

These are called the normal equations. There are three equations and three unknowns, so we can solve for $\beta_0$ , $\beta_1$ and $\beta_2$ .

Notice what changed compared with simple regression. With one input, we could solve by hand with a couple of substitutions and get a tidy formula. With more inputs, the equations are coupled: every coefficient appears in every equation, through terms like $\sum x_{i1} x_{i2}$ . That cross term measures how the inputs move together.

Setting the derivatives to zero only finds a stationary point, but here it is the minimum: SSE is a convex function of the coefficients (a sum of squared linear terms), so when the inputs are not perfectly redundant, its only stationary point is the global minimum.

6. From the normal equations to the matrix formula

At this point, we have derived three normal equations:

n\beta_0+\beta_1\sum x_{i1}+\beta_2\sum x_{i2} = \sum y_i

\beta_0\sum x_{i1} +\beta_1\sum x_{i1}^2 +\beta_2\sum x_{i1}x_{i2} = \sum x_{i1}y_i

\beta_0\sum x_{i2} +\beta_1\sum x_{i1}x_{i2} +\beta_2\sum x_{i2}^2 = \sum x_{i2}y_i

These equations are correct, but there is a problem: as soon as we add more input variables, the number of terms grows quickly.

For three features, we would have four equations. For ten features, we would have eleven equations. Writing every equation separately becomes difficult to read and work with.

This is where matrices become useful.

The important point is that the matrix formula does not introduce a new idea. It is simply a compact way of writing the same normal equations.

6.1 Put the data into a matrix

Suppose we have two input features, $x_1$ and $x_2$ , and $n$ observations.

Our regression model is:

\hat{y}_i=\beta_0+\beta_1x_{i1}+\beta_2x_{i2}

We put all observations into a matrix $X$ :

X= \begin{bmatrix} 1 & x_{11} & x_{12}\\ 1 & x_{21} & x_{22}\\ 1 & x_{31} & x_{32}\\ \vdots & \vdots & \vdots\\ 1 & x_{n1} & x_{n2} \end{bmatrix}

The first column contains $1$ 's because we have an intercept $\beta_0$ .

The coefficient vector is:

\beta= \begin{bmatrix} \beta_0\\ \beta_1\\ \beta_2 \end{bmatrix}

and the observed target values are:

y= \begin{bmatrix} y_1\\ y_2\\ y_3\\ \vdots\\ y_n \end{bmatrix}

6.2 Why does $X\beta$ represent all predictions?

Let's multiply $X$ by $\beta$ :

X\beta= \begin{bmatrix} 1 & x_{11} & x_{12}\\ 1 & x_{21} & x_{22}\\ 1 & x_{31} & x_{32}\\ \vdots & \vdots & \vdots\\ 1 & x_{n1} & x_{n2} \end{bmatrix} \begin{bmatrix} \beta_0\\ \beta_1\\ \beta_2 \end{bmatrix}

Matrix multiplication gives:

X\beta= \begin{bmatrix} \beta_0+\beta_1x_{11}+\beta_2x_{12}\\ \beta_0+\beta_1x_{21}+\beta_2x_{22}\\ \beta_0+\beta_1x_{31}+\beta_2x_{32}\\ \vdots\\ \beta_0+\beta_1x_{n1}+\beta_2x_{n2} \end{bmatrix}

But each row is exactly one prediction:

\hat{y}_i=\beta_0+\beta_1x_{i1}+\beta_2x_{i2}

Therefore:

\boxed{\hat{y}=X\beta}

This is the first major benefit of matrix notation: instead of writing every prediction separately, $X\beta$ calculates all of them at once.

6.3 Write the SSE using matrices

The prediction errors are:

y-X\beta

because $y$ contains the actual values and $X\beta$ contains the predicted values.

The ordinary sum of squared errors is:

SSE=\sum_{i=1}^{n}(y_i-\hat{y}_i)^2

We can write exactly the same thing using a vector:

\boxed{ SSE=(y-X\beta)^T(y-X\beta) }

Why?

Suppose the error vector is:

e=y-X\beta= \begin{bmatrix} e_1\\ e_2\\ \vdots\\ e_n \end{bmatrix}

Then:

e^Te= \begin{bmatrix} e_1&e_2&\cdots&e_n \end{bmatrix} \begin{bmatrix} e_1\\ e_2\\ \vdots\\ e_n \end{bmatrix}

which gives:

e_1^2+e_2^2+\cdots+e_n^2

Therefore:

e^Te=\sum_{i=1}^{n}e_i^2

So:

(y-X\beta)^T(y-X\beta)

is simply the matrix version of our original SSE formula.

6.4 Expand the matrix expression

Starting with:

SSE=(y-X\beta)^T(y-X\beta)

expand it just like ordinary algebra:

SSE = y^Ty-y^TX\beta-\beta^TX^Ty+\beta^TX^TX\beta

The two middle terms are scalars and are equal:

y^TX\beta=\beta^TX^Ty

so we can combine them:

\boxed{ SSE=y^Ty-2\beta^TX^Ty+\beta^TX^TX\beta }

Now we have the entire least-squares objective expressed using matrices.

6.5 Differentiate with respect to $\beta$

We want the values of $\beta$ that minimize SSE.

Taking the derivative with respect to the coefficient vector gives:

\frac{\partial SSE}{\partial\beta} = -2X^Ty+2X^TX\beta

At the minimum, the derivative is zero:

$2X^Ty+2X^TX\beta=0$

Move the first term to the other side:

2X^TX\beta=2X^Ty

Divide by $2$ :

\boxed{ X^TX\beta=X^Ty }

This is the matrix form of the normal equations.

But where did those three original equations go?

They are still there.

The matrix equation is simply all three equations written together.

6.6 Let's expand $X^TX$

This is the step that makes the connection clear.

Recall:

X= \begin{bmatrix} 1 & x_{11} & x_{12}\\ 1 & x_{21} & x_{22}\\ \vdots & \vdots & \vdots\\ 1 & x_{n1} & x_{n2} \end{bmatrix}

Therefore:

X^T= \begin{bmatrix} 1 & 1 & \cdots & 1\\ x_{11} & x_{21} & \cdots & x_{n1}\\ x_{12} & x_{22} & \cdots & x_{n2} \end{bmatrix}

Now multiply $X^T$ by $X$ :

X^TX= \begin{bmatrix} 1 & 1 & \cdots & 1\\ x_{11} & x_{21} & \cdots & x_{n1}\\ x_{12} & x_{22} & \cdots & x_{n2} \end{bmatrix} \begin{bmatrix} 1 & x_{11} & x_{12}\\ 1 & x_{21} & x_{22}\\ \vdots & \vdots & \vdots\\ 1 & x_{n1} & x_{n2} \end{bmatrix}

Let's calculate each entry.

The top-left entry is:

1^2+1^2+\cdots+1^2=n

The top-middle entry is:

x_{11}+x_{21}+\cdots+x_{n1} = \sum x_{i1}

The top-right entry is:

x_{12}+x_{22}+\cdots+x_{n2} = \sum x_{i2}

Continuing this multiplication gives:

\boxed{ X^TX= \begin{bmatrix} n & \sum x_{i1} & \sum x_{i2}\\ \sum x_{i1} & \sum x_{i1}^2 & \sum x_{i1}x_{i2}\\ \sum x_{i2} & \sum x_{i1}x_{i2} & \sum x_{i2}^2 \end{bmatrix} }

Notice something important: these are exactly the coefficients that appeared in our three normal equations.

6.7 Now expand $X^Ty$

We also need the right-hand side of:

X^TX\beta=X^Ty

Calculate:

X^Ty= \begin{bmatrix} 1 & 1 & \cdots & 1\\ x_{11} & x_{21} & \cdots & x_{n1}\\ x_{12} & x_{22} & \cdots & x_{n2} \end{bmatrix} \begin{bmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{bmatrix}

This gives:

\boxed{ X^Ty= \begin{bmatrix} \sum y_i\\ \sum x_{i1}y_i\\ \sum x_{i2}y_i \end{bmatrix} }

Again, these are exactly the right-hand sides of our three normal equations.

6.8 Put everything together

We now have:

X^TX= \begin{bmatrix} n & \sum x_{i1} & \sum x_{i2}\\ \sum x_{i1} & \sum x_{i1}^2 & \sum x_{i1}x_{i2}\\ \sum x_{i2} & \sum x_{i1}x_{i2} & \sum x_{i2}^2 \end{bmatrix}

and:

X^Ty= \begin{bmatrix} \sum y_i\\ \sum x_{i1}y_i\\ \sum x_{i2}y_i \end{bmatrix}

Therefore:

X^TX\beta=X^Ty

means:

\begin{bmatrix} n & \sum x_{i1} & \sum x_{i2}\\ \sum x_{i1} & \sum x_{i1}^2 & \sum x_{i1}x_{i2}\\ \sum x_{i2} & \sum x_{i1}x_{i2} & \sum x_{i2}^2 \end{bmatrix} \begin{bmatrix} \beta_0\\ \beta_1\\ \beta_2 \end{bmatrix} = \begin{bmatrix} \sum y_i\\ \sum x_{i1}y_i\\ \sum x_{i2}y_i \end{bmatrix}

Now multiply the matrices.

The first row gives:

n\beta_0 + \beta_1\sum x_{i1} + \beta_2\sum x_{i2} = \sum y_i

The second row gives:

\beta_0\sum x_{i1} + \beta_1\sum x_{i1}^2 + \beta_2\sum x_{i1}x_{i2} = \sum x_{i1}y_i

The third row gives:

\beta_0\sum x_{i2} + \beta_1\sum x_{i1}x_{i2} + \beta_2\sum x_{i2}^2 = \sum x_{i2}y_i

These are exactly our original normal equations.

So we haven't changed the mathematics at all.

We've simply compressed:

\boxed{\text{three equations}}

into:

\boxed{\text{one matrix equation}}

6.9 From the normal equations to the closed-form solution

Now we have:

X^TX\beta=X^Ty

This has the same structure as an ordinary algebra equation:

3x=12

To solve for $x$ , we divide by $3$ :

x=\frac{12}{3}

For matrices, the equivalent of "dividing by $A$ " is multiplying by its inverse $A^{-1}$ .

So multiply both sides by $(X^TX)^{-1}$ :

(X^TX)^{-1}X^TX\beta = (X^TX)^{-1}X^Ty

By definition:

(X^TX)^{-1}(X^TX)=I

where $I$ is the identity matrix.

Therefore:

I\beta=(X^TX)^{-1}X^Ty

and since:

I\beta=\beta

we finally get:

\boxed{ \beta=(X^TX)^{-1}X^Ty }

This is the closed-form solution for ordinary least squares linear regression.

6.10 What the formula is really saying

The formula:

\boxed{ \beta=(X^TX)^{-1}X^Ty }

may look complicated, but each part has a simple purpose:

$X$ contains our input data.
$X^T$ turns rows into columns so we can perform the required matrix multiplications.
$X^TX$ collects the relationships between the input variables, including the cross-products such as $x_1x_2$ .
$X^Ty$ collects the relationships between each input variable and the target.
$(X^TX)^{-1}$ solves the resulting system of equations.
$\beta$ contains the optimal coefficients.

So the entire derivation can be summarized as:

\boxed{ \text{Minimize SSE} }

\Downarrow

\boxed{ \frac{\partial SSE}{\partial\beta}=0 }

\Downarrow

\boxed{ X^TX\beta=X^Ty }

\Downarrow

\boxed{ \beta=(X^TX)^{-1}X^Ty }

The important insight is that the matrix formula did not appear out of nowhere.

We first derived the normal equations using calculus. Matrix algebra then gave us a compact way to represent and solve all of those equations simultaneously.

6.11 One practical note

Although the formula is commonly written as:

\beta=(X^TX)^{-1}X^Ty

real machine-learning libraries generally do not explicitly calculate the inverse of $X^TX$ . Direct matrix inversion can be numerically less stable and inefficient for large datasets.

Instead, implementations commonly use numerical methods such as QR decomposition or singular value decomposition (SVD) to solve:

X^TX\beta=X^Ty

or an equivalent least-squares problem.

For understanding the mathematics of linear regression, however, the closed-form expression

\boxed{\beta=(X^TX)^{-1}X^Ty}

is extremely useful because it exposes exactly how the optimal coefficients are obtained.

This is the OLS solution for any number of inputs.

The convexity argument also has a clean matrix version. The second derivative of the SSE is $2X^\top X$ , which is positive definite when the columns of $X$ are linearly independent. That guarantees a single, unique minimum.

If one input is a perfect combination of the others, for example $x_2 = 2x_1$ , then $X^\top X$ cannot be inverted and there is no unique answer. This problem is called perfect multicollinearity.

7. The whole idea in one flow

So the entire process is:

\text{Prediction equation} \rightarrow \text{Residual} \rightarrow \text{SSE} \rightarrow \text{Partial derivatives} \rightarrow \text{Set derivatives to zero} \rightarrow \boldsymbol{\beta}

Specifically:

\hat{y} = \beta_0 + \beta_1 x_1 + \beta_2 x_2

Then:

\text{SSE} = \sum_{i=1}^{n} \left(y_i - \beta_0 - \beta_1 x_{i1} - \beta_2 x_{i2}\right)^2

Differentiate with respect to each coefficient and set to zero:

\frac{\partial \text{SSE}}{\partial \beta_0} = 0, \quad \frac{\partial \text{SSE}}{\partial \beta_1} = 0, \quad \frac{\partial \text{SSE}}{\partial \beta_2} = 0

which gives the normal equations, and in matrix form:

\boxed{\boldsymbol{\beta} = \left(X^\top X\right)^{-1} X^\top \mathbf{y}}

Therefore, multiple linear regression is not simply guessing a plane.

It uses calculus to find the values of $\beta_0, \beta_1, \dots, \beta_p$ that produce the minimum possible squared error for the training data.

8. Working through the example

Let’s apply this to our salary data. First, we compute the sums we need. With $n = 8$ :

Quantity	Value
$\sum x_{i1}$	36
$\sum x_{i2}$	12
$\sum x_{i1}^2$	204
$\sum x_{i2}^2$	28
$\sum x_{i1} x_{i2}$	71
$\sum y_i$	340,000
$\sum x_{i1} y_i$	1,748,900
$\sum x_{i2} y_i$	606,700

Plugging these into the normal equations gives:

\begin{aligned} 8\beta_0 + 36\beta_1 + 12\beta_2 &= 340{,}000 \\ 36\beta_0 + 204\beta_1 + 71\beta_2 &= 1{,}748{,}900 \\ 12\beta_0 + 71\beta_1 + 28\beta_2 &= 606{,}700 \end{aligned}

In matrix form, this is $X^\top X \boldsymbol{\beta} = X^\top \mathbf{y}$ with:

X^\top X = \begin{bmatrix} 8 & 36 & 12 \\ 36 & 204 & 71 \\ 12 & 71 & 28 \end{bmatrix}, \quad X^\top \mathbf{y} = \begin{bmatrix} 340{,}000 \\ 1{,}748{,}900 \\ 606{,}700 \end{bmatrix}

Solving this system of three equations gives:

\boxed{\beta_0 \approx 19{,}881}, \quad \boxed{\beta_1 \approx 4{,}161}, \quad \boxed{\beta_2 \approx 2{,}596}

We can check the intercept with the formula from Step 1. The means are $\bar{x}_1 = 4.5$ , $\bar{x}_2 = 1.5$ and $\bar{y} = 42{,}500$ , so:

\beta_0 = 42{,}500 - 4{,}161(4.5) - 2{,}596(1.5) \approx 19{,}881

It matches.

Now let’s look at how well the fitted plane matches the data. The chart below plots each employee’s actual salary against the salary the model predicts. If the model were perfect, every point would sit on the dashed line.

Actual salary versus predicted salary for the multiple regression model, with points close to the perfect-prediction line

9. Now prediction becomes extremely simple

Our calculations give us:

\hat{y} = 19{,}881 + 4{,}161x_1 + 2{,}596x_2

Now suppose someone has 6 years of experience and 2 certifications.

We simply substitute:

x_1 = 6, \quad x_2 = 2

Therefore:

19{,}881 + (4{,}161)(6) + (2{,}596)(2)

19{,}881 + 24{,}966 + 5{,}192 = 50{,}039

That’s the prediction.

There is no mysterious “AI thinking” happening here.

The model has learned three numbers:

\boxed{\beta_0 = 19{,}881}, \quad \boxed{\beta_1 = 4{,}161}, \quad \boxed{\beta_2 = 2{,}596}

Then it applies the equation.

10. What changes when we add more inputs?

The mathematics is the same, but a few ideas deserve attention.

A coefficient depends on which other inputs are in the model.

If we fit salary using years of experience alone, the slope comes out at about 5,212 per year. With certifications added, the experience coefficient drops to about 4,161.

Why? In our data, people with more experience also tend to hold more certifications. In the simple model, the experience slope absorbs part of the certification effect. Once certifications are included as their own input, the model can separate the two effects.

This is why we say each coefficient is the effect of one input holding the others constant.

Inputs that move together cause trouble.

When two inputs are strongly correlated, $X^\top X$ is close to being non-invertible. The model may still predict well, but the individual coefficients become unstable: small changes in the data can swing them widely. This is called multicollinearity.

Inputs on different scales are hard to compare.

A coefficient of 4,161 for years and 2,596 for certifications does not tell us which input matters more, because they are measured in different units. Standardizing the inputs is a common way to make coefficients comparable.

11. So what does `model.fit(X, y)` actually do?

This is where the machine-learning terminology can make things sound more complicated than they are.

In Python, we might write:

import pandas as pd
from sklearn.linear_model import LinearRegression

df = pd.DataFrame({
    "experience":     [1, 2, 3, 4, 5, 6, 7, 8],
    "certifications": [0, 1, 0, 2, 1, 3, 2, 3],
    "salary": [24500, 30300, 32300, 41600, 43600, 52800, 53500, 61400],
})

X = df[["experience", "certifications"]]
y = df["salary"]

model = LinearRegression()

model.fit(X, y)

When we call:

model.fit(X, y)

we are essentially asking the algorithm:

“Look at these examples and find the parameters that give the best-fitting linear relationship.”

For multiple linear regression, those parameters are:

\beta_0, \quad \beta_1, \quad \dots, \quad \beta_p

After fitting the model, we can inspect them:

print(model.intercept_)   # about 19881
print(model.coef_)        # about [4161, 2596]

These are the same numbers we found by hand.

Conceptually:

model.fit(X, y)

        ↓

Find β₀, β₁, ..., βₚ

        ↓

ŷ = β₀ + β₁x₁ + ... + βₚxₚ

        ↓

Use the equation for prediction

Scikit-learn uses numerical linear-algebra methods internally rather than necessarily evaluating $\left(X^\top X\right)^{-1} X^\top \mathbf{y}$ literally. Inverting a matrix directly can be numerically unstable, so practical solvers use more stable factorizations. But mathematically, for ordinary multiple linear regression with an intercept, it is solving the same least-squares problem.

12. Why do we call this “Machine Learning”?

This is perhaps the most interesting part.

The mathematics itself is not new.

Linear regression has a long history in statistics and mathematics.

What makes it part of a machine-learning workflow is the idea that:

We provide data.
An algorithm estimates parameters from that data.
The resulting model is used to make predictions on new data.

For multiple linear regression, the learned parameters happen to be:

\beta_0, \quad \beta_1, \quad \dots, \quad \beta_p

For more complicated machine-learning models, there may be thousands, millions, or even billions of parameters.

But the basic idea remains:

Use data to estimate parameters, then use those parameters to produce predictions.

13. The important mental model

When you hear:

“I trained a multiple regression model.”

don’t imagine a mysterious black box.

For multiple linear regression, imagine this:

                 Training Data
                      │
                      ▼
          Find the best β₀, β₁, ..., βₚ
                      │
                      ▼
          ŷ = β₀ + β₁x₁ + ... + βₚxₚ
                      │
                      ▼
               New inputs x₁ ... xₚ
                      │
                      ▼
                  Prediction ŷ

The “learning” part is essentially finding the parameters.

The “prediction” part is simply using those parameters in the equation.

14. The bigger picture

Multiple linear regression is a beautiful example because it shows that going from one input to many changes the bookkeeping, not the idea.

At the surface, we might write:

model.fit(X, y)

and:

model.predict(X_new)

But underneath those convenient APIs is mathematics.

The fundamental process is:

\boxed{\text{Data} \rightarrow \text{Optimization} \rightarrow \text{Parameters} \rightarrow \text{Prediction}}

For multiple linear regression:

\boxed{\text{Optimization} \rightarrow \beta_0, \beta_1, \dots, \beta_p}

and then:

\boxed{\hat{y} = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p}

That’s really what is happening behind the scenes.

Final takeaway

The next time someone says:

“I used machine learning to predict something using several features.”

remember what is actually happening.

There is a dataset.

We find the plane that best fits the data using a mathematical optimization method such as Ordinary Least Squares.

That gives us the parameters:

\boxed{\boldsymbol{\beta} = \left(X^\top X\right)^{-1} X^\top \mathbf{y}}

Then we put them into:

\boxed{\hat{y} = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p}

And that’s the prediction.

The “machine learning” label doesn’t make the mathematics disappear. In this case, the machine is essentially helping us estimate the parameters from data. The prediction itself is just the equation.

Share this article

What Really Happens Behind the Scenes in Multiple Linear Regression?

1. The basic idea

2. What does each coefficient mean?

3. What does the intercept mean?

4. So where do the coefficients come from?

5. Deriving the coefficients using partial derivatives

Step 1: Find β0\beta_0β0​

Step 2: Find β1\beta_1β1​

Step 3: Find β2\beta_2β2​

Step 4: The normal equations

6. From the normal equations to the matrix formula

6.1 Put the data into a matrix

6.2 Why does XβX\betaXβ represent all predictions?

6.3 Write the SSE using matrices

6.4 Expand the matrix expression

6.5 Differentiate with respect to β\betaβ

6.6 Let's expand XTXX^TXXTX

6.7 Now expand XTyX^TyXTy

6.8 Put everything together

6.9 From the normal equations to the closed-form solution

6.10 What the formula is really saying

6.11 One practical note

7. The whole idea in one flow

8. Working through the example

9. Now prediction becomes extremely simple

10. What changes when we add more inputs?

11. So what does model.fit(X, y) actually do?

12. Why do we call this “Machine Learning”?

13. The important mental model

14. The bigger picture

Final takeaway

Step 1: Find $\beta_0$

Step 2: Find $\beta_1$

Step 3: Find $\beta_2$

6.2 Why does $X\beta$ represent all predictions?

6.5 Differentiate with respect to $\beta$

6.6 Let's expand $X^TX$

6.7 Now expand $X^Ty$

11. So what does `model.fit(X, y)` actually do?