As part of a self-study exercise, I am comparing various implementations of polynomial regression:

• Closed form solution
• Scipy optimize
• Sklearn
• Statsmodel

When the problem involves polynomials of degree 3 or less, no problem, all three approaches yield the same coefficients. However, when the order increases to degree 5, 10 or even 15, I find it impossible to find the correct minimum using my numpy and scipy.optimize implementations.

### Question:

Why is gradient descent, and to a certain extent the scipy.optimize algorithm, so bad a optimizing polynomial regression ?

Is this because the cost function is non convex ? Not smooth ? Due to numerical instability or collinearity ?

### Example

In my model, there is only one variable and design matrix takes the form $1,x, x^2, x^3, ..., x^n$. The data is based on a sine function with uniform noise.

#Initializing noisy non linear data
x = np.linspace(0,1,40)
noise = 1*np.random.uniform(  size = 40)
y = np.sin(x * 1.5* np.pi )
y_noise = (y + noise-1).reshape(-1,1)


### Polynomial order 3

• Closed form solution: $(X^TX)^{-1}X^Ty = \begin{bmatrix} 0.07 & 10.14 & -20,15 & 9.1 \end{bmatrix}$
• Numpy gradient descent Same coefficients with 50,000 iterations and stepsize = 1
• Scipy optimize Same coefficients using BFGS method and the first derivative (gradient)
• Sklearn: same coefficients
• Statsmodel: same coefficients

### Polynomial order 5

• Closed form solution: $(X^TX)^{-1}X^Ty = \begin{bmatrix} 0.65 & 5.82 & -17.82 & 29.10 & -35.25 & 17.08 \end{bmatrix}$
• Numpy gradient descent Smaller coefficients with 50,000 iterations and stepsize = 1: $\begin{bmatrix} 0.71 & 3.98 & -5.2 & -3.23 & -0.08 & 3.44 \end{bmatrix}$
• Scipy optimize Also smaller coefficients, of the same order as with the Numpy implementation. Using BFGS method and the first derivative (gradient): $\begin{bmatrix} 0.70 & 4.14 & -5.83 & -2.73 & 0.18 & 3.09 \end{bmatrix}$
• Sklearn: same as analytical solution
• Statsmodel: same as analytical solution

### Polynomial order 16+

All methods give different results.

As the question is quite long already, you'll find the code here

• Gradient descent is bad at optimizing functions generally. In addition, raw polynomial regression is known to suffer from poor numerical conditioning.
– Sycorax
Commented Jun 6, 2018 at 19:20
• Possible duplicate of Do we need gradient descent to find the coefficients of a linear regression model
– Sycorax
Commented Jun 6, 2018 at 19:21
• High degree polynomials can have very large or very small values which can cause numerical problems. Your gradients could be quite large and so you're failing to converge. If you standardize each polynomial separately, gradient descent will probably do a better job. Lowering step size could help too. Commented Jun 6, 2018 at 19:27
• Polynomial fits are fine, my question was related to why gradient descent (and to some extent the scipy.optimize function I used) were struggling so much to find the optimal solution). I thought that linear regression was a convex optimization problem - I could rephrase the question indeed Commented Jun 6, 2018 at 19:27
• Just to point out: In you notebook statsmodels and scikit-learn have the same parameter estimate because both use a SVD based algorithm, which is numerically more stable than inverting the moment matrix in your "analytical" solution. My guess is that a "home-made" gradient descent accumulates more floating point noise than the "tuned" LAPACK routines. Commented Jun 6, 2018 at 19:42

Is this because the cost function is non convex ? Not smooth ? Due to numerical instability or collinearity ?

This appears to be simple linear regression with a sum-of-squares loss function. If you are able to obtain a closed form solution (i.e. $X^TX$ is invertible) then that loss function is both convex and continuously differentiable (smooth). (1, 2)

Why is gradient descent, and to a certain extent the scipy.optimize algorithm, so bad a optimizing polynomial regression ?

Gradient descent is known to be both slow (compared to second-derivative methods) and sensitive to step size. I also want to second what @Sycorax and @Jonny Lomond put in the comments - this particular problem is a difficult one for GD because of the massive magnitude difference across your dimensions, and your closed form solution may also be unstable. This link has has some really fantastic material on optimization challenges and momentum-based solutions including a polynomial regression example.

A few approaches you might consider:

1. As @Jonny Lomond suggested, standardize each polynomial separately, or tune your step size.
2. Plot your loss function over iterations to determine if there are any obvious problems with your optimization. If your gradient is "overshooting", you could try using an adaptive step size (reducing it as a function of the number of iterations).
3. Use backtracking to dynamically determine a better step size at each iteration.
4. Use a momentum-based gradient method like Nesterov accelerated gradient descent. These approaches are almost as fast (in terms of convergence) as second order methods in practice.
• One thing that I would add to this answer is that the matrix $X^\top X$ for several polynomial terms has very high condition number, so numerical instability can be a problem, depending on the values $X_{ij}$.
– Sycorax
Commented Jun 8, 2018 at 14:33

Thanks for your responses, after investigating the shape of the cost function and the behaviour of the gradient descent algorithm here are my findings (which won't surprise any one but some self-learners might find this useful)

### 1) The cost function exhibits a very "flat" bottom

Plotting the convergence of the cost function for various polynomial orders and step sizes shows that the gradient descent algorithm converges very rapidly at first, and then slows down significantly. Here is a plot for $X = [x, x^2]$ but the behaviour is the same for higher orders

My intuition is that a gradient descent algorithm which automatically increases the step size when the cost function is flat would perform much better.

### 2) Since the cost function looks like a flat valley, the starting point matters a lot

In fact, initializing at $[0,0]$ was not a particularly good idea because the value is very close to the bottom of the valley already. Hence the gradient descent struggles to reach to global minimum.

Initializing at random values and comparing results would improve this

### 3) Scipy.optimize algorithms are doing just fine

The 'BFGS' algorithm is in fact very good at finding the global minimum. The issue was that the default tolerance value was too large and the algorithm terminated before reaching the global minimum. Setting the option: 'gtol': 1e-10 leads to convergence in a few hundred iterations (based on 1st order derivative only)

### 4) Numerical instability

Indeed as @Jonny Lomond hinted, very small $x$ values led to extreme numbers for high order polynomials, so I have truncated values close to zero. This improved the behaviour of the algorithms for polynomials order 15 and more

Code here