Normal Equation for linear regression

3
votes

I have the following X and y matrices:

enter image description here

for which I want to calculate the best value for theta for a linear regression equation using the normal equation approach with:

theta = inv(X^T * X) * X^T * y

the results for theta should be : [188.400,0.3866,-56.128,-92.967,-3.737]

I implement the steps with:

X=np.matrix([[1,1,1,1],[2104,1416,1534,852],[5,3,3,2],[1,2,2,1],[45,41,30,36]])
y=np.matrix([460,232,315,178])

XT=np.transpose(X)

XTX=XT.dot(X)

inv=np.linalg.inv(XTX)

inv_XT=inv.dot(XT)

theta=inv_XT.dot(y)

print(theta)

But I dont't get the desired results. Instead it throws an error with:

Traceback (most recent call last): File "C:/", line 19, in theta=inv_XT.dot(y) ValueError: shapes (4,5) and (1,4) not aligned: 5 (dim 1) != 1 (dim 0)

What am I doing wrong?

2
pythonnumpymachine-learninglinear-regression
I'd advise against using np.inv as that explicitly computes the inverse of a matrix. This is slower numerically less stable than simply solving the linear system and computing (X'X)^{-1}X' directly. To do this, write np.linalg.solve(XTX, X) (after making sure that XTX actually is X'X and not XX' as MaxU mentioned that you computed). – Yngve Moe

2 Answers

3
votes

I think you have messed up dimensions a little bit. Your X is actually a XT and XT is X.

Try this:

In [163]: X=np.matrix([[1,1,1,1],[2104,1416,1534,852],[5,3,3,2],[1,2,2,1],[45,41,30,36]]).T

In [164]: y=np.matrix([460,232,315,178])

In [165]: X
Out[165]:
matrix([[   1, 2104,    5,    1,   45],
        [   1, 1416,    3,    2,   41],
        [   1, 1534,    3,    2,   30],
        [   1,  852,    2,    1,   36]])

In [166]: XT = X.T

In [167]: np.linalg.inv(XT @ X) @ XT @ y.T
Out[167]:
matrix([[243.4453125 ],
        [ -0.47787476],
        [268.609375  ],
        [  3.1328125 ],
        [ -5.83056641]])

UPDATE: this approach gives values that are closer to your desired values:

In [197]: (np.linalg.inv(X @ X.T) @ X).T @ y.T
Out[197]:
matrix([[182.27200269],
        [  0.34497234],
        [-38.43393186],
        [-82.90625955],
        [ -3.84484213]])

UPDATE2: how to create a correct matrix initially:

In [217]: np.array([[1, 2104, 5, 1, 45],
     ...:  [1, 1416, 3, 2, 41],
     ...:  [1, 1534, 3, 2, 30],
     ...:  [1, 852, 2, 1, 36]])
     ...:
Out[217]:
array([[   1, 2104,    5,    1,   45],
       [   1, 1416,    3,    2,   41],
       [   1, 1534,    3,    2,   30],
       [   1,  852,    2,    1,   36]])
0
votes

I have solved the problem by using numpy.linalg.pinv() that is the "pseudo-inverse" instead of numpy.linalg.inv() for the inverting of the matrix since the ducumentation sais:

"The pseudo-inverse of a matrix A, denoted A^+, is defined as: “the matrix that ‘solves’ [the least-squares problem] Ax = b,” i.e., if \bar{x} is said solution, then A^+ is that matrix such that \bar{x} = A^+b."

and solving the least-squares problem is exactly what I want to achieve in the context of linear regression.

Consequently the code is:

X=np.matrix([[1,2104,5,1,45],[1,1416,3,2,40],[1,1534,3,2,30],[1,852,2,1,36]])
y=np.matrix([[460],[232],[315],[178]])

XT=X.T
XTX=XT@X

inv=np.linalg.pinv(XTX)

theta=(inv@XT)@y
print(theta)

[[188.40031946]
 [  0.3866255 ]
 [-56.13824955]
 [-92.9672536 ]
 [ -3.73781915]]

Edit: There is also the possibility of regularization to get rid of the problem of NoN-invertibility by changing the normal-equation to:

theta = (XT@X + lambda*matrix)^(-1)@XT@y where lambda is a real number and called regularization parameter and matrix is a (n+1 x n+1) dimensional matrix of the shape:

 0 0 0 0 ... 0 0 
 0 1 0 0 ... 0 0 
 0 0 1 0 ... 0 0
 0 0 0 1 ... 0 0
 .
 .
 .
 0 0 0 0 0 0 0 1

That is a eye() matrix with the element [0,0] set to 0

More about the concept of regularization can be read here