[hw1] Adding code for problem 6

1 files changed, 74 insertions, 25 deletions

diff --git a/hw1/main.py b/hw1/main.py
index 50f1d54..3813fa7 100644
--- a/hw1/main.py
+++ b/hw1/main.py
@@ -1,8 +1,12 @@
 import numpy as np
 import scipy.linalg
 from scipy.optimize import minimize
-from math import sqrt
+from math import sqrt, pi
 import time
+import matplotlib.pyplot as plt
+import seaborn
+
+seaborn.set_style("white")
 
 data = np.loadtxt('CASP.csv', delimiter=',', skiprows=1)
 
@@ -43,22 +47,17 @@ def map_estimate(X, y, sigma=1.0, tau=10.):
                                          check_finite=False)
 
 
-def problem_4():
-    w = map_estimate(X_train, y_train)
-    predictions = np.dot(X_test, w)
-    rmse = np.linalg.norm(predictions - y_test) / sqrt(len(y_test))
-    print w
-    print rmse
+def train_qr(X, y):
+    return map_estimate(X, y)
 
 
-def log_posterior(w, sigma=1.0, tau=10):
-    return (1 / (2 * sigma) * np.linalg.norm(y_train - np.dot(X_train, w)) ** 2
+def log_posterior(w, X, y, sigma=1.0, tau=10):
+    return (1 / (2 * sigma) * np.linalg.norm(y - np.dot(X, w)) ** 2
             + 1 / (2 * tau) * np.linalg.norm(w) ** 2)
 
 
-def gradient(w, sigma=1.0, tau=10):
-    return (- 1 / sigma * np.dot(X_train.T, y_train - np.dot(X_train, w))
-            + 1 / tau * w)
+def gradient(w, X, y, sigma=1.0, tau=10):
+    return (- 1 / sigma * np.dot(X.T, y - np.dot(X, w)) + 1 / tau * w)
 
 
 def verify_gradient(eps=1e-10):
@@ -66,21 +65,71 @@ def verify_gradient(eps=1e-10):
     for i in range(X_train.shape[1]):
         epsilon = np.zeros(X_train.shape[1])
         epsilon[i] = eps
-        print (log_posterior(w + epsilon)
-               - log_posterior(w - epsilon)) / (2 * eps)
-        print gradient(w)[i]
+        print (log_posterior(w + epsilon, X_train, y_train)
+               - log_posterior(w - epsilon, X_train, y_train)) / (2 * eps)
+        print gradient(w, X_train, y_train)[i]
+
+
+def rmse(X, y, w):
+    predictions = np.dot(X, w)
+    return np.linalg.norm(predictions - y) / sqrt(len(y))
+
+
+def train_lbfgs(X, y):
+    x0 = np.zeros(X.shape[1])
+
+    def f(w):
+        return log_posterior(w, X, y)
+
+    def g(w):
+        return gradient(w, X, y)
+
+    return minimize(f, x0, jac=g, method='L-BFGS-B',
+                    options={'maxiter': 100}).x
 
 
-def problem_5():
-    x0 = np.zeros(X_train.shape[1])
-    w = minimize(log_posterior, x0, jac=gradient, method='L-BFGS-B',
-                 options={'maxiter': 100}).x
-    predictions = np.dot(X_test, w)
-    rmse = np.linalg.norm(predictions - y_test) / sqrt(len(y_test))
-    print w
-    print rmse
+def map_features(X, d):
+    A = np.random.normal(0, 1, size=(d, X.shape[1]))
+    b = np.random.uniform(0, 2 * pi, d)
+    return np.cos(np.dot(A, X.T).T + b)
 
 
 if __name__ == "__main__":
-    problem_4()
-    problem_5()
+
+    # problem 4
+    wqr = train_qr(X_train, y_train)
+    print wqr
+    print rmse(X_test, y_test, wqr)
+
+    # problem 5
+    wlbfgs = train_lbfgs(X_train, y_train)
+    print wlbfgs
+    print rmse(X_test, y_test, wlbfgs)
+
+    # problem 6
+    lqr = []
+    llbfgs = []
+    for d in [100, 200, 400, 600]:
+        X = map_features(X_train, d)
+        X_t = map_features(X_test, d)
+
+        t_start = time.time()
+        wqr = train_qr(X, y_train)
+        rqr = time.time() - t_start
+
+        t_start = time.time()
+        wlbfgs = train_lbfgs(X, y_train)
+        rlbfgs = time.time() - t_start
+
+        lqr.append((rqr, rmse(X_t, y_test, wqr)))
+        llbfgs.append((rlbfgs, rmse(X_t, y_test, wlbfgs)))
+
+    plt.figure(figsize=(8, 8))
+    x, y = zip(*lqr)
+    plt.plot(x, y, "-ro", label="QR")
+    x, y = zip(*llbfgs)
+    plt.plot(x, y, "-bo", label="L-BFGS")
+    plt.xlabel("Time (s)")
+    plt.ylabel("RMSE")
+    plt.legend()
+    plt.savefig("plot.pdf", bbox_inches="tight")