What algorithm is used in linear regression?

Question

I usually hear about "ordinary least squares". Is that the most widely used algorithm used for linear regression? Are there reasons to use a different one?

Answer 1

Regarding the question in the title, about what is the algorithm that is used:

In a linear algebra perspective, the linear regression algorithm is the way to solve a linear system $\mathbf{A}x=b$ with more equations than unknowns. In most of the cases there is no solution to this problem. And this is because the vector $b$ doesn't belong to the column space of $\mathbf{A}$, $C(\mathbf{A})$.

The best straight line is the one that makes the overall error $e=\mathbf{A}x-b$ as small as it takes. And is convenient to think as small to be the squared length, $\lVert e \rVert^2$, because it's non negative, and it equals 0 only when $b\in C(\mathbf{A})$.

Projecting (orthogonally) the vector $b$ to the nearest point in the column space of $\mathbf{A}$ gives the vector $b^*$ that solves the system (it's components lie on the best straight line) with the minimum error.

$\mathbf{A}^T\mathbf{A}\hat{x}=\mathbf{A}^Tb \Rightarrow \hat{x}=(\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^Tb$

and the projected vector $b^*$ is given by: $$ b^*=\mathbf{A}\hat{x}=\mathbf{A}(\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^Tb $$

Perhaps the least squares method is not exclusively used because that squaring overcompensates for outliers.

Let me give a simple example in R, that solves the regression problem using this algorithm:

library(fBasics)

reg.data <- read.table(textConnection("
   b      x
  12      0
  10      1
   8      2
  11      3
   6      4
   7      5
   2      6
   3      7
   3      8 "), header = TRUE)

attach(reg.data)

A <- model.matrix(b ~ x)

# intercept and slope
inv(t(A) %*% A) %*% t(A) %*% b

# fitted values - the projected vector b in the C(A)
A %*% inv(t(A) %*%A ) %*% t(A) %*% b

# The projection is easier if the orthogonal matrix Q is used, 
# because t(Q)%*%Q = I
Q <- qr.Q(qr(A))
R <- qr.R(qr(A))

# intercept and slope 
best.line <- inv(R) %*% t(Q) %*% b

# fitted values 
Q %*% t(Q) %*% b

plot(x,b,pch=16)
abline(best.line[1],best.line[2])

Answer 2

To answer the letter of the question, "ordinary least squares" is not an algorithm; rather it is a type of problem in computational linear algebra, of which linear regression is one example. Usually one has data $\{(x_1,y_1),\dots,(x_m,y_m)\}$ and a tentative function ("model") to fit the data against, of the form $f(x)=c_1 f_1(x)+\dots+c_n f_n(x)$. The $f_j(x)$ are called "basis functions" and can be anything from monomials $x^j$ to trigonometric functions (e.g. $\sin(jx)$, $\cos(jx)$) and exponential functions ($\exp(-jx)$). The term "linear" in "linear regression" here does not refer to the basis functions, but to the coefficients $c_j$, in that taking the partial derivative of the model with respect to any of the $c_j$ gives you the factor multiplying $c_j$; that is, $f_j(x)$.

One now has an $m\times n$ rectangular matrix $\mathbf A$ ("design matrix") that (usually) has more rows than columns, and each entry is of the form $f_j(x_i)$, $i$ being the row index and $j$ being the column index. OLS is now the task of finding the vector $\mathbf c=(c_1\,\dots\,c_n)^\top$ that minimizes the quantity $\sqrt{\sum\limits_{j=1}^{m}\left(y_j-f(x_j)\right)^2}$ (in matrix notation, $\|\mathbf{A}\mathbf{c}-\mathbf{y}\|_2$ ; here, $\mathbf{y}=(y_1\,\dots\,y_m)^\top$ is usually called the "response vector").

There are at least three methods used in practice for computing least-squares solutions: the normal equations, QR decomposition, and singular value decomposition. In brief, they are ways to transform the matrix $\mathbf{A}$ into a product of matrices that are easily manipulated to solve for the vector $\mathbf{c}$.

George already showed the method of normal equations in his answer; one just solves the $n\times n$ set of linear equations

$\mathbf{A}^\top\mathbf{A}\mathbf{c}=\mathbf{A}^\top\mathbf{y}$

for $\mathbf{c}$. Due to the fact that the matrix $\mathbf{A}^\top\mathbf{A}$ is symmetric positive (semi)definite, the usual method used for this is Cholesky decomposition, which factors $\mathbf{A}^\top\mathbf{A}$ into the form $\mathbf{G}\mathbf{G}^\top$, with $\mathbf{G}$ a lower triangular matrix. The problem with this approach, despite the advantage of being able to compress the $m\times n$ design matrix into a (usually) much smaller $n\times n$ matrix, is that this operation is prone to loss of significant figures (this has something to do with the "condition number" of the design matrix).

A slightly better way is QR decomposition, which directly works with the design matrix. It factors $\mathbf{A}$ as $\mathbf{A}=\mathbf{Q}\mathbf{R}$, where $\mathbf{Q}$ is an orthogonal matrix (multiplying such a matrix with its transpose gives an identity matrix) and $\mathbf{R}$ is upper triangular. $\mathbf{c}$ is subsequently computed as $\mathbf{R}^{-1}\mathbf{Q}^\top\mathbf{y}$. For reasons I won't get into (just see any decent numerical linear algebra text, like this one), this has better numerical properties than the method of normal equations.

One variation in using the QR decomposition is the method of seminormal equations. Briefly, if one has the decomposition $\mathbf{A}=\mathbf{Q}\mathbf{R}$, the linear system to be solved takes the form

$$\mathbf{R}^\top\mathbf{R}\mathbf{c}=\mathbf{A}^\top\mathbf{y}$$

Effectively, one is using the QR decomposition to form the Cholesky triangle of $\mathbf{A}^\top\mathbf{A}$ in this approach. This is useful for the case where $\mathbf{A}$ is sparse, and the explicit storage and/or formation of $\mathbf{Q}$ (or a factored version of it) is unwanted or impractical.

Finally, the most expensive, yet safest, way of solving OLS is the singular value decomposition (SVD). This time, $\mathbf{A}$ is factored as $\mathbf{A}=\mathbf{U}\mathbf \Sigma\mathbf{V}^\top$, where $\mathbf{U}$ and $\mathbf{V}$ are both orthogonal, and $\mathbf{\Sigma}$ is a diagonal matrix, whose diagonal entries are termed "singular values". The power of this decomposition lies in the diagnostic ability granted to you by the singular values, in that if one sees one or more tiny singular values, then it is likely that you have chosen a not entirely independent basis set, thus necessitating a reformulation of your model. (The "condition number" mentioned earlier is in fact related to the ratio of the largest singular value to the smallest one; the ratio of course becomes huge (and the matrix is thus ill-conditioned) if the smallest singular value is "tiny".)

This is merely a sketch of these three algorithms; any good book on computational statistics and numerical linear algebra should be able to give you more relevant details.

Answer 3

The wiki link: Estimation Methods for Linear Regression gives a fairly comprehensive list of estimation methods including OLS and the contexts in which alternative estimation methods are used.

Answer 4

It is easy to get confused between definitions and terminology. Both terms are used, sometimes interchangeably. A quick lookup on Wikipedia should help:

• ordinary least squares
• lnear regression

Ordinary Least Squares (OLS) is a method used to fit linear regression models. Because of the demonstrable consistency and efficiency (under supplementary assumptions) of the OLS method, it is the dominant approach. See the articles for further leads.

Answer 5

I tend to think of 'least squares' as a criterion for defining the best fitting regression line (i.e., that which makes the sum of 'squared' residuals 'least') and the 'algorithm' in this context as the set of steps used to determine the regression coefficients that satisfy that criterion. This distinction suggests that it is possible to have different algorithms that would satisfy the same criterion.

I'd be curious to know whether others make this distinction and what terminology they use.

Answer 6

An old book, yet one I find myself repeatedly turning to, is

Lawson, C.L. and Hanson, R.J. Solving Least Squares Problems, Prentice-Hall, 1974.

It contains a detailed and very readable discussion of some of the algorithms that previous answers have mentioned. You might want to look at it.