Linear Regression - Solutions¶

The solutions are written by previous instructors of this course.

In [77]:
import numpy as np
import scipy as sc
import scipy.linalg as la
from matplotlib import pyplot as plt 
%matplotlib inline

Task 1¶

In [78]:
D = np.array([(5,110), (7,114), (9,118), (4,108)])
A = np.vander(D[:,0], 3) 
A
Out[78]:
array([[25,  5,  1],
       [49,  7,  1],
       [81,  9,  1],
       [16,  4,  1]])
In [79]:
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
Out[79]:
array([[ 1.,  1.,  1.,  1.],
       [ 1.,  6.,  5., 30.],
       [ 1.,  7.,  6., 42.],
       [ 1.,  7.,  7., 49.]])

Task 2¶

In [80]:
A = np.array([[0,1],[2,1],[2,1]])
b = np.array([1,1,2]).reshape(-1,1)
In [81]:
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)
In [82]:
leastsq = la.lstsq(A,b)
np.allclose(coeff,leastsq[0])
Out[82]:
True
In [83]:
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
In [84]:
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))
Out[84]:
[<matplotlib.lines.Line2D at 0x7fb2781ed030>]

Task 3¶

In [85]:
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
In [86]:
gauss_seidel(A.T@A,A.T@b)
Out[86]:
array([[0.25],
       [1.  ]])

Task 4¶

  • Gram-Schmidt (QR decomp): we do m inner products of size n, O(mn^2)
  • Gauss: O(n^3)
  • Inversion: O(n^3)

Task 5¶

In [87]:
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)
In [88]:
Q,R = la.qr(A, mode="economic")
coeff = la.solve_triangular(R, Q.T @ b)
In [89]:
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))
Out[89]:
[<matplotlib.lines.Line2D at 0x7fb2731eebc0>]

Task 6¶

In [90]:
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))

Task 7¶

In [91]:
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
In [92]:
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
Out[92]:
(array([0.84938352, 0.84938352, 0.16235147]), 7788.329428639043)
In [93]:
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)
Out[93]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7fb273095450>
In [94]:
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
Out[94]:
(array([ 2.71604846, -1.86666494,  2.71604846, -1.86666494,  1.04922444]),
 7788.347892034986)
In [95]:
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)
Out[95]:
<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.

Task 8¶

In [96]:
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")
In [97]:
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)