title: Linear least squares in scipy - accuracy of QR factorization vs other methods
body: <p>I have tried solving a linear least squares problem Ax = b in scipy using the following methods:</p>
  
  <pre><code>x = numpy.linalg.inv(A.T.dot(A)).dot(A.T).dot(b) #Usually not recommended
  </code></pre>
  
  <p>and</p>
  
  <pre><code>x = numpy.linalg.lstsq(A, b)
  </code></pre>
  
  <p>both give almost identical results. I also tried manually using the QR algorithm to do so ie:</p>
  
  <pre><code>Qmat, Rmat = la.qr(A)
  
  bpr = dot(Qmat.T,b)
  n=len(bpr)
  x = np.zeros(n)
  for i in xrange(n-1, -1,-1):
      x[i] = bpr[i]
      for j in xrange(i+1, n):
          x[i] -= Rmat[i, j]*x[j]
      x[i] /= Rmat[i,i]
  </code></pre>
  
  <p>This method, however, gives very inaccurate results (errors on the order of 1e-2). Have I made a n00b mistake with the code or maths? Or, is the an issue with the method, or scipy itself?</p>
  
  <p>My numpy version is 1.6.1 (the mkl compiled version from <a href="http://www.lfd.uci.edu/~gohlke/pythonlibs/" rel="nofollow">http://www.lfd.uci.edu/~gohlke/pythonlibs/</a>), with Python 2.7.3 on x86_64.</p>