The best fit, and what "best" was defined to mean

The problem with no exact answer

You have 500 houses, four features each, and one price. You want X w = y, where X is 500×4 and w is the four coefficients you are looking for. There are 500 equations and 4 unknowns. Unless the data is astonishingly obliging, no w satisfies all 500. There is no exact answer.

So you settle for the closest one, and the definition of "closest" is a choice. Least squares chooses to minimise the sum of squared residuals:

minimise  ||X w - y||^2  =  sum over rows of (prediction - actual)^2

That is a decision, not a law. Choosing squares rather than absolute values means one house that is ₹50 lakh out counts as much as a hundred houses ₹5 lakh out, because 50² = 2,500 and 100 × 5² = 2,500. Whether that is right depends on whether one big mistake really is as bad as a hundred small ones for your purpose. Frequently it is not, and that is the argument for absolute-error regression instead.

The geometry, which makes the formula obvious

Here is the picture worth carrying. X w — as w varies over all possibilities — traces out the column space of X: a 4-dimensional flat surface sitting inside 500-dimensional space. Your target y is a point in that 500-dimensional space, and almost certainly not on the surface.

The closest point on a flat surface to an outside point is found by dropping a perpendicular. So the best Xw is the projection of y onto the column space, and the residual y − Xw must be perpendicular to that surface — which means perpendicular to every column of X:

X^T (y - X w) = 0

Rearrange:

X^T X w = X^T y            <- the normal equations
w = (X^T X)^-1 X^T y       <- the closed form

That is the whole derivation, and it comes from one geometric fact: the error must be orthogonal to everything the model could have produced. The word "normal" in "normal equations" means perpendicular, not typical.

Worked, small

Fit y = w x through three points (1, 2), (2, 4.1), (3, 5.8), no intercept.

X = [1; 2; 3]        y = [2; 4.1; 5.8]
X^T X = 1 + 4 + 9 = 14
X^T y = 1*2 + 2*4.1 + 3*5.8 = 2 + 8.2 + 17.4 = 27.6
w = 27.6 / 14 = 1.971

Predictions: 1.971, 3.943, 5.914. Residuals: +0.029, +0.157, −0.114. Their dot product with the column [1, 2, 3] is 0.029 + 0.314 − 0.343 = 0, up to rounding. The residual really is perpendicular to the column space. The geometry is not a metaphor.

What to run in practice

Do not code the closed form. Use np.linalg.lstsq, which solves it via QR decomposition or SVD, is numerically far better behaved, and works even when X^T X is singular:

import numpy as np
w, residuals, rank, sv = np.linalg.lstsq(X, y, rcond=None)

Note what it returns alongside the answer: the rank of X and its singular values. Those tell you whether your features were linearly dependent, which is the failure the closed form hides.

Reading R-squared without being fooled

The standard summary statistic is

R^2 = 1 - (sum of squared residuals) / (sum of squared deviations from the mean of y)

It says what fraction of the variance in y your model accounts for, compared with predicting the mean every time. R² = 0 means you did no better than the mean; R² = 1 means perfect fit.

Two honest warnings. First, R² never decreases when you add a feature, even a column of random numbers, so comparing models by R² always favours the bigger one. Adjusted R² penalises this, imperfectly. Second, R² is computed on the data you fitted; the only version worth quoting is computed on data the model has not seen.

The assumptions you have quietly made

Least squares always returns an answer, so it is worth naming what has to be true for the answer to be meaningful.

• Linearity. You have assumed y is a straight-line function of the features. Fit a line through a curve and you get a line, with no complaint from the arithmetic.
• Independent errors. With time-series data, consecutive residuals are usually correlated, and the coefficient standard errors will be too small — your significance tests become overconfident.
• Roughly constant error size. If the errors grow with the prediction, as they do for prices and counts, the fit is dominated by the large end. Taking logs of the target often fixes this in one line.
• No influential outliers. Squaring makes a single wild point capable of dragging the whole fit. Always plot residuals against predictions; the plot shows all four of these problems and takes ten seconds.

The rule to keep

Least squares projects your data onto what your model can represent, and calls the perpendicular distance the error. If your model cannot represent the truth, the projection is still computed, still optimal, and still wrong — optimality is relative to a model class you chose.