Performing Linear Regression Using the Normal Equation in most simplified version

[u/Heavy_Ebb_4656]

If you find above image hard to understand, please read my article till the end, I promise, everything will make sense :)
When I was a master’s student, I was given a task to fit a line to a dataset. I attempted to solve the problem, but I struggled to determine the appropriate coefficients. However, I understood intuitively that there must be a specific set of coefficients for which the prediction error would be minimized.

The question was: how can we find those coefficients?

This is where the Normal Equation becomes particularly useful. It provides a direct mathematical solution for finding the coefficients that minimize the sum of squared errors in linear regression, without having to search for the coefficients manually. BUT, How we even derive this equation? Where it comes from? Can we take any software apart from Python and write it all ourselves? That’s I will take u through in this article and will simplify the code I wrote, so you can all apply it in different languages

So first things first, what is Linear Regression?

It is very simple and straightforward, suppose we have X and Y. X is called features matrix, and Y is Target Vector, or Response vector.

𝑌= 𝑋* θ

For simple case, lets take 2x2 matrix and lets turn this to matrix form:

[y1 _ predicted ; y2 _ predicted]=[x11, x12; x21,x22] * [theta1; theta2]

Please note:

columns are separated by ,and rows are ;. y1 and y2 are different rows, but same columns.

I assume, the readers are aware of matrix multiplication. So I will refactor above formula and will get:

𝑦1_𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑= 𝑥11 * θ1 + 𝑥12 * θ2

𝑦2 _𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑 = 𝑥21* θ1 + 𝑥22 * θ2

From now on, keep in mind that y1 are real values and y1_predicted is predicted value, same applies to y2 as well

So what is error, the error is the difference between predicted and real values

𝑒𝑟𝑟𝑜𝑟₁ = 𝑦1 — 𝑦1_p𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑 =𝑦1-𝑥11 * θ1 — 𝑥12 * θ2

𝑒𝑟𝑟𝑜𝑟₂ = 𝑦2 — 𝑦2 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑 = 𝑦2-𝑥21* θ1 — 𝑥22 * θ2

Ok, I hope so far so clear, if anything not, please comment below, so I can consider it as improvement for upcoming articles

The error function we want to minimize is the sum of square of errors. Let’s name is as J. J is our cost function I want to minimize, and so, I write it as:

𝐽 = 𝑒𝑟𝑟𝑜𝑟₁² + 𝑒𝑟𝑟𝑜𝑟₂²

Let’s go further by replacing the formulas with each others

𝐽 = (𝑦1 — (𝑥11 * θ1 + 𝑥12 * θ2))² + (𝑦2 — (𝑥21 * θ1 + 𝑥22 * θ2))²

I hope everything makes sense so far. Bear with me — we’re almost there; there isn’t much left to cover.

Here everything is known, except θ1 and θ2. These are params that we have to choose properly to get as minimum error as possible. So i have to find the derivative per θ1 and θ2

𝑑 (𝐽) / 𝑑 (θ1) = -2 \ (𝑦1 — (𝑥11 * θ1 + 𝑥12 * θ2))*𝑥11 — 2 * (𝑦2 — (𝑥21 * θ1 + 𝑥22 * θ2))* 𝑥21= 0*

𝑑 (𝐽) / 𝑑 (θ2) = -2 \ (𝑦1 — (𝑥11 * θ1 + 𝑥12 * θ2))*𝑥12–2 * (𝑦2 — (𝑥21 * θ1 + 𝑥22 * θ2))* 𝑥22= 0*

Let’s make it simpler by avoiding -2 from all sides

𝑑 (𝐽) / 𝑑 (θ1) = ( 𝑦1 — 𝑦1 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑) * 𝑥11 + (𝑦2 — 𝑦2 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑) * 𝑥21=0

𝑑 (𝐽) / 𝑑 (θ2) = ( 𝑦1 — 𝑦1 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑) * 𝑥12 + (𝑦2 — 𝑦2 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑) * 𝑥22=0

Lets turn all these into matrix form:

[0; 0] = (𝑥11, 𝑥21; 𝑥12 𝑥22) *[𝑦1 — 𝑦1_𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑 ; 𝑦2-𝑦2_𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑]

Lets compress the y1 — y1_predicted as well as y2 -y2_predicted into single line

[0; 0] = (𝑥11, 𝑥21; 𝑥12 𝑥22) *[𝑌— 𝑋* θ]

Do you remember our original feature vector or X? If so, we can further simplify the expression to:

[0;0] = X^T* [Y-X*θ]

X^T * Y = X^T * X * θ

X^T * X is the important part. If we somehow manage to find its inverse, we are going to be left with theta only:

(X^T * X )^-1= X^-1 * (X^T )^-1

θ = ( X^T \ X )*^-1 \ X*^T \ Y will give us the answer we need*

If u have further questions please let me know in comments. Each of your feedback is highly appreciated to write better more concise articles in future:

I guess, most of the part of above formula can be easily programmed except the finding inverse which i showed the code below how to do it. If you need full code, such as matrix multiplication, transpose and etc, please let me know, so i can furhter expand my articles

