# coding: utf-8 # # Least Squares using the SVD # In[2]: import numpy as np import numpy.linalg as la import scipy.linalg as spla # In[19]: # tall and skinny w/nullspace np.random.seed(12) A = np.random.randn(6, 4) b = np.random.randn(6) A[3] = A[4] + A[5] A[1] = A[5] + A[1] A[2] = A[3] + A[1] A[0] = A[3] + A[1] # ### Part I: Singular least squares using QR # # Let's see how successfully we can solve the least squares problem **when the matrix has a nullspace** using QR: # In[4]: Q, R = la.qr(A) # In[5]: R.round(3) # We can choose x_qr[3] as we please: # In[6]: x_qr = np.zeros(A.shape[1]) # In[20]: x_qr[3] = 0 # In[21]: QTbnew = Q.T.dot(b)[:3,] - R[:3, 3] * x_qr[3] x_qr[:3] = spla.solve_triangular(R[:3,:3], QTbnew, lower=False) # Let's take a look at the residual norm and the norm of x_qr: # In[22]: R.dot(x_qr)-Q.T.dot(b)[:4] # In[23]: la.norm(A.dot(x_qr)-b, 2) # In[24]: la.norm(x_qr, 2) # Choose a different x_qr[3] and compare residual and norm of x_qr. # -------------- # ### Part II: Solving least squares using the SVD # Now compute the SVD of $A$: # In[25]: U, sigma, VT = la.svd(A) # Make a matrix Sigma of the correct size: # In[26]: Sigma = np.zeros(A.shape) Sigma[:4,:4] = np.diag(sigma) # And check that we've actually factorized A: # In[27]: (U.dot(Sigma).dot(VT) - A).round(4) # Now define Sigma_pinv as the "pseudo-"inverse of Sigma, where "pseudo" means "don't divide by zero": # In[28]: Sigma_pinv = np.zeros(A.shape).T Sigma_pinv[:3,:3] = np.diag(1/sigma[:3]) Sigma_pinv.round(3) # Now compute the SVD-based solution for the least-squares problem: # In[29]: x_svd = VT.T.dot(Sigma_pinv).dot(U.T).dot(b) # In[30]: la.norm(A.dot(x_svd)-b, 2) # In[31]: la.norm(x_svd) # * What do you observe about $\|\text{x_svd}\|_2$ compared to $\|\text{x_qr}\|_2$? # * Is $\|\text{x_svd}\|_2$ compared to $\|\text{x_qr}\|_2$?