Linear Regression

Linear regression is a method for modeling the relationship between two scalar values: the input variable x and the output variable y. The model assumes that $y$ is a linear function or a weighted sum of the input variable. $y = f(x)$
Or, stated with the coefficients.
$y = b_0 + b_1 x_1$

The model can also be used to model an output variable given multiple input variables called multivariate linear regression $y = b_0 + (b_1.x_1) + (b_2.x_2) + ....+(b_n.x_n)$

The objective of creating a linear regression model is to find the values for the coefficient values ($b$) that minimize the error in the prediction of the output variable $y$.

 
Matrix Formulation of Linear Regression

 
Linear regression can be stated using Matrix notation; for example:
$y = X.b$

 
Where $X$ is the input data and each column is a data feature, $b$ is a vector of coefficients and $y$ is a vector of output variables for each row in $X$.

 

Reformulated, the problem becomes a system of linear equations where the $b$ vector values are unknown. This type of system is referred to as over determined because there are more equations than there are unknowns, i.e. each coefficient is used on each row of data. It is a challenging problem to solve analytically because there are multiple inconsistent solutions, e.g.multiple possible values for the coefficients. Further, all solutions will have some error because there is no line that will pass nearly through all points, therefore the approach to solving the equations must be able to handle that. The way this is typically achieved is by finding a solution where the values for $b$ in the model minimize the squared error. This is called linear least squares.

In matrix notation, this problem is formulated using the so-named normal equation:

$X^T.y=X^T.X.b$

This can be re-arranged in order to specify the solution for $b$ as:

$b= (X^T.X)^{-1}.X^T.y$

 
This can be solved directly, although given the presence of the matrix inverse can be numerically challenging or unstable. 

Solve via Inverse

The first approach is to attempt to solve the regression problem directly using the matrix inverse.

That is, given $X$, what are the set of coefficients $b$ that when multiplied by $X$ will give $y$. As we saw in a previous section, the normal equations define how to calculate $b$ directly.

$b= (X^T.X)^{-1}.X^T.y$

This can be calculated directly in NumPy using the inv() function for calculating the matrix inverse.

 
$b = inv(X^T.dot(X)).dot(X^T).dot(y)$

 
Once the coefficients are calculated, we can use them to predict outcomes given $X$.
$\hat{y} = X.dot(b)$

Following is the example code
**********************************************
# direct solution to linear least squares
from numpy import array
from numpy.linalg import inv
from matplotlib import pyplot
# define dataset
data = array([
[0.05, 0.12],
[0.18, 0.22],
[0.31, 0.35],
[0.42, 0.38],
[0.5, 0.49]])
# split into inputs and outputs
X, y = data[:,0], data[:,1]
X = X.reshape((len(X), 1))
# linear least squares
b = inv(X.T.dot(X)).dot(X.T).dot(y)
print("coefficient obtained")
print(b)
# predict using coefficients
yhat = X.dot(b)
print("Predictions")
print(yhat)
print("actual values")
print(y)
# plot data and predictions
pyplot.scatter(X, y)
pyplot.plot(X, yhat, color='red')
pyplot.show()
***************************************
coefficient obtained
[1.00233226]
Predictions
[0.05011661 0.18041981 0.310723 0.42097955 0.50116613]
actual values
[0.12 0.22 0.35 0.38 0.49]


Solve via QR Decomposition
The $QR$ decomposition is an approach of breaking a matrix down into its constituent elements.
$A = Q . R$
Where $A$ is the matrix that we wish to decompose, $Q$ is a orthogonal matrix with the size $m \times m$, and $R$ is an upper triangle matrix with the size $m \times n$. The $QR$ decomposition is a popular approach for solving the linear least squares equation. The coefficients can be found using the $Q$ and $R$ elements as follows:
$y = X.b$
$y=Q.R.b$
$b=(QR)^{-1} y$
$b=R^{-1}Q^T.y$

The approach still involves a matrix inversion, but in this case only on the simpler $R$ matrix.The $QR$ decomposition can be found using the qr() function in NumPy. The follwing program demonstrate this.

*************************************************************

