Normal Equation Implementation in Python / Numpy

Question

I've written some beginner code to calculate the co-efficients of a simple linear model using the normal equation.

# Modules
import numpy as np

# Loading data set
X, y = np.loadtxt('ex1data3.txt', delimiter=',', unpack=True)

data = np.genfromtxt('ex1data3.txt', delimiter=',')

def normalEquation(X, y):
    m = int(np.size(data[:, 1]))

    # This is the feature / parameter (2x2) vector that will
    # contain my minimized values
    theta = []

    # I create a bias_vector to add to my newly created X vector
    bias_vector = np.ones((m, 1))

    # I need to reshape my original X(m,) vector so that I can
    # manipulate it with my bias_vector; they need to share the same
    # dimensions.
    X = np.reshape(X, (m, 1))

    # I combine these two vectors together to get a (m, 2) matrix
    X = np.append(bias_vector, X, axis=1)

    # Normal Equation:
    # theta = inv(X^T * X) * X^T * y

    # For convenience I create a new, tranposed X matrix
    X_transpose = np.transpose(X)

    # Calculating theta
    theta = np.linalg.inv(X_transpose.dot(X))
    theta = theta.dot(X_transpose)
    theta = theta.dot(y)

    return theta

p = normalEquation(X, y)

print(p)

Using the small data set found here:

http://www.lauradhamilton.com/tutorial-linear-regression-with-octave

I get the co-efficients: [-0.34390603; 0.2124426 ] using the above code instead of: [24.9660; 3.3058]. Could anyone help clarify where I am going wrong?

you have your X and Y around the wrong way from the example! If I reverse them, I get the answers you suggest — jeremycg, Oct 05 '17 at 13:24

score 3 · Answer 1 · answered Oct 05 '17 at 16:19

3

You can implement normal equation like below:

import numpy as np

X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

X_b = np.c_[np.ones((100, 1)), X]  # add x0 = 1 to each instance
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)

X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]  # add x0 = 1 to each instance
y_predict = X_new_b.dot(theta_best)
y_predict

answered Oct 05 '17 at 16:19

Batuhan B

1,835
4
29
39

You can also use [`add_constant`](http://www.statsmodels.org/dev/generated/statsmodels.tools.tools.add_constant.html) for this. – Brad Solomon Oct 06 '17 at 13:14
1

Could you please explain why you add x0 = 1 to each instance? – Dory Nguyen Nov 27 '22 at 13:36

score 2 · Answer 2 · answered Feb 01 '19 at 01:18

This assumes X is an m by n+1 dimensional matrix where x_0 always = 1 and y is a m-dimensional vector.

import numpy as np

step1 = np.dot(X.T, X)
step2 = np.linalg.pinv(step1)
step3 = np.dot(step2, X.T)
theta = np.dot(step3, y) # if y is m x 1.  If 1xm, then use y.T

score 1 · Accepted Answer · answered Oct 05 '17 at 13:29

1

Your implementation is correct. You've only swapped X and y (look closely how they define x and y), that's why you get a different result.

The call normalEquation(y, X) gives [ 24.96601443 3.30576144] as it should.

answered Oct 05 '17 at 13:29

Maxim

52,561
27
155
209

Oh, the shame! Thank you both for your responses. – PS94 Oct 05 '17 at 13:50
1

Thanks for the advice @Maxim – PS94 Oct 07 '17 at 14:34

score 0 · Answer 4 · edited May 12 '20 at 12:27

0

Here is the normal equation in one line:

theta = np.dot(np.linalg.inv(np.dot(X.T,X)),np.dot(X.T,Y))

edited May 12 '20 at 12:27

FreshD

2,914
2
23
34

answered Apr 10 '19 at 19:23

SyedHashimReja

29
2

Normal Equation Implementation in Python / Numpy

4 Answers4