title: Least Squares: Python
body: <p>I am trying to implement least squares:</p>
  
  <p>I have:  $y=\theta\omega$</p>
  
  <p>The least square solution is \omega=(\theta^{T}\theta)^{-1}\theta^{T}y</p>
  
  <p>I tryied: </p>
  
  <pre><code>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
  </code></pre>
  
  <p><strong>The problem is that this method becomes quickly unstable</strong>
  (for small problems its okay)</p>
  
  <p>I realized that, when I compared the result to this least square calculation:</p>
  
  <pre><code> 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)
  </code></pre>
  
  <p>Compare both methods:
  I tried to fit data with a polynomial of degree 12 [1, x,x^2,x^3,x^4...,x^12]</p>
  
  <p>First method:</p>
  
  <p><a href="https://i.stack.imgur.com/NqZ37.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/NqZ37.png" alt="enter image description here"></a></p>
  
  <p>Second method:</p>
  
  <p><a href="https://i.stack.imgur.com/RU9B5.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/RU9B5.png" alt="enter image description here"></a></p>
  
  <p><strong>Do you know why the first method diverges for large polynomials ?</strong></p>
  
  <p>P.S. I only added "import numpy as np" for your convinience, if you want to test the functions.</p>