# numpy.linalg.lstsq underdetermined case

I would like to understand what I am doing wrong here. I am trying to perform polynomial regression by minimizing the least squares, ||Au-y||^2, where y is the given data and A is the matrix where the i-th line holds [1, x_i, (x_i)^2, ... (x_i)^n-1].

If n is larger than the number of data points, the problem is underdetermined, and I expect the numpy.linalg.lstsq() routine to give any of the infinitely possible solutions. But, as you can see, I don't get a solution at all.


import matplotlib.pyplot as plt
import numpy as np

x = np.array([-6 ,1, 2, 3, 4])     # x data
y = np.array([2, -3, 4, 20, -10])  # y data

A = []
n = 50      # polynomial degree
for i in range(0,n):   # create A matrix of proper form: i-th line is [1, x_i, x_i**2, ...]
A.append(x**i)
A=np.array(A)
A=A.T

u=np.linalg.lstsq(A,y, rcond=None)  # solve underdetermined problem

x_test=np.linspace(-6, 5, 100)   # create more x values for plotting
B=[]
for i in range(0,n):   # same as before, their power matrix
B.append(x_test**i)
B=np.array(B)
B=B.T

plt.plot(x_test, B@u)   # plot polynom
plt.scatter(x,y)   # plot data
plt.show()

print(A@u-y)    # this should be zero vector??
$$$$

• What is A@u and B@u? Have you tried reducing the polynomial degree? – jeremy Nov 14 at 23:59
• It is the the matrix-vector product, nothing special. Reducing, the degree leads to better results. but it should not matter, the thing should be solvable for any polynomial of degree higher than the amount of data points. – SchroedingersLion Nov 15 at 12:07
• What is A@u and B@u? – jeremy Nov 15 at 23:44
• I already answered: It is the matrix vector product: A multiplied by u, and B multiplied by u. – SchroedingersLion Nov 20 at 21:20

I ran your code and got the below plot as well as an almost zero vector of:

[-1.19209290e-07, -2.42083020e-09, -8.94069672e-08,  2.23517418e-07, -6.33299351e-08]


This means numpy successfully curve-fit a line to your data. An order of 6 works a lot better, though:

[-2.79396772e-09 -2.04281037e-14 -1.27897692e-13  4.33431069e-12 -2.27373675e-11]
` numpy won't give you an infinite number of solutions. It will give you an optimally curve-fit solution to your data for the given polynomial order.

• But the first point of the first plot does not lie on the curve. That's what I don't get. Also numbers of 10^-7 in your vector is not what I would call "almost" zero. It should be exactly zero (in floating point arithmetics at least of the order of 10^-15), shouldn't it? – SchroedingersLion Nov 20 at 21:19
• The magnitude of the errors will be different for every order of polynomial. The underlying algorithms in numpy are just doing a best fit. All numpy is telling you is "this is the best 50th power polynomial curve for your data". With this kind of problem, it's important not to overfit the data with a high power polynomial. If the data is random, you are just trying to say "This is my trend". I think you are dealing with the floating point accuracy limit of the estimation algorithms in numpy. You might get lower values if you went and changed the types in Python to be longs. – jeremy Nov 21 at 20:37
• You can try even and odd powers for the polynomial and you might get different results for the number of points the polynomial goes through. – jeremy Nov 21 at 20:50
• My point was for the case where the polynomial degree is larger than the number of data points. In this case, there are analytic solutions so that the polynomial hits all data points exactly. Therefore, aside from floating point errors, I would not have an idea why there are still numbers of order 10^-9 in the error vector. – SchroedingersLion Nov 21 at 21:47