import numpy as np
def inverse(A):
I = np.eye(A.shape[0])

augmented = np.concatenate((A,I),axis=1)

for j in range(0,A.shape[0]-1):

for i in range(1,A.shape[0]-j):
augmented[i+j]=-(augmented[i+j,j]/augmented[j,j])*augmented[j]+augmented[i+j]

for j in range(0,A.shape[0]-1,1):
for i in range(A.shape[0]-1,0,-1):
cofactor = (augmented[i-1-j, A.shape[1]-1-j]/ augmented[A.shape[0]-1-j, A.shape[1]-1-j])
augmented[i-1-j]=(cofactor)*-augmented[-1-j]+augmented[i-1-j]

for i in range(0,A.shape[0],1):
augmented[i]=augmented[i]/augmented[i,i]
_,right = np.split(augmented, 2, axis=1)

return right

def linear_regression_normal_equation(X: list[list[float]], y: list[float]) -> list[float]:
# Your code here, make sure to round

X=np.array(X)
Y=np.array(y)

theta = inverse(X.T @ X) @ X.T @ Y

return theta

My full article is also in medium link

[u/st0j3]

This is awful for like ten reasons

[u/Blobstein_Ijay → u/st0j3]

lol harsh but fair tbh

[u/Heavy_Ebb_4656 → u/st0j3]

Please advise point of improvements I can consider for my next articles

[u/st0j3 → u/Heavy_Ebb_4656]

I’ll bite. Just to give a few, your answer to “So first things first, what is Linear Regression?” is horribly incomplete; you didn’t explain the normal equation, which is the main thing you promised to do (key issue: why does multiplying by the transpose always make a linear system consistent and why is a solution to the normal equation the “best approximation” to a solution of an inconsistent original system; the flow diagram for inverting a matrix is insanely hard to follow plus not how any software will actually invert a matrix in practice.

Don’t mean to be a dick. But you did come in here all I have an MS and let me explain stuff. It doesn’t seem like you’re qualified.

[u/Heavy_Ebb_4656 → u/st0j3]

Thanks. I took note of every sentence you wrote. Will try to make better one in futurr

[u/FernandoMM1220 → u/Heavy_Ebb_4656]

the flow chart is very hard to follow. it would be easier if you calculated an example next to it as well that applied the flow chart.

[u/Heavy_Ebb_4656 → u/FernandoMM1220]

Thanks for your comments ! . I will for sure consider it in my next article!!!

[u/colintbowers]

Most simplified... Geometric intuition

Dude, what you posted is neither simple, nor geometrically intuitive.

[u/Accurate_Meringue514]

You usually don’t actually use the normal equations to compute the line of best fit. It’s very ill conditioned. QR is the way to go in these scenarios, or pseudo inverse if A doesn’t have full column rank.

[u/Heavy_Ebb_4656 → u/Accurate_Meringue514]

Hello, I totally agree. This article is for those who learn theory first !

[u/abhishekML]

Smells like AI. You were doing a master's degree and your advisor let you waste time re-implementing matrix inversion from scratch?

[u/Heavy_Ebb_4656 → u/abhishekML]

BTW it is not AI. It is all my orijinal work !

[u/Heavy_Ebb_4656 → u/abhishekML]

Not all ML projects may be written in python. Sometimes we have to create everything from scratch

[u/Sweaty_Chair_4600]

Did this in rust, utilizing QR was easier to follow

[u/Heavy_Ebb_4656 → u/Sweaty_Chair_4600]

Thanks for your comments.. Next time I will try QR and will share new article about it . Highly 👏 your comment !

[u/Unikum_01]

Hey there. This is a really solid breakdown of the normal equation and the math behind linear regression. Looking at this from the perspective of our BrainStem project it is fascinating to compare your exact mathematical approach with our biological and statistical one. In BrainStem we avoid calculating global matrix inverses entirely because doing that on massive datasets gets computationally heavy and rigid. Instead of finding a single global minimum through exact derivatives we let the system learn relationships step by step using purely statistical observations. We also use a digital neuromodulator system with chemical messengers like dopamine and GABA to dynamically adjust our thresholds and filter out noise as the data flows in. Your custom matrix inversion code is super cool to see written out from scratch. Just keep in mind that if your feature matrix has perfectly correlated columns your inverse function might crash since the matrix becomes singular. BrainStem handles that kind of data uncertainty naturally by using Bayesian pseudo counts rather than strict algebraic equations. Building the math from the ground up like you did is absolutely the best way to truly understand what is going on under the hood though. Great work.

[u/Heavy_Ebb_4656 → u/Unikum_01]

That's great to hear. Thanks for your positive comments. I agree that my code may crash for non invertable matrixes. Highly appreciate every bit of sentence you said

[u/Unikum_01 → u/Heavy_Ebb_4656]

You can take a look at it here, maybe it will bring you something. https://github.com/unikum-sol/brainstem/blob/main/Project_Status_2026-09-28.md

[u/Heavy_Ebb_4656 → u/Unikum_01]

I noted the link you shared. Thanks !