The solutions are written by previous instructors of this course.
import numpy as np
import scipy as sc
import scipy.linalg as la
from matplotlib import pyplot as plt
%matplotlib inline
D = np.array([(5,110), (7,114), (9,118), (4,108)])
A = np.vander(D[:,0], 3)
A
array([[25, 5, 1], [49, 7, 1], [81, 9, 1], [16, 4, 1]])
D = np.array([(1,1,6), (6,5,8), (7,6,9), (7,7,12)])
b0 = np.ones((len(D), 1))
b12 = D[:, 0:2]
b3 = D[:, 0]*D[:, 1]
A = np.hstack((b0, b12, b3[...,None]))
A
array([[ 1., 1., 1., 1.], [ 1., 6., 5., 30.], [ 1., 7., 6., 42.], [ 1., 7., 7., 49.]])
A = np.array([[0,1],[2,1],[2,1]])
b = np.array([1,1,2]).reshape(-1,1)
def qr_solve(A, b):
# from slides we have: Rz = Q^Ty
Q,R = la.qr(A, mode="economic")
coeff = la.solve_triangular(R, Q.T @ b)
return coeff
coeff = qr_solve(A,b)
leastsq = la.lstsq(A,b)
np.allclose(coeff,leastsq[0])
True
pl = sc.poly1d(coeff[:,0])
E = sum((b[:,0]-pl(A[:,0]))**2) # sum of squared errors, ie, error on the training set
print(E, leastsq[1], la.norm(b-A@coeff,ord=2)**2)
0.5 [0.5] 0.5000000000000001
plt.plot(A[:,0],b,'o')
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
x = np.linspace(0,2,50)
plt.plot(x,pl(x))
[<matplotlib.lines.Line2D at 0x7fb2781ed030>]
def gauss_seidel(A, b, tolerance=1e-10, max_iterations=10000):
x = np.zeros_like(b, dtype=np.double)
for k in range(max_iterations):
x_old = x.copy()
for i in range(A.shape[0]):
x[i] = (b[i] - np.dot(A[i,:i], x[:i]) - np.dot(A[i,(i+1):], x_old[(i+1):])) / A[i ,i]
if np.linalg.norm(x - x_old, ord=np.inf) / np.linalg.norm(x, ord=np.inf) < tolerance:
break
return x
gauss_seidel(A.T@A,A.T@b)
array([[0.25], [1. ]])
D=np.load("../assets/housing.npy")
A=np.column_stack([D[:,0],np.ones(D.shape[0])]).reshape(-1,2)
b=np.array(D[:,1],).reshape(-1,1)
Q,R = la.qr(A, mode="economic")
coeff = la.solve_triangular(R, Q.T @ b)
pl = np.poly1d(coeff[:,0])
plt.plot(A[:,0],b,'o')
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
x = np.linspace(0,16,50)
plt.plot(x,pl(x))
[<matplotlib.lines.Line2D at 0x7fb2731eebc0>]
plt.figure(figsize=(13,3))
b = D[:,1].reshape(-1,1)
k = [3,6,9,12]
for i in range(4):
plt.subplot(1,4,i+1)
A=np.vander(D[:,0],k[i]+1)
coeff = la.lstsq(A, b)[0]
pl = np.poly1d(coeff[:,0])
x = np.linspace(0,16,50)
plt.plot(D[:,0],b,'o')
plt.plot(x,pl(x))
np.random.seed(5)
m=50
x1 = np.linspace(0, 1, m)
x2 = np.linspace(0, 1, m)
y = x1**2-7*x1+2*x2**2+5*x2+np.random.exponential(1,m) # the unknown target function
A=np.column_stack([x1,x2,np.ones(m)]).reshape(-1,3)
coeff = la.lstsq(A,y)[0]
coeff, la.norm(b-A@coeff,ord=2) # root squared errors
(array([0.84938352, 0.84938352, 0.16235147]), 7788.329428639043)
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x1, x2, y, marker='o')
X1, X2 = np.meshgrid(x1, x2)
Y = coeff[0]*X1 + coeff[1]*X2 + coeff[2]
ax.plot_surface(X1, X2, Y)
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7fb273095450>
A=np.column_stack([x1**2, x1, x2**2, x2, np.ones(m)]).reshape(-1,5)
coeff = la.lstsq(A,y)[0]
coeff, la.norm(b-A@coeff,ord=2) # root squared errors
(array([ 2.71604846, -1.86666494, 2.71604846, -1.86666494, 1.04922444]), 7788.347892034986)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x1, x2, y, marker='o')
Y = coeff[0]*X1**2 + coeff[1]*X1 +coeff[2]*X2**2 + coeff[3]*X2 + coeff[4]
ax.plot_surface(X1, X2, Y)
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7fb272ddd5a0>
Looking at the root squared erros on the training set, the first model seems better than the second. However the final decision should be taken on the basis of a test set of points different from the training set.
def plot_ellipse(a, b, c, d, e):
"""Plot an ellipse of the form ax^2 + bx + cxy + dy + ey^2 = 1."""
theta = np.linspace(0, 2*np.pi, 200)
cos_t, sin_t = np.cos(theta), np.sin(theta)
A = a*(cos_t**2) + c*cos_t*sin_t + e*(sin_t**2)
B = b*cos_t + d*sin_t
r = (-B + np.sqrt(B**2 + 4*A)) / (2*A)
plt.plot(r*cos_t, r*sin_t, lw=2)
plt.gca().set_aspect("equal", "datalim")
D = np.load("../assets/ellipse.npy")
x, y = D[:, 0], D[:, 1]
A = np.column_stack((x**2, x, x*y, y, y**2))
b = np.ones(len(D))
a, b, c, d, e = qr_solve(A, b)
plt.plot(x,y,'o')
plot_ellipse(a, b, c, d, e)