# Why does ordinary least squares have to be linear in the parameters?

I've been looking into linear regression, and on the wikipedia page it says:

"In contrast, non-linear least squares problems generally must be solved by an iterative procedure"

This got me thinking more about OLS, and the differences between it and non-linear regression methods. Mores specifically, why equations that are non-linear in their parameters can't also be solved using the OLS assumption that $$y=\beta x$$ where $$\beta =(X^TX)^{-1}X^Ty$$.

So i guess my question is:

What is it about the process of solving OLS that requires the parameters to be linear? What would happen if they were non-linear and we tried to solve using OLS?

• You can minimize $L(\beta)$ by solving $\nabla L(\beta) = 0$. Usually this is a nonlinear system of equations, which makes it difficult to solve. But if $L(\beta) = (1/2) \| X \beta - y \|^2$, then $\nabla L(\beta) = X^T (X \beta - y)$, so $\nabla L(\beta) = 0$ is equivalent to $X^T X \beta = X^T y$. This is a linear system of equations that can be solved with Gaussian elimination! That is the special thing that makes linear regression easy. – littleO Jun 1 at 5:48

What is it about the process of solving OLS that requires the parameters to be linear?

Because equations which are nonlinear in their parameters can't be written as $$y=X\beta$$. OLS estimates $$\beta$$ in the equation

$$y = X\beta +\epsilon.$$

This is a linear relationship, so when we say that $$\hat{\beta} = (X^\top X)^{-1}X^\top y$$ is the optimal estimator of $$\beta$$, what we mean is that it's optimal in the sense that it minimizes $$\|y - X\beta\|_2^2$$. Minimizing $$\|y - X\beta\|_2^2$$ is only important if this objective is meaningful for your task; particularly, if the task isn't linear in these parameters, then the fit may be poor.

However, one reason that OLS is so flexible is that if you can find a way to represent your data in a linear way, then it is linear in the parameters, otherwise known as basis expansion.

A textbook example of a change of basis is using a polynomial basis, so you have $$X_\text{polynomial} = [1, x, x^2, x^3, \dots, x^p]$$. The model $$X_\text{polynomial}\beta$$ is linear in its parameters, but viewed as a function of $$x$$, it's a nonlinear polynomial.

What would happen if they were non-linear and we tried to solve using OLS?

It won't work very well!

This data's deterministic component is given by $$y = \beta_0 + \beta_1 \sin (\beta_2 x + \beta_3)$$ which is not linear in $$\beta$$, the parameter vector to be estimated, because you can't write this in the form $$y=X\beta$$. I also add small, independent 0-mean Gaussian noise to each observation.

If we do the naive thing and assume that our output $$y$$ is a linear function of $$x$$, then we find a poor fit, in the sense that there is a large discrepancy between the estimated line (red) and the true function (blue). The model finds that the best linear approximation is a decreasing line, completely ignoring the sinusoidal behavior.

One way to try to improve the fit is to re-express $$x$$. Since this looks like something sinusoidal, we might try a sine function. This gives the design matrix $$X_\text{sine}=[1, \sin(x)]$$. This give a flatter line, but it's still not a satisfying model. Even though the model and the desired function are both sine waves, we're implicitly using $$\beta_0 + \beta_1 \sin(1 \times x + 0)$$ to approximate $$y = \beta_0 + \beta_1 \sin (\beta_2 x + \beta_3).$$ This is not a good approximation, because we've fixed $$\beta_2=1$$ and $$\beta_3=0$$, so the further the true values are from these assumed values, the poorer this approximation will be.

What we really need is a way to recover all of the parameters in the function $$y = \beta_0 + \beta_1 \sin (\beta_2 x + \beta_3),$$ but this is a nonlinear estimation task, so we need to use the appropriate tools to accommodate the nonlinearity of the $$\beta$$s. Nonlinear least squares is one method to achieve this, among many others.

# Code

set.seed(13)
N <- 1000

x <- runif(N, -pi, pi)
f <- function(x) pi + 2 * sin(4 * x)
y <- f(x) + rnorm(N,sd=0.5)

model <- lm(y ~ x)

png("~/Desktop/nonlinear.png")
plot(x,y,col="grey")
abline(model, col="red", lwd=2, lty="dashed")
lines(sort(x), f(sort(x)), lwd=2, col="blue")
dev.off()

model2 <- lm(y ~ sin(x) )

png("~/Desktop/nonlinear2.png")
plot(x,y,col="grey")
abline(model2, col="red", lwd=2, lty="dashed")
lines(sort(x), f(sort(x)), lwd=2, col="blue")
dev.off()

• In this particular case, I’d change it to b1 sin (b3x) + b2 cos (b3 x), which gives exactly the same numbers but has one Non-linear parameter only. – gnasher729 Jun 1 at 14:29
• That's a useful observation for a person carrying out NLLS on this toy problem, but it doesn't do much to explain why OLS is a poor substitute for NLLS. – Sycorax Jun 1 at 17:25