Why additional regressors in OLS increase SSR?

Question

I have two nested OLSes:

>>> f1 = 'Y ~ np.log(X + 1) + X + np.log(X) + np.exp(-X) + np.sqrt(X) + np.power(X, 4)'
>>> f2 = 'Y ~ np.log(X + 1) + X + np.log(X) + np.exp(-X) + np.sqrt(X) + np.power(X, 4) + np.power(X, -1) + np.power(X, -2) + np.power(X, 2) + np.power(X, 3) + np.exp(X) + np.power(X, -3)'

When analysing this dataset (17M) with statsmodels I have discovered that the bigger model has larger sum of squared residuals (SSR):

>>> ssr1 = sm.ols(formula=f1, data=DF).fit().ssr
>>> ssr2 = sm.ols(formula=f2, data=DF).fit().ssr
>>> ssr1 - ssr2
-1037.4640076523647

Same happens with sklearn.

As OLS is expected to minimize SSR, I am surprised that inclusion of extra regressors increased it. As regressors of f2 are a superset of those of f1, ssr2 may be at least as small as ssr1 by simply setting coefficients of all extra parameters to 0.

So why OLS fails to minimize SSR? Is it some numerical artefact which may be ignored or something more serious?

VIF for f1:

np.log(X + 1)   13567.4057637
X               11541.1910286
np.log(X)       446.977604585
np.exp(-X)      38.674353333
np.sqrt(X)      53077.4595776
np.power(X, 4)  34.2445023554

VIF for f2:

np.log(X + 1)   28942.1737842
X               180304.260711
np.log(X)       121.384763923
np.exp(-X)      22.2770600233
np.sqrt(X)      319544.359507
np.power(X, 4)  4606.23688654
np.power(X, -1) 113.759877386
np.power(X, -2) 2452.98781276
np.power(X, 2)  19895.5099527
np.power(X, 3)  24609.1808767
np.exp(X)       3.45225568589
np.power(X, -3) 1605.22900297

Answer 1

The direct problem is that your data matrix has a ridiculously high condition number, and we're probably not even finding the solution to the least squares problem. Up to numerical precision, the data matrix is rank deficient, you're not successfully solving the least squares problems, and the properties of a least squares solution don't apply.

As @Josef in the comment suggests, you have an ill-posed problem because of massively unequal scaling. Standardizing the variables improves matters quite a lot.

Even with better scaling, it's quite problematic to run a kitchen sink regression with neither clear motivation for the right hand side variables nor regularlization to reduce overfitting. I'd have serious concerns as to model performance on new, independent test data.

Discussion

Let $X$ be an $n$ by $k$ matrix representing your data. A classic assumption to run OLS is no collinearity, that $X$ is rank $k$. If your data lies in a $s < k$ dimensional subspace, things are going to go haywire without regularization.

What does this mean once we move to computers and actual, numerical computing?

A quick diagnostic is to examine the condition number. Wolfram's Mathworld says

The base-b logarithm of C is an estimate of how many base-b digits are lost in solving a linear system with that matrix. In other words, it estimates worst-case loss of precision.

With

X2 = np.matrix([np.log(DF.X + 1), DF.X, np.log(DF.X), np.exp(-DF.X), np.sqrt(DF.X), np.power(DF.X, 4), np.power(DF.X, -1), np.power(DF.X, -2), np.power(DF.X, 2), np.power(DF.X, 3), np.exp(DF.X), np.power(DF.X, -3)])
X2_condition_number = np.linalg.cond(X2)

You find that the condition number is about 7.7e+15, which means you can lose 15-16 digits of precision when you solve a linear system with X2. Double precision floating point only has 15-17 digits of precision! The first digit of the solution returned by numpy could be wrong.

The solution isn't to use more advanced numerical techniques, it's that numpy is being asked to solve a ridiculous, collinear least squares problem. Such systems are ridiculously sensitive to small changes.

Example of a system with a ridiculous condition number

The matrix on the left has a condition number of 2e+15.

$$ \begin{bmatrix} 1 & 10000000.0 \\ 1 & 10000000.1 \end{bmatrix} \begin{bmatrix}b_1 \\ b_2 \end{bmatrix} = \begin{bmatrix} 2 \\ 1\end{bmatrix}$$

The first column and the 2nd column are almost collinear.

Ways forward

• Standardize variables.

As @Josef described in the comments, you can have ill-conditioned problems because of widely unequal scaling.

X3 = (X2 - X2.mean(axis=1)) / X2.std(axis=1)

X3 now has a condition number of 2.9e+3. Another common scaling in machine learning is max-min scaling.

Other ways forward include:

• Drop regressors that are collinear with existing regressors (i.e. stop doing a kitchen sink regression).
• Add regularlization, eg. do ridge regression or LASSO.

If you're trying to fit a flexible curve, there are better things to do such as low order polynomials with splines.

Answer 2

I would like to make a specific point about polynomials where high condition number is a side-effect for improper scaling or choosing badly scaled basis functions.

One of the NIST cases to check the numerical reliability of econometric software packages uses a badly scaled polynomial of order 10. The smallest eigenvalue is around zero and makes linear regression numerically unstable. http://jpktd.blogspot.com/2012/03/numerical-accuracy-in-linear-least.html

However, when we estimate the polynomial with numpy polynomial fitting, then it does not have any problems at all because it rescales the values for the polynomial before computing the polynomial basis functions, i.e. uses a fixed domain which is by default [-1, 1].

Similarly using orthogonal polynomials like Chebychev (on fixed domain) does not have the same problem as standard polynomials.

Related: Spline basis: I was running an example with a very large number of truncated power splines. By construction those are highly collinear, and computing standard errors, i.e. inverting the moment matrix failed in most cases. Using B-splines with the same number of knots has no multicollinearity problem or problem of high condition number by the way the basis functions are defined.

Aside: statsmodels and scikit-learn use SVD for estimating linear regression parameters. This removes singular vectors with small singular values. However, the singular value decomposition will put too much emphasis on variables with large fluctuations or magnitudes and so cannot really correct for badly scaled columns in the design matrix.