The model

Multiple linear regression uses multiple features:

ŷ = w1·x1 + w2·x2 + ... + wk·xk + bŷ = w1·x1 + w2·x2 + ... + wk·xk + b

  flowchart LR
  X1[x1] --> M[Linear Model]
  X2[x2] --> M
  Xk[xk] --> M
  M --> Y[Prediction ŷ]

Interpreting coefficients

If all else is equal:

• wkwk tells how much the target changes when feature xkxk increases by 1.

But be careful:

• if features are correlated, coefficient interpretation becomes tricky (multicollinearity)

The Normal Equation (closed-form solution)

Scikit-learn’s LinearRegressionLinearRegression doesn’t guess-and-check its way to the best weights — it can solve for them directly. Hands-On Machine Learning calls this the Normal Equation, a closed-form equation that gives you θθ straight away:

θ̂ = (Xᵀ X)⁻¹ Xᵀ yθ̂ = (Xᵀ X)⁻¹ Xᵀ y

• θ̂θ̂ is the value of θθ that minimizes the MSE cost function
• XX is the training feature matrix (with a leading column of 1s for the bias)
• yy is the vector of target values

import numpy as np
 
# Generate linear-looking data: y = 4 + 3x + noise
np.random.seed(42)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)
 
# add x0 = 1 to every instance
X_b = np.c_[np.ones((100, 1)), X]
 
# theta = (X^T X)^-1 X^T y
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
print("theta:", theta_best.ravel())
 
# predict for two new instances
X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]
y_predict = X_new_b.dot(theta_best)
print("predictions:", y_predict.ravel())

import numpy as np
 
# Generate linear-looking data: y = 4 + 3x + noise
np.random.seed(42)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)
 
# add x0 = 1 to every instance
X_b = np.c_[np.ones((100, 1)), X]
 
# theta = (X^T X)^-1 X^T y
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
print("theta:", theta_best.ravel())
 
# predict for two new instances
X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]
y_predict = X_new_b.dot(theta_best)
print("predictions:", y_predict.ravel())

The values should land close to [4, 3][4, 3] — the function used to generate the data — but not exactly, because of the added noise.

Scikit-learn example

import numpy as np
from sklearn.linear_model import LinearRegression
 
# Example: [size_sqft, bedrooms, age]
X = np.array([
    [800, 2, 10],
    [1000, 3, 5],
    [1200, 3, 20],
    [1500, 4, 7],
])
 
y = np.array([180, 240, 220, 320])
 
model = LinearRegression()
model.fit(X, y)
print("coefficients:", model.coef_)
print("intercept:", model.intercept_)

import numpy as np
from sklearn.linear_model import LinearRegression
 
# Example: [size_sqft, bedrooms, age]
X = np.array([
    [800, 2, 10],
    [1000, 3, 5],
    [1200, 3, 20],
    [1500, 4, 7],
])
 
y = np.array([180, 240, 220, 320])
 
model = LinearRegression()
model.fit(X, y)
print("coefficients:", model.coef_)
print("intercept:", model.intercept_)

Common issues

• multicollinearity: features carry overlapping signal
• scaling: if using regularization, scale inputs

Visualize it

With two features, hθ(x)hθ(x) is no longer a line — it’s a plane sitting over the (x1, x2)(x1, x2) grid. The amber mesh below is that plane; each blue dot is a training instance floating above or below it, and the red stem is its residual. Watch the Normal Equation solve for the plane from scratch each cycle: it starts flat (the “always predict the mean” baseline) and tilts into the least-squares fit as θθ is solved, and the residual stems visibly shrink as it locks into place.

Mini-checkpoint

• Which features are strongly correlated?
• Consider removing or combining them (feature engineering) if needed.

🧪 Try It Yourself

Exercise 1 – Add the Bias Column

Exercise 2 – Solve the Normal Equation

Exercise 3 – Fit Multiple Features with Scikit-Learn

Next

Continue to Polynomial Regression — keep using a linear model, but expand the features themselves to fit curved, non-linear data.