17
$\begingroup$

I saw this list here and couldn't believe there were so many ways to solve least squares. The "normal equations" on Wikipedia seemed to be a fairly straight forward way: $$ {\displaystyle {\begin{aligned}{\hat {\alpha }}&={\bar {y}}-{\hat {\beta }}\,{\bar {x}},\\{\hat {\beta }}&={\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}}\end{aligned}}} $$

So why not just use them? I assumed there must be a computational or precision issue given that in the first link above Mark L. Stone mentions that SVD or QR are popular methods in statistical software and that the normal equations are "TERRIBLE from reliability and numerical accuracy standpoint". However, in the following code, the normal equations are giving me accuracy to ~12 decimal places when compared to three popular python functions: numpy's polyfit; scipy's linregress; and scikit-learn's LinearRegression.

What's more interesting is that the normal equation method is the fastest when n = 100000000. Computational times for me are: 2.5s for linregress; 12.9s for polyfit; 4.2s for LinearRegression; and 1.8s for the normal equation.

Code:

import numpy as np
from sklearn.linear_model import LinearRegression
from scipy.stats import linregress
import timeit

b0 = 0
b1 = 1
n = 100000000
x = np.linspace(-5, 5, n)
np.random.seed(42)
e = np.random.randn(n)
y = b0 + b1*x + e

# scipy                                                                                                                                     
start = timeit.default_timer()
print(str.format('{0:.30f}', linregress(x, y)[0]))
stop = timeit.default_timer()
print(stop - start)

# numpy                                                                                                                                      
start = timeit.default_timer()
print(str.format('{0:.30f}', np.polyfit(x, y, 1)[0]))
stop = timeit.default_timer()
print(stop - start)

# sklearn                                                                                                                                    
clf = LinearRegression()
start = timeit.default_timer()
clf.fit(x.reshape(-1, 1), y.reshape(-1, 1))
stop = timeit.default_timer()
print(str.format('{0:.30f}', clf.coef_[0, 0]))
print(stop - start)

# normal equation                                                                                                                            
start = timeit.default_timer()
slope = np.sum((x-x.mean())*(y-y.mean()))/np.sum((x-x.mean())**2)
stop = timeit.default_timer()
print(str.format('{0:.30f}', slope))
print(stop - start) 
$\endgroup$
  • $\begingroup$ The answers are quite exaggerated. It's not so horrible if you just avoid explicitly calculating the inverse. $\endgroup$ – mathreadler Apr 27 '18 at 10:01
  • 3
    $\begingroup$ A few notes on speed: you are only looking at a single covariate, so the cost of matrix inversion is essentially 0. If you look at a few thousand covariates, that will change. Second, because you only have a single covariate, data-munging is what actually takes a lot of your time in the packaged competitors (but this should only scale linearly, so not a big deal). The normal equations solution does not do data munging, so it's faster, but has no bells and whistles attached with it's results. $\endgroup$ – Cliff AB Apr 27 '18 at 19:20
22
$\begingroup$

For the problem $Ax \approx b$, forming the Normal equations squares the condition number of $A$ by forming $A^TA$. Roughly speaking $log_{10}(cond)$ is the number of digits you lose in your calculation if everything is done well. And this doesn't really have anything to do with forming the inverse of $A^TA$. No matter how $A^TAx = A^Tb$ is solved, you've already lost $log_{10}(cond(A^TA)) = 2*log_{10}(cond(A))$ digits of accuracy. I.e., forming the Normal equations has doubled the number of digits of accuracy lost, right off the bat.

If the condition number is small (one is the best possible), it doesn't matter much. If the condition number = $10^8$ and you use a stable method such as QR or SVD, you may have about 8 digits of accuracy in double precision. If you form the Normal equations, you've squared the condition number to $10^{16}$, and you have essentially no accuracy in your answer.

Sometimes you get away with the Normal equations, and sometimes you don't.

$\endgroup$
  • 2
    $\begingroup$ A simpler way to see it (if you don't know/care about condition numbers) is that you're (essentially) multiplying something by itself ("squaring" it), meaning that you can expect to lose about half your bits of precision. (This should be more obvious if A is a scalar, and it should be easy to see that making A a matrix doesn't really change the underlying problem.) $\endgroup$ – Mehrdad Apr 27 '18 at 5:41
  • $\begingroup$ Besides the differences in accuracy, is there also a big speed difference between QR and normal equations? because in the latter case you might be solving (X'X)-1*X'Y, which is slow because of the inverse? I ask because Im not sure how QR works, so maybe theres something there that is just as slow as inverting a matrix. Or is the only point of consideration the accuracy loss? $\endgroup$ – Simon Apr 27 '18 at 5:43
  • 4
    $\begingroup$ @Simon Well, when you solve the normal equations you never actually form the inverse matrix, which is too slow. You do have to form the matrix $A^TA$ and the vector $A^Tb$ though, and then you solve the system. $\endgroup$ – Ant Apr 27 '18 at 8:05
8
$\begingroup$

If you only have to solve this one variable problem, then go ahead and use the formula. There's nothing wrong with it. I could see you writing a few lines of code in ASM for an embedded device, for instance. In fact, I used this kind of solution in some situations. You don't need to drag large statistical libraries just to solve this one little problem, of course.

The numerical instability and performance are issues of larger problems and general setting. If you solve multivariate least squares etc. For a general problem you wouldn't use this, of course.

$\endgroup$
0
$\begingroup$

No modern statistical package would solve a linear regression with the normal equations. The normal equations exist only in the statistical books.

The normal equations shouldn't be used as computing the inverse of matrix is very problematic.

Why use gradient descent for linear regression, when a closed-form math solution is available?

... although direct normal equation is available. Notice that in normal equation one has to invert a matrix. Now inverting a matrix costs O(N3) for computation where N is the number of rows in X matrix i.e. the observations. Moreover, if the X is ill conditioned then it will create computational errors in estimation...

$\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.