University:
Stony Brook University
Course:
CSE 114 | Introduction to Object-Oriented Programming
Academic year:

2019
Views:

124

Pages:

9
Author:

Maverick Martin

Part 3: The cost of solving the normal equations This notebook helps you explore the execution time cost of solving the normal equations, This notebook only has one exercise, but it is not graded. So, you should complete it for your own edification. In [1]: import numpy as np import matplotlib.pyplot as plt %matplotlib inline Scalability with the problem size To start, here is some code to help generate synthetic problems of a certain size, namely, number of predictors. The , where is the number of observations and comes from our usual dummy coefficient for a non-zero intercept. We will also implement a linear least squares solver, estimate_coeffs(), that simply calls Numpy's lstsq() routine. In [2]: def generate_model (n): """Returns a set of (random) n+1 linear model coefficients.""" return np.random.rand (n+1, 1) def generate_data (m, theta, sigma=1.0/(2**0.5)): """ Generates 'm' noisy observations for a linear model whose predictor (non-intercept) coefficients are given in 'theta'. Decrease 'sigma' to decrease the amount of noise. """ assert (type (theta) is np.ndarray) and (theta.ndim == 2) and (theta.shape[1] == 1) n = len (theta) X = np.random.rand (m, n) X[:, 0] = 1.0 y = X.dot (theta) + sigma*np.random.randn (m, 1) return (X, y) def estimate_coeffs(X, y): """ Solves X*theta = y by a linear least squares method. """ result = np.linalg.lstsq (X, y, rcond = None) theta = result[0] return theta In [3]: # Demo the above routines for a 2-D dataset. m = 50 theta_true = generate_model (1) (X, y) = generate_data (m, theta_true, sigma=0.1) print ("Dimensions of X:", X.shape) print ("Dimensions of theta_true:", theta_true.shape) print ("Dimensions of y:", y.shape) print ("Condition number of X: ", np.linalg.cond (X)) print ("True model coefficients:", theta_true.T) theta = estimate_coeffs (X, y) the print ("Estimated model coefficients:", theta.T) fig = plt.figure() ax1 = fig.add_subplot(111) ax1.plot (X[:, 1], y, 'b+') # Noisy observations ax1.plot (X[:, 1], X.dot (theta), 'r*') # Fit ax1.plot (X[:, 1], X.dot (theta_true), 'go') # True solution Dimensions of X: (50, 2) Dimensions of theta_true: (2, 1) Dimensions of y: (50, 1) Condition number of X: 4.1056142618693565 True model coefficients: [[0.74734942 0.76050905]] Estimated model coefficients: [[0.73296025 0.74693774]] Out[3]: [] Benchmark varying . Let's benchmark the time to compute running time scale with ? when the dimension of each point is fixed but the number In [4]: # Benchmark, as 'm' varies: n = 32 # dimension M = [100, 1000, 10000, 100000, 1000000] times = [0.] * len (M) for (i, m) in enumerate (M): theta_true = generate_model (n) (X, y) = generate_data (m, theta_true, sigma=0.1) t = %timeit -o estimate_coeffs (X, y) times[i] = t.best 804 µs ± 7.11 3.13 ms ± 931 10.8 ms ± 899 107 ms ± 13.7 1.39 s ± 47.6 µs µs µs ms ms per per per per per loop loop loop loop loop (mean (mean (mean (mean (mean ± ± ± ± ± std. std. std. std. std. dev. dev. dev. dev. dev. of of of of of 7 7 7 7 7 runs, runs, runs, runs, runs, 1000 loops each) 100 loops each) 100 loops each) 10 loops each) 1 loop each) In [5]: t_linear = [times[0]/M[0]*m for m in M] fig = plt.figure() ax1 = fig.add_subplot(111) ax1.loglog (M, times, 'bo') ax1.loglog (M, t_linear, 'r--') ax1.set_xlabel ('m (number of observations)') fig.suptitle ('Running time (fixed number of predictors)') Out[5]: Text(0.5,0.98,'Running time (fixed number of predictors)') of points varies. How does the Exercise 0 (ungraded). Now fix the number of observations but vary the dimension . How does time scale with ? Complete the benchmark code below to find out. In particular, given the array N[:], compute an array, times[:], such that times[i] is the running time for a problem of size m (N[i]+1). Hint: You can adapt the preceding benchmark. Also, note that the code cell following the one immediately below will plot your results against and . In [6]: N = [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024] m = 5000 times = [0.] * len (N) # Implement a benchmark to compute the time, # `times[i]`, to execute a problem of size `N[i]`. for (i, n) in enumerate (N): ### BEGIN SOLUTION theta_true = generate_model (n) (X, y) = generate_data (m, theta_true, sigma=0.1) t = %timeit -o estimate_coeffs (X, y) times[i] = t.best ### END SOLUTION 347 µs ± 238 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each) 520 µs ± 372 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each) 998 µs ± 75.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) 1.96 ms ± 171 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) 5.82 ms ± 454 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) 16.7 ms ± 1.78 ms per loop (mean ± std. dev. of 7 runs, 100 loops each) 40.1 ms ± 7.43 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) 97.9 ms ± 12.7 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) 249 ms ± 366 µs per loop (mean ± std. dev. of 7 runs, 1 loop each) 876 ms ± 178 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) In [7]: t_linear = [times[0]/N[0]*n for n in N] t_quadratic = [times[0]/N[0]/N[0]*n*n for n in N] fig = plt.figure() ax1 = fig.add_subplot(111) ax1.loglog (N, times, 'bo') ax1.loglog (N, t_linear, 'r--') ax1.loglog (N, t_quadratic, 'g--') ax1.set_xlabel ('n (number of predictors)') fig.suptitle ('Running time (fixed number of observations)') Out[7]: Text(0.5,0.98,'Running time (fixed number of observations)') Thus, the empirical scaling appears to be pretty reasonable, being roughly linear in that in practical regression problems. Fin! If you've gotten this far without errors, your notebook is ready to submit. . And while being quadratic in sounds bad, one expects (or hopes!) Part 4: "Online" linear regression When you are trying to fit a model to data and you get to see all of the data at once, we refer to the problem as an offline or batch problem, and you would try to use certain algorithms to compute the fit that can take advantage of the fact that you have a lot of available data. But what if you only get to see one or a few data points at a time? In that case, you might want to get an initial model from whatever data you've got, and gradually improve the model as you see new data points. In this case, we refer to the problem as being an online problem. The goal of this notebook is to introduce you to online algorithms. You'll start by reviewing the offline linear regression problem, and then look at its online variant. The neat thing about the online method is that you can derive it using all the tools you already have at your disposal, namely, multivariate calculus. In [1]: import numpy as np import matplotlib.pyplot as plt %matplotlib inline Review: Offline or batch linear regression Let's start with a quick review of the linear regression problem: given a response vector, , and a data matrix ---whose rows are observations and columns are variables---the problem is to find the best linear model, , where is the vector of best-fit model parameters that we wish to compute. Computing it using a conventional batch linear least squares method has an asymptotic running time of . To start, here is some code to help generate synthetic problems of a certain size, namely, number of predictors. The comes from our usual dummy coefficient for a non-zero intercept. , where is the number of observations and the In [2]: def generate_model (n): """Returns a set of (random) n+1 linear model coefficients.""" return np.random.rand (n+1, 1) def generate_data (m, theta, sigma=1.0/(2**0.5)): """ Generates 'm' noisy observations for a linear model whose predictor (non-intercept) coefficients are given in 'theta'. Decrease 'sigma' to decrease the amount of noise. """ assert (type (theta) is np.ndarray) and (theta.ndim == 2) and (theta.shape[1] == 1) n = len (theta) X = np.random.rand (m, n) X[:, 0] = 1.0 y = X.dot (theta) + sigma*np.random.randn (m, 1) return (X, y) def estimate_coeffs (X, y): """ Solves X*theta = y by a linear least squares method. """ result = np.linalg.lstsq (X, y, rcond=None) theta = result[0] return theta In [3]: def rel_diff(x, y, ord=2): """ Computes ||x-y|| / ||y||. Uses 2-norm by default; override by setting 'ord'. """ return np.linalg.norm (x - y, ord=ord) / np.linalg.norm (y, ord=ord) An online algorithm The empirical scaling of linear least squares appears to be pretty good, being roughly linear in or at worst quadratic in . But there is still a downside in time and storage: each time there is a change in the data, you appear to need to form the data matrix all over again and recompute the solution from scratch, possibly touching the entire data set again! This begs the question, is there a way to incrementally update the model coefficients whenever a new data point, or perhaps a small batch of new data points, arrives? Such a procedure would be considered incremental or online, rather than batched or offline. Setup: Key assumptions and main goal. In the discussion that follows, assume that you only get to see the observations one-at-a-time. Let denote the current observation. (Relative to our previous notation, this tuple is just element We will use to denote a row of since we previously used of and row to denote column of of . . That is, where the first form is our previous "columns-view" representation and the second form is our "rows-view." Additionally, assume that, at the time the -th observation arrives, you start with a current estimate of the parameters, reason you need to refer to element of that vector, use below, further assume that you throw out . You will then compute a new estimate, once you have using , which is a vector. If for whatever and . For the discussion . As for your goal, recall that in the batch setting you start with all the observations, . From this starting point, you may estimate the linear regression model's parameters, , by solving . In the online setting, you compute estimates one at a time. After seeing all observations in , your goal is to compute an . An intuitive (but flawed) idea. Indeed, there is a technique from the signal processing literature that we can apply to the linear regression problem, known as the least mean square (LMS) algorithm. Before describing it, let's start with an initial idea. Suppose that you have a current estimate of the parameters, , when you get a new sample, Ideally, this error would be zero. So, let's ask if there exists a correction, . The error in your prediction will be, , such that Then, you could compute a new estimate of the parameter by . This idea has a major flaw, which we will discuss below. But before we do, please try the following exercise. Mental exercise (no points). Verify that the following choice of would make the preceding equation true. Refining (or rather, "hacking") the basic idea: The least mean square (LMS) procedure. The basic idea sketched above has at least one major flaw: the choice of might allow you to correctly predict from and the new estimate , but there is no guarantee that this new estimate preserves the quality of predictions made at all previous iterations! There are a number of ways to deal with this problem, which includes carrying out an update with respect to some (or all) previous data. However, there is also a simpler "hack" that, though it might require some parameter tuning, can be made to work in practice. That hack is as follows. Rather than using as computed above, let's compute a different update that has a "fudge" factor, : A big question is how to choose . There is some analysis out there that can help. We will just state the results of this analysis without proof. Let be the largest eigenvalue of . The result is that as the number of samples condition will eventually converge to the best least-squares estimator of , that is, the estimate of , any choice of that satisfies the following you would have gotten by solving the linear least squares problem with all of the data. This condition is not very satisfying, because you cannot really know until you've seen all the data, whereas we would like to apply this procedure online as the data arrive. Nevertheless, in practice you can imagine hybrid schemes that, given a batch of data points, use the QR fitting procedure to get a starting estimate for as well as to estimate a value of to use for all future updates. Summary of the LMS algorithm. To summarize, the algorithm is as follows: Choose any initial guess such as Choose any initial guess, , such as For each observation , do the update: . , where . Trying out the LMS idea Now you should implement the LMS algorithm and see how it behaves. To start, let's generate an initial 1-D problem (2 regression coefficients, a slope, and an intercept), and solve it using the batch procedure. Recall that we need a value for , for which we have an upper-bound of need to do something different. . Let's cheat by computing it explicitly, even though in practice we would In [4]: m = 100000 n = 1 theta_true = generate_model(n) (X, y) = generate_data(m, theta_true, sigma=0.1) print("Condition number of the data matrix:", np.linalg.cond(X)) theta = estimate_coeffs(X, y) e_rel = rel_diff(theta, theta_true) print("Relative error:", e_rel) Condition number of the data matrix: 4.400944337141808 Relative error: 0.0016564742343831158 In [5]: LAMBDA_MAX = max(np.linalg.eigvals(X.T.dot(X))) print(LAMBDA_MAX) 126970.43769457837 Exercise 1 (5 points). Implement the online LMS algorithm in the code cell below where indicated. It should produce a final parameter estimate, theta_lms, as a column vector. In addition, the skeleton code below uses rel_diff() to record the relative difference between the estimate and the true vector, storing the -th relative difference in rel_diffs[k]. Doing so will allow you to see the convergence behavior of the method. Lastly, to help you out, we've defined a constant in terms of that you can use for . In practice, you would only maintain the current estimate, or maybe just a few recent estimates, rather than all of them. Since we want to inspect these vectors later, go ahead and store them all. In [6]: PHI = 1.99 / LAMBDA_MAX # Fudge factor rel_diffs = np.zeros((m+1, 1)) theta_k = np.zeros((n+1)) for k in range(m): rel_diffs[k] = rel_diff(theta_k, theta_true) # Implement the online LMS algorithm. # Take (y[k], X[k, :]) to be the k-th observation. ### BEGIN SOLUTION x_k = X[k, :] r_k = y[k] - x_k.T.dot(theta_k) delta_k = PHI * r_k * x_k theta_k = theta_k + delta_k ### END SOLUTION theta_lms = theta_k rel_diffs[m] = rel_diff(theta_lms, theta_true) Let's compare the true coefficients against the estimates, both from the batch algorithm and the online algorithm. The values of the variables below might change if the notebooks are re-run from start. In [7]: print (theta true T) In [7]: print (theta_true.T) print (theta.T) print (theta_lms.T) print("\n('Passed' -- this cell appears to run without error, but we aren't checking the solution.)") [[0.87103594 0.46153527]] [[0.87172037 0.46005275]] [0.74984566 0.40192802] ('Passed' -- this cell appears to run without error, but we aren't checking the solution.) Let's also compute the relative differences between each estimate Theta[:, k] and the true coefficients theta_true, measured in the two-norm, to see if the estimate is converging to the truth. In [8]: plt.plot(range(len(rel_diffs)), rel_diffs) Out[8]: [] You should see it converging, but not especially quickly. Finally, if the dimension is n=1, let's go ahead and do a sanity-check regression fit plot. The plot can change if the notebooks are re-run from start. In [9]: STEP = int(X.shape[0] / 500) if n == 1: fig = plt.figure() ax1 = fig.add_subplot(111) ax1.plot(X[::STEP, 1], y[::STEP], 'b+') # blue - data ax1.plot(X[::STEP, 1], X.dot(theta_true)[::STEP], 'r*') # red - true ax1.plot(X[::STEP, 1], X.dot(theta)[::STEP], 'go') # green - batch ax1.plot(X[::STEP, 1], X.dot(theta_lms)[::STEP], 'mo') # magenta - pure LMS else: print("Plot is multidimensional; I live in Flatland, so I don't do that.") Exercise 2 (ungraded, optional). We said previously that, in practice, you would probably do some sort of hybrid scheme that mixes full batch updates (possibly only initially) and incremental updates. Implement such a scheme and describe what you observe. You might observe a different plot each time the cell is re-run. In [10]: # Setup problem and compute the batch solution m = 100000 n = 1 theta_true = generate_model(n) (X, y) = generate_data(m, theta_true, sigma=0.1) theta_batch = estimate_coeffs(X, y) # Your turn, below: Implement a hybrid batch-LMS solution # assuming you observe the first few data points all at # once, and then see the remaining points one at a time. ### BEGIN SOLUTION def lms(X, y, theta_0=None, PHI=1e-6): """Implements the LMS algorithm for the system X theta = y.""" print(X.shape) (m, n) = X.shape if theta_0 is None: theta_k = np.zeros((n)) else: theta_k = theta_0 theta_k = np.reshape(theta_k, (n, 1)) for k in range(m): x_k = X[k:k+1, :].T r_k = y[k] - x_k.T.dot(theta_k) delta_k = PHI * r_k * x_k theta_k = theta_k + delta_k return theta_k # Pure LMS solution START = 100 lambdas = np.linalg.eigvals(X[:START, :].T.dot(X[:START, :])) lambda_max = max(lambdas) lambda_min = min(lambdas) phi = 2.0 / (lambda_max + lambda_min) # See: https://en.wikipedia.org/wiki/Least_mean_squares_filter theta_lms = lms(X, y, PHI=phi) # Batch on first START observations + LMS on the rest: theta_batch_START = estimate_coeffs(X[:START, :], y[:START]) print(theta_batch_START) theta_hybrid = lms(X[START:, :], y[START:], theta_0=theta_batch_START, PHI=phi) print(theta_hybrid) def calc_ssqr(theta, X, y): r = X.dot(theta) - y return r.T.dot(r) print("=== Sum of squared errors ===") print("True:", calc_ssqr(theta_true, X, y)) print("Batch:", calc_ssqr(theta_batch, X, y)) print("Batch-{}:".format(START), calc_ssqr(theta_batch_START, X, y)) print("LMS(phi={}):".format(phi), calc_ssqr(theta_lms, X, y)) print("Batch (k <= %d) + LMS:" % START, calc_ssqr(theta_hybrid, X, y)) print("=== Relative error in the coefficients (theta) ===") print("Batch:", rel_diff(theta_batch, theta_true)) print("Batch-{}:".format(START), rel_diff(theta_batch_START, theta_true)) print("LMS(phi={}):".format(phi), rel_diff(theta_lms, theta_true)) print("Batch (k <= %d) + LMS:" % START, rel_diff(theta_hybrid, theta_true)) STEP = int(X.shape[0] / 100) if n == 1: fig = plt.figure() ax1 = fig.add_subplot(111) ax1.plot(X[::STEP, 1], y[::STEP], 'b+') # blue - data ax1.plot(X[::STEP, 1], X.dot(theta_true)[::STEP], 'r*') # red - true ax1.plot(X[::STEP, 1], X.dot(theta_batch)[::STEP], 'go') # green - batch ax1.plot(X[::STEP, 1], X.dot(theta_lms)[::STEP], 'mo') # magenta - pure LMS ax1.plot(X[::STEP, 1], X.dot(theta_hybrid)[::STEP], 'yx') # yellow - hybrid else: print ("Plot is multidimensional; I live in Flatland, so I don't do that.") ### END SOLUTION (100000, 2) [[0.54400102] [0.33553129]] (99900, 2) [[0.51075904] [0.39327854]] === Sum of squared errors === True: [[994.76733987]] Batch: [[994.76638474]] Batch-100: [[1040.65200804]] LMS(phi=0.014605068287428212): [[999.28578087]] Batch (k <= 100) + LMS: [[999.28578087]] === Relative error in the coefficients (theta) === Batch: 0.00025460403740161595 Batch-100: 0.12047813595445639 LMS(phi=0.014605068287428212): 0.018194600645439905 Batch (k <= 100) + LMS: 0.018194600645439905 Fin! If you've gotten this far without errors, your notebook is ready to submit.