# QR decomposition solution to linear least squares
from numpy import array
from numpy.linalg import inv
from numpy.linalg import qr
from matplotlib import pyplot
# define dataset
data = array([
[0.05, 0.12],
[0.18, 0.22],
[0.31, 0.35],
[0.42, 0.38],
[0.5, 0.49]])
# split into inputs and outputs
X, y = data[:,0], data[:,1]
X = X.reshape((len(X), 1))
# linear least squares using qr decomposition
q,r=qr(X)
b = inv(r).dot(q.T).dot(y)
print("coefficent obtained")
print(b)
# predict using coefficients
yhat = X.dot(b)
print("Predictions")
print(yhat)
print("actual values")
print(y)
# plot data and predictions
pyplot.scatter(X, y)
pyplot.plot(X, yhat, color='red')
pyplot.show()
**************************************************
coefficient obtained
[1.00233226]
Predictions
[0.05011661 0.18041981 0.310723 0.42097955 0.50116613]
actual values
[0.12 0.22 0.35 0.38 0.49]

**************************************************
Solve via SVD and Pseudoinverse

The pseudo inverse is the generalization of the matrix inverse for square matrices to rectangular matrices where the number of rows and columns are not equal. It is also called the Moore-Penrose Inverse after two independent discoverers of the method or the Generalized Inverse.Matrix inversion is not defined for matrices that are not square. When A has more columns than rows, then solving a linear equation using the pseudo inverse provides one of the many possible solutions.

The pseudo inverse is of matrix A is calculated using the singular value decomposition of A:

The pseudo inverse of 

$A = V . D .U^T$

where $V$ and $U$ is the matrix obtained from SVD and $D$ is the inverse of the diagonal matrix $Ʃ.D$ can be obtained by calculating the reciprocal of each non zero element in Ʃ and taking the transpose if the original matrix was rectangular.

The pseudo inverse provides one way of solving the linear regression equation, specifically when there are more rows than there are columns, which is often the case.NumPy provides the function pinv() for calculating the pseudo inverse of a rectangular matrix.
*********************************************************
# SVD solution via pseudo inverse to linear least squares
from numpy import array
from numpy.linalg import inv
from numpy.linalg import pinv
from matplotlib import pyplot
# define dataset
data = array([
[0.05, 0.12],
[0.18, 0.22],
[0.31, 0.35],
[0.42, 0.38],
[0.5, 0.49]])
# split into inputs and outputs
X, y = data[:,0], data[:,1]
X = X.reshape((len(X), 1))
# linear least squares using qr SVD pseudo inverse
b = pinv(X).dot(y)
print("coefficient obtained")
print(b)
# predict using coefficients
yhat = X.dot(b)
print("Predictions")
print(yhat)
print("actual values")
print(y)
# plot data and predictions
pyplot.scatter(X, y)
pyplot.plot(X, yhat, color='red')
pyplot.show()


*****************************************************

Solve via Convenience Function

The pseudo inverse via SVD approach to solving linear least squares is the defacto standard.This is because it is stable and works with most datasets. NumPy provides a convenience function named lstsq() that solves the linear least squares function using the SVD approach.

The function takes as input the $X$ matrix and $y$ vector and returns the $b$ coefficients as well as residual errors, the rank of the provided $X$ matrix and the singular values. The example below demonstrate the lstsq() function on the test dataset.

# least squares via convenience function
from numpy import array
from numpy.linalg import lstsq
from matplotlib import pyplot
# define dataset
data = array([
[0.05, 0.12],
[0.18, 0.22],
[0.31, 0.35],
[0.42, 0.38],
[0.5, 0.49]])
# split into inputs and outputs
X, y = data[:,0], data[:,1]
X = X.reshape((len(X), 1))
# linear least squares using qr SVD pseudo inverse
b,residuals,rank,s = lstsq(X,y,rcond=None)
print("coefficent obtained")
print(b)
# predict using coefficients
yhat = X.dot(b)
print("Predictions")
print(yhat)
print("actual values")
print(y)
# plot data and predictions
pyplot.scatter(X, y)
pyplot.plot(X, yhat, color='red')
pyplot.show()