CODEWITHSUNDEEP
pythonpython-3.xnumpyleast-squares
user7086216
user7086216

Reputation:

Least Squares: Python

I am trying to implement least squares:

I have: $y=\theta\omega$

The least square solution is \omega=(\theta^{T}\theta)^{-1}\theta^{T}y

I tryied: 

import numpy as np    
def least_squares1(y, tx):
        """calculate the least squares solution."""
        w = np.dot(np.linalg.inv(np.dot(tx.T,tx)), np.dot(tx.T,y))

        return w

The problem is that this method becomes quickly unstable (for small problems its okay)

I realized that, when I compared the result to this least square calculation:

 import numpy as np 
 def least_squares2(y, tx):
        """calculate the least squares solution."""
        a = tx.T.dot(tx)
        b = tx.T.dot(y)
        return np.linalg.solve(a, b)

Compare both methods: I tried to fit data with a polynomial of degree 12 [1, x,x^2,x^3,x^4...,x^12]

First method:

enter image description here

Second method:

enter image description here

Do you know why the first method diverges for large polynomials ?

P.S. I only added "import numpy as np" for your convinience, if you want to test the functions.

Upvotes: 0

Views: 3236

Answers (2)

user1911226
user1911226

Reputation:

I am going to improve on what was said before. I answered this yesterday. The problem with higher order polynomials is something called Runge's phenomena. The reason why the person resorted orthogonal polynomials which are known as Hermite polynomials is that they attempt to get rid of the Gibbs phenomenon which is an adverse oscillatory effect when Fourier series methods are applied to non-periodic signals.

You can sometimes improve under the conditioning be resorting to regularizing methods if the matrix is low rank as I did in the other post. Other parts may be due to smoothness properties of the vector.

Upvotes: 0

dmuir
dmuir

Reputation: 4431

There are three points here:

One is that it is generally better (faster, more accurate) to solve linear equations rather than to compute inverses.

The second is that it's always a good idea to use what you know about a system of equations (e.g. that the coefficient matrix is positive definite) when computing a solution, in this case you should use numpy.linalg.lstsq

The third is more specifically about polynomials. When using monomials as a basis, you can end up with a very poorly conditioned coefficient matrix, and this will mean that numerical errors tend to be large. This is because, for example, the vectors x->pow(x,11) and x->pow(x,12) are very nearly parallel. You would get a more accurate fit, and be able to use higher degrees, if you were to use a basis of orthogonal polynomials, for example https://en.wikipedia.org/wiki/Chebyshev_polynomials or https://en.wikipedia.org/wiki/Legendre_polynomials

Upvotes: 3

Related Questions