# coding: utf-8 # # Solving Least-Squares Problems # In[24]: import numpy as np import numpy.linalg as la import scipy.linalg as spla # In[25]: m = 6 n = 4 A = np.random.randn(m, n) b = np.random.randn(m) # Let's try solving that as a linear system using `la.solve`: # In[26]: la.solve(A, b) # OK, let's do QR-based least-squares then. # In[27]: Q, R = la.qr(A) # What did we get? Full QR or reduced QR? # In[28]: Q.shape # In[29]: R.shape # Is that a problem? # * Do we really need the bottom part of $R$? (A bunch of zeros) # * Do we really need the far right part of $Q$? (=the bottom part of $Q^T$) # # ----------------- # OK, so find the minimizing $x$: # In[39]: x = spla.solve_triangular(R, Q.T.dot(b), lower=False) # We predicted that $\|Ax-b\|_2$ would be the same as $\|Rx-Q^Tb\|_2$: # In[45]: la.norm(A.dot(x)-b, 2) # In[47]: la.norm(R.dot(x) - Q.T.dot(b)) # -------------- # Heh--*reduced* QR left out the right half of Q. Let's try again with complete QR: # In[59]: Q2, R2 = la.qr(A, mode="complete") # In[60]: x2 = spla.solve_triangular(R[:n], Q.T[:n].dot(b), lower=False) # In[63]: la.norm(A.dot(x)-b, 2) # In[64]: la.norm(R2.dot(x2) - Q2.T.dot(b)) # Did we get the same `x` both times? # In[69]: x - x2 # Finally, let's compare against the normal equations: # In[70]: x3 = la.solve(A.T.dot(A), A.T.dot(b)) # In[71]: x3 - x # In[ ]: