# 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$?