42
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ @hxd, barring any special structure in the design matrix, these are all $O(mn^2)$ algorithms, differing only in the constant factor. The decompositional approach is a good habit inherited from the tradition of numerical linear algebra. $\endgroup$ – J. M. is not a statistician Jun 15 '16 at 18:35
  • $\begingroup$ @hxd, and that's why my answer was tailored to be an exposition of the algorithms involved. If you have questions not covered by this thread, consider asking a new question. $\endgroup$ – J. M. is not a statistician Jun 15 '16 at 18:40
32
$\begingroup$

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 = T)

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])
$\endgroup$
  • $\begingroup$ I get an err could not find inv?! $\endgroup$ – hhh Dec 3 '11 at 9:16
  • 1
    $\begingroup$ Load the fBasics package. finzi.psych.upenn.edu/R/library/fBasics/html/matrix-inv.html $\endgroup$ – George Dontas Dec 3 '11 at 13:28
  • 5
    $\begingroup$ Is there a reason for using inv from fBasics when it's just a synonym for solve? Wouldn't it be better for the answer to not require a dependency on external packages if it is unnecessary? $\endgroup$ – Dason Dec 4 '11 at 20:46
  • $\begingroup$ @George I love the clear answer, However, I think the original question was asking algorithms, and QR is just one of it. How about LU, SVD, Cholesky decomposition? Also, in R, the method for lm is QR, there are reasons for that, could you explain why? $\endgroup$ – Haitao Du Jun 16 '16 at 15:26
  • $\begingroup$ @GeorgeDontas Note that it might be that $A^T A$ isn't invertible. As explained in this answer, one way to deal with that is removing from $A$ columns that are linear combinations of other columns. $\endgroup$ – Oren Milman Aug 25 '18 at 7:13
70
+50
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ Nice explanation! $\endgroup$ – Mike Spivey Oct 26 '11 at 4:23
  • $\begingroup$ How do you compute R^{-1} Q^T y if A is not square? Do you drop the zero rows in R? $\endgroup$ – bhan Apr 14 '16 at 18:32
  • 1
    $\begingroup$ @bhan, I am assuming the "economy" (or "thin") variant of QR decomposition, where $\mathbf R$ is square and $\mathbf Q$ has the same dimensions as the design matrix. Something for you to do: look up the difference between "full QR" and "thin QR". $\endgroup$ – J. M. is not a statistician Apr 14 '16 at 18:42
  • $\begingroup$ @J.M. any recommendations on " good book on computational statistics and numerical linear algebra"? really want to learn more. $\endgroup$ – Haitao Du Jun 15 '16 at 19:06
  • 1
    $\begingroup$ @hxd, off the top of my head: Monahan for computational statistics, and Golub/van Loan for numerical linear algebra. $\endgroup$ – J. M. is not a statistician Jun 15 '16 at 23:03
6
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Doesn't address the question (the page doesn't even mention QR) $\endgroup$ – user603 Dec 18 '12 at 13:46
4
$\begingroup$

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 (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.

$\endgroup$
  • $\begingroup$ Right, that's why I consider OLS to be an "algorithm" used in linear regression... $\endgroup$ – Belmont Aug 18 '10 at 15:22
  • 3
    $\begingroup$ Ordinary least squares is an estimator. There are a variety of algorithms for computing the estimate: usually some sort of orthogonal matrix decomposition, such as QR, is used. See en.wikipedia.org/wiki/… $\endgroup$ – Simon Byrne Aug 18 '10 at 15:34
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ By algorithm, I mean roughly the software implementation used for fitting a line to model the mean of a distribution. $\endgroup$ – Belmont Aug 18 '10 at 15:18
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.