Sample Dataset (LSD) Part 1

Sample solutions Part 0: Sample dataset (LSD) In 1968, Wagner Agahajanian, and Bing conducted a study to determine whether you could improve a student's math test scores using lysergic acid diethylamide, also known as "LSD." Here is the original data sources. The code cell below downloads the file from an alternative location, for compatibility with the Azure Notebook or Vocareum platforms you are using. Raw data: http://www.stat.ufl.edu/~winner/data/lsd.dat (http://www.stat.ufl.edu/~winner/data/lsd.dat) Data file description: http://www.stat.ufl.edu/~winner/data/lsd.txt (http://www.stat.ufl.edu/~winner/data/lsd.txt) In [ ]: from pandas import read_fwf from IPython.display import display import import import import requests os hashlib io def on_vocareum(): return os.path.exists('.voc') if on_vocareum(): URL_BASE = "https://cse6040.gatech.edu/datasets/" DATA_PATH = "../resource/asnlib/publicdata/" else: URL_BASE = "https://github.com/cse6040/labs-fa17/raw/master/datasets/" DATA_PATH = "" def download(file, local_dir="", url_base=URL_BASE, checksum=None): local_file = "{}{}".format(local_dir, file) if not os.path.exists(local_file): url = "{}{}".format(url_base, file) print("Downloading: {} ...".format(url)) r = requests.get(url) with open(local_file, 'wb') as f: f.write(r.content) if checksum is not None: with io.open(local_file, 'rb') as f: body = f.read() body_checksum = hashlib.md5(body).hexdigest() assert body_checksum == checksum, \ "Downloaded file '{}' has incorrect checksum: '{}' instead of '{}'".format(local_file, body_checksum, checksum) print("'{}' is ready!".format(file)) datasets = {'lsd.dat': '4c119057baf86cff8da03d825d7ce141'} for filename, checksum in datasets.items(): download(filename, local_dir=DATA_PATH, url_base=URL_BASE, checksum=checksum) print("\n(All data appears to be ready.)") Let's take a look at the data, first as a table and then using a scatter plot. In [ ]: df = read_fwf('{}lsd.dat'.format(DATA_PATH), colspecs=[(0, 4), (7, 13)], names=['lsd_concentration', 'exam_score']) display(df) In [ ]: from matplotlib.pyplot import scatter, xlabel, title, plot %matplotlib inline scatter(df['lsd_concentration'], df['exam_score']) xlabel ('LSD Tissue Concentration') title ('Shocking news: Math scores degrade with increasing LSD!'); Fitting a model Exercise 0 (2 points). Complete the function below so that it computes arrays y[:] for the responses and x[:] for the predictor. and for the univariate model, , given observations stored as NumPy Use the linear regression formulas derived in class. In [ ]: def linreg_fit(x, y): """Returns (alpha, beta) s.t. y ~ alpha*x + beta.""" from numpy import ones m = len(x) ; assert len(y) == m ### BEGIN SOLUTION u = ones(m) alpha = x.dot(y) - u.dot(x)*u.dot(y)/m alpha /= x.dot(x) - (u.dot(x)**2)/m beta = u.dot(y - alpha*x)/m ### END SOLUTION return (alpha, beta) # Compute the coefficients for the LSD data: x, y = df['lsd_concentration'], df['exam_score'] alpha, beta = linreg_fit(x, y) print("alpha:", alpha) print("beta:", beta) In [ ]: # Test cell: `linreg_fit_test` x, y = df['lsd_concentration'], df['exam_score'] alpha, beta = linreg_fit(x, y) r = alpha*x + beta - y ssqr = r.dot(r) ssqr_ex = 253.88132881 from numpy import isclose assert isclose(ssqr, ssqr_ex, rtol=.01), "Sum-of-squared residuals is {} instead of {}.".format(ssqr, ssqr_ex) print("\n(Passed!)") In [ ]: from numpy import linspace, floor, ceil # Two points make a line: x_fit = linspace(floor(x.min()), ceil(x.max()), 2) y_fit = alpha*x_fit + beta scatter(x, y, marker='o') plot(x_fit, y_fit, 'r--') xlabel('LSD Tissue Concentration') title('Best-fit linear model'); Fin! If you've gotten this far without errors, your notebook is ready to submit. Part 1: Gradients example This notebook is designed to illustrate the concept of the gradient. Let be a two-dimensional vector, and let Exercise 0 (1 point). Implement the Python function, f(x0, x1), so that it computes In [1]: def f(x0, x1): ### BEGIN SOLUTION return x0*x0 + x1*x1 ### END SOLUTION In [2]: # Test cell: `f_test` from numpy import sin, cos, vectorize, isclose from numpy.random import randn f_vec = vectorize(f) theta = randn(1000) assert all(isclose(f_vec(sin(theta), cos(theta)), 1.0)) print("\n(Passed!)") (Passed!) The gradient Let's create a mesh of coordinate values: In [3]: from numpy import linspace, meshgrid x0 = linspace(-2, 2, 11) x1 = linspace(-2, 2, 11) X0, X1 = meshgrid(x0, x1) In [4]: print("X0:\n", X0) print("X1:\n", X1) X0: [[-2. [-2. [-2. [-2. [-2. [-2. [-2. [-2. [-2. [-2. [-2. X1: [[-2. [-1.6 [-1.2 [-0.8 [-0.4 [ 0. [ 0.4 [ 0.8 [ 1.2 [ 1.6 [ 2. -1.6 -1.2 -0.8 -0.4 -1.6 -1.2 -0.8 -0.4 -1.6 -1.2 -0.8 -0.4 -1.6 -1.2 -0.8 -0.4 -1.6 -1.2 -0.8 -0.4 -1.6 -1.2 -0.8 -0.4 -1.6 -1.2 -0.8 -0.4 -1.6 -1.2 -0.8 -0.4 -1.6 -1.2 -0.8 -0.4 -1.6 -1.2 -0.8 -0.4 -1.6 -1.2 -0.8 -0.4 -2. -1.6 -1.2 -0.8 -0.4 0. 0.4 0.8 1.2 1.6 2. -2. -1.6 -1.2 -0.8 -0.4 0. 0.4 0.8 1.2 1.6 2. -2. -1.6 -1.2 -0.8 -0.4 0. 0.4 0.8 1.2 1.6 2. -2. -1.6 -1.2 -0.8 -0.4 0. 0.4 0.8 1.2 1.6 2. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -2. -1.6 -1.2 -0.8 -0.4 0. 0.4 0.8 1.2 1.6 2. 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 -2. -1.6 -1.2 -0.8 -0.4 0. 0.4 0.8 1.2 1.6 2. 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 -2. -1.6 -1.2 -0.8 -0.4 0. 0.4 0.8 1.2 1.6 2. 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 -2. -1.6 -1.2 -0.8 -0.4 0. 0.4 0.8 1.2 1.6 2. 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 -2. -1.6 -1.2 -0.8 -0.4 0. 0.4 0.8 1.2 1.6 2. 2. ] 2. ] 2. ] 2. ] 2. ] 2. ] 2. ] 2. ] 2. ] 2. ] 2. ]] -2. ] -1.6] -1.2] -0.8] -0.4] 0. ] 0.4] 0.8] 1.2] 1.6] 2. ]] Apply f() to each of the points in the mesh: In [5]: Z = f_vec(X0, X1) print("Z:\n", Z) Z: [[8. 6.56 5.44 4.64 4.16 4. 4.16 4.64 5.44 6.56 8. ] [6.56 5.12 4. 3.2 2.72 2.56 2.72 3.2 4. 5.12 6.56] . [5.44 [4.64 [4.16 [4. [4.16 [4.64 [5.44 [6.56 [8. 4. 3.2 2.72 2.56 2.72 3.2 4. 5.12 6.56 2.88 2.08 1.6 1.44 1.6 2.08 2.88 4. 5.44 2.08 1.28 0.8 0.64 0.8 1.28 2.08 3.2 4.64 1.6 0.8 0.32 0.16 0.32 0.8 1.6 2.72 4.16 1.44 0.64 0.16 0. 0.16 0.64 1.44 2.56 4. 1.6 0.8 0.32 0.16 0.32 0.8 1.6 2.72 4.16 2.08 1.28 0.8 0.64 0.8 1.28 2.08 3.2 4.64 2.88 2.08 1.6 1.44 1.6 2.08 2.88 4. 5.44 4. 3.2 2.72 2.56 2.72 3.2 4. 5.12 6.56 5.44] 4.64] 4.16] 4. ] 4.16] 4.64] 5.44] 6.56] 8. ]] Plot Z[:, :] as a three-dimensional surface: In [6]: from mpl_toolkits.mplot3d import Axes3D from matplotlib.pyplot import figure, xlabel, ylabel %matplotlib inline fig = figure(figsize=(10, 5)) ax3d = fig.add_subplot(121, projection='3d') ax3d.plot_wireframe(X0, X1, Z) xlabel('x0') ylabel('x1') ax2d = fig.add_subplot(122) cp = ax2d.contour(X0, X1, Z) xlabel('x0') ylabel('x1') Out[6]: Text(0,0.5,'x1') The gradient of with respect to is Exercise 1 (1 point). Implement a function, grad_f(x0, x1), that implements the gradient gradient for this shown above. It should return a pair of values since the has two components. In [7]: def grad_f(x0, x1): ### BEGIN SOLUTION return (2*x0, 2*x1) ### END SOLUTION In [8]: # Test cell: `grad_f_test` grad_f_vec = vectorize(grad_f) z = randn(5) gx, gy = grad_f_vec(z, -z) assert all(isclose(gx*0.5, z)) and all(isclose(gy*(-0.5), z)), "Your function might have a bug..." print("\n(Passed!)") (Passed!) Visualizing the gradient Let's generate and plot : In [9]: dX0, dX1 = grad_f(X0, X1) from matplotlib.pyplot import contour, quiver, axis cp = contour (X0, X1, Z) quiver(X0, X1, dX0, dX1, scale=40, headwidth=5) xlabel('x0') ylabel('x1') axis('square') Out[9]: (-2.0, 2.0, -2.0, 2.0) Fin! If you've gotten this far without errors, your notebook is ready to submit. Part 2: Algorithms for the linear least squares problem Recall the linear regression problem: given a data matrix, , and responses , we wish to determine the model parameters that minimizes . This problem is also known as the linear least squares problem. Numpy has a function, np.linalg.lstsq() (https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html), that will compute linear least squares solutions for you. However, the purpose of this notebook is to give you a sense of how lstsq() works. So, instead of using it as a black box, you'll study how one might implement lstsq() using two different numerical algorithms. You may rightly ask, why bother with such details? Here are three reasons it's worth looking more closely. 1. It's helpful to have some deeper intuition for how one formalizes a mathematical problem and derives a computational solution, in case you ever encounter a problem that does not exactly fit what a canned library can do for you. 2. If you have ever used a statistical analysis package, it's likely you have encountered "strange" numerical errors or warnings. Knowing how problems are derived can help you understand what might have gone wrong. We will see an example below. 3. Because data analysis is quickly evolving, it's likely that new problems and new models will not exactly fit the template of existing models. Therefore, it's possible you will need to derive a new model or know how to talk to someone who can derive one for you. Implementation note. In this notebook, we ask that you use the following convention: any column vector should be explicit. That means its shape should have two dimensions where the column dimension equals one (1). Exercise 0 (ungraded). Inspect the following code cell and make sure you understand the difference between two conventions for storing a vector, namely, as a one-dimensional array versus as a two-dimensional array (matrix) where the number of columns equals one (1). When you are asked to produce a vector, we will generally ask you to follow the second convention (z_colvec). In [1]: import numpy as np # By default, Numpy vectors constructed from a list are 1-D # and have neither a row nor a column orientation. z_array = np.array([1.0, 2.0, 3.0]) # By contrast, we want you to ensure your vectors are # column vectors. z_colvec = np.reshape(z_array, (len(z_array), 1)) print("`z_array`:\n\n", z_array, "\n\n==> shape:", z_array.shape) print("\n") print("`z_colvec`:\n\n", z_colvec, "\n\n==> shape:", z_colvec.shape) `z_array`: [1. 2. 3.] ==> shape: (3,) `z_colvec`: [[1.] [2.] [3.]] ==> shape: (3, 1) Before beginning, run this code cell to load some of the key modules you'll need. In [2]: # Data import import import import and computation numpy as np scipy as sp scipy.linalg pandas as pd # Viz from IPython.display import display, Math from matplotlib.pyplot import figure, subplot, xlim, ylim from matplotlib.pyplot import scatter, axis, xlabel, ylabel, title, plot %matplotlib inline # Some functions we'll use later to display results def show_cond_fancy(x, name, opt=''): """Display a condition number in 'fancy' format (using LaTeX).""" def sci_to_latex(x, fmt='{:.2e}'): s_raw = fmt.format(x) s, e = s_raw.split('e') return s + r'\times 10^{{{}}}'.format(int(e)) from IPython.display import Math x_s = sci_to_latex(x) display(Math(r'\kappa({}){} \approx {}'.format(name, opt, x_s))) def show_2vecs_tibble(x, y, xname='x', yname='y', error=False): """Display two column vectors side-by-side in a tibble.""" assert type(x) is np.ndarray and x.ndim >= 2 and x.shape[1] == 1 assert type(y) is np.ndarray and y.ndim >= 2 and y.shape[1] == 1 assert x.shape == y.shape x_df = pd.DataFrame(x, columns=[xname]) y_df = pd.DataFrame(y, columns=[yname]) df = pd.concat([x_df, y_df], axis=1) if error: df['error'] = x - y display(df) # Display (X, y) problem as a tibble def make_data_tibble(X, y=None): df = pd.DataFrame(X, columns=['x_{}'.format(i) for i in range(X.shape[1])]) if y is not None: y_df = pd.DataFrame(y, columns=['y']) df = pd.concat([y_df, df], axis=1) return df # From: https://stackoverflow.com/questions/17129290/numpy-2d-and-1d-array-to-latex-bmatrix def nparray_to_bmatrix(a): """Returns a LaTeX bmatrix""" assert len(a.shape) <= 2, 'bmatrix can at most display two dimensions' lines = str(a).replace('[', '').replace(']', '').splitlines() rv = [r'\begin{bmatrix}'] rv += [' ' + ' & '.join(l.split()) + r'\\' for l in lines] rv += [r'\end{bmatrix}'] return '\n'.join(rv) # Stash this function for later: SAVE_LSTSQ = np.linalg.lstsq # You may ignore this line, which some test cells will use Notation and review Here is a quick summary of how we can formulate and approach the linear regression problem. For a more detailed derivation, see these accompanying notes (https://courses.edx.org/asset-v1:GTx+CSE6040x+1T2018+type@asset+block@nb12-notes-linreg.html). Your data consists of variables, observations and variables. One of these variables is the response variable, , which you want to predict from the other . You wish to fit a linear model of the following form to these data, where is the set of unknown coefficients. Your modeling task is to choose values for these coefficients that "best fit" the data. If we further define a set of dummy variables, , associated with the parameter, then the model can be written more compactly in matrix notation as where we will refer to as the (input) data matrix. Visually, you can also arrange the observations into a tibble like this one: y x0 x1 xn-1 xn 1.0 1.0 1.0 1.0 1.0 This tibble includes an extra column (variable), , whose entries are all equal to 1.0. Synthetic problem generator. For the exercises in this notebook, we will generate synthetic data. The function, gen_problem(m, n), will return a triple, X, y, theta, which are an m x (n+1) data matrix X, a response vector y, and the "true" model parameters theta. We will then run two different numerical algorithms that estimate theta from X and y, and see how their answers compare against the true value. Note 1. The problem generator constructs the data matrix such that each entry is . This structure makes it an instance of a Vandermonde matrix (https://en.wikipedia.org/wiki/Vandermonde_matrix), which arises when fitting a polynomial to data. The "true" parameter vector is set to all ones, and computed simply by summing the rows. Note 2. Although our usual convention is to make the last column all ones, the Vandermonde matrix has its first column set to all ones. This ordering is not important in this problem, but it does mean one would interpret convention. as the intercept rather than , which will be our usual In [3]: def gen_problem(m, n): from numpy import arange, tile, cumprod, insert, ones # 1 + x + x^2 + ... + x^n, x = 0:m X = np.empty((m, n+1)) x_col = arange(m).reshape((m, 1)) # 0, 1, 2, ..., m-1 X[:, 0] = 1.0 X[:, 1:] = tile(x_col, reps=(1, n)) X[:, 1:] = cumprod(X[:, 1:], axis=1) theta = ones((n+1, 1)) y = np.sum(X, axis=1).reshape((m, 1)) return X, y, theta print("Sample generated problem:") m, n = 10, 2 X, y, theta_true = gen_problem(m, n) display(Math(r'X = {}, \quad y = {} \quad \implies \quad \theta^* = {}'.format(nparray_to_bmatrix(X), nparray_to_bmatrix(y), nparray_to_bmatrix(theta_true)))) Sample generated problem: We are interested primarily in overdetermined systems, meaning has more rows than columns, i.e., , as shown above. That's because typically we have more observations (data points, or rows) than predictors (variables or columns). For such problems, there is generally no unique solution. Therefore, to identify some solution, we need to ask for the "best" fit and say what we mean by "best." For linear regression, the usual definition of "best" is minimizing the sum-of-squared residual error: Solving this minimization problem is equivalent to solving a special system known as the normal equations, So, our computational task is to solve this problem. Algorithm 1: Direct solution of the normal equations The preceding calculation immediately suggests the following algorithm to estimate 1. Form 2. Form 3. Solve . Given and : . This object is sometimes called the Gram matrix (https://en.wikipedia.org/wiki/Gramian_matrix) or Gramian of . for . . But, is this a "good" algorithm? There are at least three dimensions along which we might answer this question. 1. Is it accurate enough? 2. Is it fast enough? 3. Is it memory-efficient enough? Let's examine these questions by experiment. Exercise 1 (3 points). Implement a function, solve_neq(X, y) that implements Algorithm 1. It should return a Numpy vector containing the model parameter estimates. Recall the steps of the algorithm as previously outlined: 1. Form the Gramian of 2. Form 3. Solve . for , . . Your algorithm should carry out these steps. For the third step, use Scipy's routine, scipy.linalg.solve() (https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#solving-linear-system). It has an option that allows you to indicate that positive definite, which will be true of for our synthetic problem. is symmetric and The code cell will run your function to compute a set of parameter estimates. It will store these in a variable named theta_neq, which we will refer to later. In [4]: def solve_neq(X, y): ### BEGIN SOLUTION C = X.T.dot(X) b = X.T.dot(y) theta_est = sp.linalg.solve(C, b, sym_pos=True) return theta_est ### END SOLUTION theta_neq = solve_neq(X, y) print("Your implementation's solution versus the true solution:") show_2vecs_tibble(theta_neq, theta_true, xname='theta_neq', yname='theta_true', error=True) Your implementation's solution versus the true solution: theta_neq theta_true error 0 1.0 1.0 -2.065015e-14 1 1.0 1.0 4.218847e-15 2 1.0 1.0 0.000000e+00 In [5]: # Test cell: `solve_neq_test` try: del np.linalg.lstsq solve_neq(X, y) except NameError as n: if re.findall('lstsq', n.args[0]): print("*** Double-check that you did not try to use `lstsq()`. ***") raise n except AttributeError as a: if re.findall('lstsq', a.args[0]): print("*** Double-check that you did not try to use `lstsq()`. ***") raise a finally: np.linalg.lstsq = SAVE_LSTSQ assert type(theta_neq) is np.ndarray, "`theta_neq` should be a Numpy array, but isn't." assert theta_neq.shape == (n+1, 1), "`theta_neq.shape` is {} instead of {}.".format(theta_neq.shape, (n+1, 1)) assert (np.abs(theta_neq - theta_true) <= 1e-12).all(), \ "Your `theta_neq` does not match the true solution, `theta_true`." print("\n(Passed!)") (Passed!) Exercise 2 (1 point). Write a function to calculate the residual norm, Although we are minimizing , for this exercise your function should return . . In [6]: def calc_residual_norm(X, y, theta): ### BEGIN SOLUTION from numpy.linalg import norm return norm(X.dot(theta) - y, ord=2) ### END SOLUTION r_norm_neq = calc_residual_norm(X, y, theta_neq) print("\nThe squared residual norm:", r_norm_neq) The squared residual norm: 3.766978705903275e-14 In [7]: # Test cell: `calc_residual_norm_test` r_norm_neq = calc_residual_norm(X, y, theta_neq) assert 1e-16 <= np.abs(r_norm_neq) <= 1e-12 print ("\n(Passed.)") (Passed.) Sources of error We said before that one question we should ask about our algorithm is whether it is "accurate enough." But what does that mean? Exercise 3 (ungraded). For any modeling problem, there will be several sources of error. Describe at least three such sources. Answer. Here are some possibilities. 1. There will be errors in the inputs. That is, the data itself may only represent measurements of a certain accuracy. 2. There will be errors in the model. That is, the model is only an approximation of the underlying phenomena. 3. There will be errors in the algorithm. That is, you may implement an algorithm that can only approximately estimate the parameters of the model. 4. There will be roundoff errors. Recall that floating-point arithmetic necessarily represents all values finitely, which means you may lose accuracy every time you do an arithmetic operation. Perturbations. One way to understand error in a numerical computation is to consider how sensitive the computed solution is to perturbations to the input. That is, suppose we change by an amount . We can then ask by how much the computed model parameters change. If they change by a lot, then our method for computing them may be overly sensitive to perturbations. Instead, we might prefer one method over another one that is more sensitive to changes. Let's see how Algorithm 1 fares under small perturbations. But first, we'll need a method to generate a random perturbation of a certain maximize size. Exercise 4 (2 points). Implement a function that returns an matrix whose entries are uniformly randomly distributed in the interval, for a given value of . Hint: Check out Numpy's module for generating (pseudo)random numbers: numpy.random (https://docs.scipy.org/doc/numpy/reference/routines.random.html) In [8]: def random_mat (m, n, eps): ### BEGIN SOLUTION return np.random.random((m, n)) * eps ### END SOLUTION print(random_mat(3, 2, 1e-3)) [[0.00059339 0.00057468] [0.00092299 0.00063703] [0.0003482 0.00091574]] In [9]: # Test cell: `rand_eps_test` Z = random_mat (5, 3, 1e-2) assert Z.shape == (5, 3) assert ((Z >= 0) & (Z <= 1e-2)).all() print("\n(Passed.)") (Passed.) Exercise 5 (2 points). Use your random_mat() function to write another function, perturb_system(X, y, eps), that creates two "perturbations" to the system defined by X and y. 1. Let be the first perturbation. It should have the same dimensions as , and its entries should lie in the interval . The value of eps. 2. The second is , a small perturbation to the response variable, . Its entries should also lie in the same interval, Your function should return a perturbed system, and , , as a pair. In [10]: def perturb_system(X, y, eps): ### BEGIN SOLUTION Delta_X = random_mat(X.shape[0], X.shape[1], 2*eps) - eps Delta_y = random_mat(y.shape[0], y.shape[1], 2*eps) - eps return X + Delta_X, y + Delta_y ### END SOLUTION EPSILON = 0.1 X_perturbed, y_perturbed = perturb_system(X, y, EPSILON) Delta_X = X_perturbed - X Delta_y = y_perturbed - y display(Math(r'\Delta X = {}, \quad \Delta y = {}'.format(nparray_to_bmatrix(Delta_X[:5, :]), nparray_to_bmatrix(Delta_y[:5])))) In [11]: # Test cell: `delta_X_test` Delta X = X perturbed - X is given by Delta_y = y_perturbed - y assert Delta_X.shape == X.shape, "`Delta_X` has shape {} instead of {}.".format(Delta_X.shape, X.shape) assert (np.abs(Delta_X) <= EPSILON).all(), "The perturbation lies outside the interval, [-{}, {}]".format(EPSILON, EPSILON) assert Delta_y.shape == y.shape, "`Delta_y` has shape {} instead of {}.".format(Delta_y.shape, y.shape) assert (np.abs(Delta_y) <= EPSILON).all(), "The perturbation lies outside the interval, [-{}, {}]".format(EPSILON, EPSILON) print ("\n(Passed.)") (Passed.) Sensitivity of Algorithm 1 Let's now run the following code, which uses your code from above to perform a "sensitivity experiment." In particular, the function run_perturbation_trials() will repeatedly perturb the system and measure the resulting change to the estimated All of the estimated . are stored in an array, Thetas_neq. Each column k of Thetas_neq, or Thetas_neq[:, k], is one of the calculated estimates under a random perturbation of the system. The size of the random perturbation is set, by default, to eps=0.01. Recall that our synthetic problem consists of numerical values that are all greater than or equal to one, so this perturbation may be regarded as fairly small. In [12]: def run_perturbation_trials(solver, X, y, eps=0.01, trials=100): Thetas = np.zeros((X.shape[1], trials)) # Store all computed thetas for t in range(trials): X_p, y_p = perturb_system(X, y, eps) Thetas[:, t:t+1] = solver(X_p, y_p) return Thetas Thetas_neq = run_perturbation_trials(solve_neq, X, y) print("Unperturbed solution:") print(theta_neq) print("First few perturbed solutions (columns):") print(Thetas_neq[:, :5]) Unperturbed solution: [[1.] [1.] [1.]] First few perturbed solutions (columns): [[0.99028632 0.99231283 1.02243806 1.00139275 1.01358837] [1.00529885 1.00407045 0.99118553 0.99889269 0.9952955 ] [0.99954928 0.99960384 1.00074977 1.00007595 1.0003454 ]] Here is a quick plot of the that shows two coordinates of the true parameters (red star), compared to all perturbed estimates (blue points). We would have more confidence in the algorithm's computed solutions if it did not appear to be too sensitive to changes in the input. Since may have more than two coordinates, the code below shows the first two coordinates. In [13]: # # # # Makes a 2-D scatter plot of given theta values. If the thetas have more than two dimensions, only the first and last are displayed by default. (Override by setting ax and ay.) def scatter_thetas(Thetas, theta_true=None, ax=0, ay=-1, xylim=None, title=None): import matplotlib.pyplot as plt assert type(Thetas) is np.ndarray and Thetas.shape[0] >= 2 scatter(Thetas[ax, :], Thetas[ay, :]) xlabel('{}-coordinate'.format(ax if ax >= 0 else Thetas.shape[0]+ax)) ylabel('{}-coordinate'.format(ay if ay >= 0 else Thetas.shape[0]+ay)) if xylim is not None: axis(xylim) else: axis('equal') if theta_true is not None: assert type(theta_true) is np.ndarray and theta_true.shape[0] >= 2 and theta_true.shape[1] == 1 scatter(theta_true[ax], theta_true[ay], marker='*', color='red', s=15**2) if title is not None: plt.title(title) def calc_lims(x, buffer=0.1): xmin = x min() xmin = x.min() xmax = x.max() dx = (xmax - xmin) * buffer return xmin-dx, xmax+dx scatter_thetas(Thetas_neq, theta_true=theta_true, ax=0, ay=2) axis('square'); You should observe that the change in the estimates are of the same order as the perturbation. So for this example system, the algorithm seems reliable enough. Stress-testing Algorithm 1 This experiment suggests all is fine. But what should we expect to happen? We've prepared another notebook (https://courses.edx.org/asset-v1:GTx+CSE6040x+1T2018+type@asset+block@nb12-notes-cond.html) that walks you through an analysis of solving linear systems. It turns out you can estimate how hard it is to solve a linear system using a measure called the condition number. We can denote the condition number of solving a system by is. where is the matrix. The larger this number is, the more sensitive the problem In Numpy, there is a condition number estimator that will tell us approximately what the condition number is for a given matrix. Let's compare : and In [14]: cond_X = np.linalg.cond(X) cond_XTX = np.linalg.cond(X.T.dot(X)) assert 1. <= cond_X <= 3e3 assert 1. <= cond_XTX <= 6e6 show_cond_fancy(cond_X, 'X') show_cond_fancy(cond_XTX, 'X^T X') show_cond_fancy(cond_X**2, 'X', opt='^2') Ill-conditioning. As it happens, is roughly the square of . So, by forming explicitly and then trying to solve a system based on it, we make the problem more difficult. Indeed, if the problem is ill-conditioned enough, this algorithm based on directly constructing the normal equations will produce very different results even under small changes, and we call the algorithm unstable. In this particular example, the condition numbers are not very "big." You would be more concerned if the condition numbers were close to machine epsilon. In double-precision, recall that , so the values shown above is nothing to be worried about. , where is But what if we had a "hard" problem, that is, one whose condition number is large? The synthetic data generator allows us to create such a problem by making the problem bigger. Let's try that next. In [15]: # Generate a "hard" problem m_hard, n_hard = 100, 6 X_hard, y_hard, theta_hard_true = gen_problem(m_hard, n_hard) df_hard = make_data_tibble(X_hard, y_hard) print("First few rows of data:") df_hard.head() print("True parameter estimates:\n{}".format(theta_hard_true)) cond_X_hard = np.linalg.cond(X_hard) cond_XTX_hard = np.linalg.cond(X_hard.T.dot(X_hard)) name_X_hard = 'X_h' show_cond_fancy(cond_X_hard, name_X_hard) show_cond_fancy(cond_XTX_hard, '{}^T {}'.format(name_X_hard, name_X_hard)) First few rows of data: True parameter estimates: [[1.] [1.] [1.] [1.] [1.] [1.] [1.]] These condition numbers are much larger. So, let's run the same sensitivity experiment as before, and see how the estimate varies for the hard problem. How does it compare to the well-conditioned case? In [16]: Thetas_hard_neq = run_perturbation_trials(solve_neq, X_hard, y_hard) scatter_thetas(Thetas_hard_neq, theta_true=theta_hard_true, ax=0, ay=2) print("Residual norm for one of the trials:") theta_hard_neq_example = np.random.randint(Thetas_hard_neq.shape[1]) calc_residual_norm(X_hard, y_hard, theta_hard_neq_example) Residual norm for one of the trials: Out[16]: 72708670496575.33 Observe that the computed estimates can be relatively far from the true value, even getting the sign completely wrong in the case of the . Algorithm 2: QR decomposition A different method for solving an overdetermined systems is to use a tool from linear algebra known as the QR decomposition (https://en.wikipedia.org/wiki/QR_decomposition). Here is how we can use QR. If has linearly independent columns, then we would first factor the matrix into the product orthogonal matrix and is an invertible upper-triangular matrix. (These dimensions assume .) That , the identity matrix; being upper-triangular means all of its entries below the main diagonal are zero. Next, observe that the normal equations can be transformed if we substitute , where is orthogonal means that : Lastly, because is triangular, solving a system is "easy" using (backward) substitution. Consider the following (http://www.purplemath.com/modules/systlin6.htm)): example (taken from here is an Because it is upper-triangular, you can see right away that . Lastly, So, to summarize, a different algorithm to solve 1. Compute . Then, going to the equation above it, . using QR would look like the following: . 2. Form the modified right-hand side, 3. Use back-substitution to solve . . Conditioning. What about the sensitivity of this algorithm? Given , we only need to solve linear systems involving . Therefore, it's that will determine the stability of the algorithm. So if is comparable to , then the algorithm should be as stable as one can expect any algorithm to be. Exercise 6 (1 point). Use numpy.linalg.qr() (https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.qr.html) to compute the QR decomposition of (precomputed above as the variable, X). Store the In [17]: and factors in two variables named Q and R. print(X[:5], "\n ...\n") ### BEGIN SOLUTION Q, R = np.linalg.qr(X) ### END SOLUTION # Print the dimensions of your result print("Q:", Q.shape, "\n") print("R:", R.shape, "==") print(R) [[ 1. [ 1. [ 1. [ 1. [ 1. ... 0. 0.] 1. 1.] 2. 4.] 3. 9.] 4. 16.]] Q: (10, 3) R: (3, 3) == [[ -3.16227766 -14.23024947 -90.12491331] [ 0. 9.08295106 81.74655956] [ 0. 0. 22.97825059]] In [18]: # Test cell: `qr_test` assert type(Q) is np.ndarray, "`Q` is not a Numpy array but should be." assert type(R) is np.ndarray, "`R` is not a Numpy array but should be." assert Q.shape == (m, n+1), "`Q` has the wrong shape: it's {} rather than {}.".format(Q.shape, (m, n+1)) assert R.shape == (n+1, n+1), "`R` has the wrong shape: it's {} rather than {}.".format(R.shape, (m, n+1)) for i in range(R.shape[0]): for j in range(i): assert np.isclose(R[i][j], 0.0), "R[{}][{}] == {} instead of 0!".format(i, j, R[i][j]) QTQ = Q.T.dot(Q) assert np.isclose(QTQ, np.eye(Q.shape[1])).all(), "Q^T Q is not nearly the identity matrix, as it should be." assert np.isclose(X, Q.dot(R)).all(), "QR is not sufficiently close in values to X!" print("\n(Passed!)") (Passed!) Condition number of In [19]: . Let's check the condition number of cond_R = np.linalg.cond(R) show_cond_fancy(cond_X, 'X') show_cond_fancy(cond_XTX, 'X^T X') show_cond_fancy(cond_R, 'R') empirically, to verify that it is comparable to . Exercise 7 (3 points). Implement a function, solve_qr(X, y), which uses the QR-based algorithm to estimate . To solve the triangular system, use Scipy's specialized function, available as sp.linalg.solve_triangular() (https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_triangular.html). In [20]: import scipy.linalg def solve_qr(X, y): ### BEGIN SOLUTION Q, R = np.linalg.qr(X) b = Q.T.dot(y) theta = sp.linalg.solve_triangular(R, b) # Solves R u = b return theta ### END SOLUTION theta_qr = solve_qr(X, y) print("Comparing your QR solution to the true solution:") show_2vecs_tibble(theta_qr, theta_true, xname='theta_qr', yname='theta_true', error=True) print("Residual norm:") calc_residual_norm(X, y, theta_qr) Comparing your QR solution to the true solution: theta_qr theta_true error 0 1.0 1.0 -1.154632e-14 1 1.0 1.0 3.552714e-15 2 1.0 1.0 -2.220446e-16 Residual norm: Out[20]: 1.5434895314732317e-14 In [21]: # Test cell: `solve_qr_test` import re try: del np.linalg.lstsq solve_qr(X, y) except NameError as n: if re.findall('lstsq', n.args[0]): print("*** Double-check that you did not try to use `lstsq()`. ***") raise n except AttributeError as a: if re.findall('lstsq', a.args[0]): print("*** Double-check that you did not try to use `lstsq()`. ***") raise a finally: np.linalg.lstsq = SAVE_LSTSQ assert np.isclose(theta_qr, theta_true).all(), "Your QR-based solution should be closer to the true solution." print("\n(Passed!)") (Passed!) Is QR more stable? Let's run the same perturbation experiments on the "hard" regression problem and see the result. In [22]: Thetas_hard_qr = run_perturbation_trials(solve_qr, X_hard, y_hard) # Plot side-by-side against normal equations method def compare_scatter_thetas(T0, title0, T1, title1, ax=0, ay=1, **kwargs): xmin, xmax = calc_lims(np.array([Thetas_hard_neq[ax, :], Thetas_hard_qr[ax, :]])) ymin, ymax = calc_lims(np.array([Thetas_hard_neq[ay, :], Thetas_hard_qr[ay, :]])) xylim = [xmin, xmax, ymin, ymax] figure(figsize=(12, 4)) subplot(1, 2, 1) scatter thetas(T0 title title0 ax ax ay ay xylim xylim **kwargs) scatter_thetas(T0, title=title0, ax=ax, ay=ay, xylim=xylim, **kwargs) subplot(1, 2 ,2) scatter_thetas(T1, title=title1, ax=ax, ay=ay, xylim=xylim, **kwargs) compare_scatter_thetas(Thetas_hard_neq, 'Normal equations', Thetas_hard_qr, 'QR', ax=0, ay=-1, theta_true=theta_hard_true) print("Sample estimate for one of the trials:") theta_hard_neq_example = Thetas_hard_neq[:, np.random.randint(Thetas_hard_neq.shape[1])] theta_hard_qr_example = Thetas_hard_qr[:, np.random.randint(Thetas_hard_qr.shape[1])] msg = "- {}-based method: theta^T =\n\t{}" print(msg.format("Gramian", theta_hard_neq_example.T)) print(msg.format("QR", theta_hard_qr_example.T)) Sample estimate for one of the trials: - Gramian-based method: theta^T = [0.64091071 1.16266565 0.98363847 1.00065254 0.99998782 1.00000011 1. ] - QR-based method: theta^T = [0.98958625 0.99839116 1.00028587 0.99998727 1.00000024 1. 1. ] You should observe that the QR-based method does, indeed, produce estimates much closer to the true value despite the problem's high condition number. Performance tradeoff. Although QR produces more reliable results, there can be a performance tradeoff, as the following quick test should show. In [23]: print("=== Performance of the normal equations-based algorithm ===") %timeit solve_neq(X_hard, y_hard) print("\n=== Performance of the QR-based algorithm ===") %timeit solve_qr(X_hard, y_hard) === Performance of the normal equations-based algorithm === 87.7 µs ± 1.46 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) === Performance of the QR-based algorithm === 146 µs ± 2.17 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) Summary comment. The intent of this notebook was to help you appreciate some of the reliability and performance issues involved in the design of numerical algorithms. The key ideas are as follows. 1. The conditioning of a problem, which is a measure of how hard it is to solve. For solving linear systems, the condition number of the matrix produced by numpy.linalg.cond() is this measure. 2. For the specific problem of linear least squares, a method based on QR factorization is better than forming the Gram matrix, which squares the condition number. Although we did not discuss the theory of QR better, it exemplifies a common pattern or trick in linear algebraic problems: converting the original problem (or matrix) into one whose structure is easier to solve. 3. There is a potential tradeoff between reliability of an algorithm and speed. Many of these algorithmic and engineering ideas have already been worked out, and are encapsulated in Numpy's built-in lstsq() solver (https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html) or its Scipy equivalent (https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html). So when you need to solve a linear system in a future notebook, you can use this function. Fin! If you've gotten this far without errors, your notebook is ready to submit.