$28
1. Python and Numpy Warm-up Exercises [20 pts]: This part of the homework is intended to help you practice your linear algebra and how to implement linear algebraic and statistical operations in Python using numpy (to which we refer in the code below as np). For each of the problems below, write a method (e.g., problem 1a) that returns the answer for the corresponding problem. In all problems, you may assume that the dimensions of the matrices and/or vectors that are given as input are compatible for the requested mathematical operations. You do not need to perform error-checking. Note 1: In mathematical notation we usually start indices with j = 1. However, in numpy (and many other programming settings), it is more natural to use 0-based array indexing. When answering the questions below, do not worry about “translating” from 1-based to 0-based indexes. For example, if the (i, j)th element of some matrix is requested, you can simply write A[i,j]. Note 2: To represent and manipulate vectors and matrices, please use numpy’s array class (not the matrix class). Note 3: While the difference between a row vector and a column vector is important when doing math, numpy does not care about this difference as long as the array is 1-D. This means, for example, that if you want to compute the inner product between two vectors x and y, you can just write x.dot(y) without needing to transpose the x. If x and y are 2-D arrays, however, then it does matter whether they are row-vectors or column-vectors, and hence you might need to transpose accordingly. (a) Given two matrices A and B, compute and return an expression for A + B. [ 0 pts ] Answer : While it is completely valid to use np.add(A, B), this is unnecessarily verbose; you really should make use of the “syntactic sugar” provided by Python’s/numpy’s operator overloading and just write: A + B. Similarly, you should use the more compact (and arguably more elegant) notation for the rest of the questions as well. (b) Given matrices A, B, and C, compute and return AB − C (i.e., right-multiply matrix A by matrix B, and then subtract C). Use dot or np.dot. [ 1 pts ] (c) Given matrices A, B, and C, return A⊙B+C⊤, where ⊙ represents the element-wise (Hadamard) product and ⊤ represents matrix transpose. In numpy, the element-wise product is obtained simply with *. [ 1 pts ] (d) Given column vectors x and y, compute the inner product of x and y (i.e., x ⊤y). [ 1 pts ] (e) Given matrix A and integer i, return the sum of all the entries in the ith row whose column index is even, i.e., P j:j is even Aij . Do not use a loop, which in Python can be very slow. Instead use the np.sum function. [ 2 pts ] (f) Given matrix A and scalars c, d, compute the arithmetic mean over all entries of A that are between c and d (inclusive). In other words, if S = {(i, j) : c ≤ Aij ≤ d}, then compute 1 |S| P (i,j)∈S Aij . Use np.nonzero along with np.mean. [ 2 pts ] (g) Given an (n×n) matrix A and integer k, return an (n×k) matrix containing the right-eigenvectors of A corresponding to the k eigenvalues of A with the largest absolute value. Use np.linalg.eig. [ 2 pts ] 1 (h) Given a column vector (with n components) x, an integer k, and positive scalars m, s, return an (n × k) matrix, each of whose columns is a sample from multidimensional Gaussian distribution N (x + mz, sI), where z is column vector (with n components) containing all ones and I is the identity matrix. Use either np.random.multivariate normal or np.random.randn. [ 2 pts ] (i) Given a matrix A with n rows, return a matrix that results from randomly permuting the columns (but not the rows) in A. [ 2 pts] (j) Z-scoring: Given a vector x, return a vector y such that each yi = (xi − x)/σ, where x is the mean (use np.mean) of the elements of x and σ is the standard deviation (use np.std). [ 2 pts ] (k) Given an n-vector x and a non-negative integer k, return a n × k matrix consisting of k copies of x. You can use numpy methods such as np.newaxis, np.atleast 2d, and/or np.repeat. [ 2 pts ] (l) Given a k × n matrix X = x (1) . . . x (n) and a k × m matrix Y = y (1) . . . y (m) , compute an n × m matrix D = d11 . . . d1m . . . dn1 . . . dnm consisting of all pairwise L2 distances dij = ∥x (i) − y (j)∥2. In this problem you may not use loops. Instead, you can avail yourself of numpy objects & methods such as np.newaxis, np.atleast 3d, np.repeat, np.swapaxes, etc. (There are various ways of solving it.) Hint: from X (resp. Y), construct a 3-d matrix that contains multiple copies of each of the vectors in X (resp. Y); then subtract these 3-d matrices. [ 3 pts ] 2. Training 2-Layer Linear Neural Networks with Stochastic Gradient Descent [25 pts]: (a) Train an age regressor that analyzes a (48 × 48 = 2304)-pixel grayscale face image and outputs a real number ˆy that estimates how old the person is (in years). The training and testing data are available here: • https://s3.amazonaws.com/jrwprojects/age_regression_Xtr.npy • https://s3.amazonaws.com/jrwprojects/age_regression_ytr.npy • https://s3.amazonaws.com/jrwprojects/age_regression_Xte.npy • https://s3.amazonaws.com/jrwprojects/age_regression_yte.npy Your prediction model g should be a 2-layer linear neural network that computes ˆy = g(x; w) = x ⊤w + b, where w is the vector of weights and b is the bias term. The cost function you should optimize is fMSE(w, b) = 1 2n Xn i=1 (ˆy (i) − y (i) ) 2 where n is the number of examples in the training set Dtr = {(x (1), y(1)), . . . ,(x (n) , y(n) )}, each x (i) ∈ R 2304 and each y (i) ∈ R. To optimize the weights, you should implement stochastic gradient descent (SGD). Note: you must complete this problem using only linear algebraic operations in numpy – you may not use any off-the-shelf linear regression or neural network software, as that would defeat the purpose. There are several different hyperparameters that you will need to optimize: • Mini-batch size ˜n. • Learning rate ϵ. • Number of epochs. In order not to cheat (in the machine learning sense) – and thus overestimate the performance of the network – it is crucial to optimize the hyperparameters only on a validation set. (The training set would also be acceptable but typically leads to worse performance.) To create 2 a validation set, simply set aside a fraction (e.g., 20%) of the age regression Xtr.npy and age regression ytr.npy to be the validation set; the remainder (80%) of these data files will constitute the “actual” training data. While there are fancier strategies (e.g., Bayesian optimization) that can be used for hyperparameter optimization, it’s often effective to just use a grid search over a few values for each hyperparameter. In this problem, you are required to explore systematically (e.g., using nested for loops) at least 2 different values for each hyperparameter. Performance evaluation: Once you have tuned the hyperparameters and optimized the weights and bias term so as to minimize the cost on the validation set, then: (1) stop training the network and (2) evaluate the network on the test set. Report both the training and the test fMSE in the PDF document you submit, as well as the training cost values for the last 10 iterations of gradient descent (just to show that you actually executed your code). 3. Gradient descent: what can go wrong? [30 pts] Please enter your code in a file named homework1 problem3.py, with one Python function (e.g., problem3a) for each subproblem. (a) [10 pts]: The graph below plots the function f(x) that is defined piece-wise as: f(x) = −x 3 : x < −0.1 −3x/100 − 1/500 : −0.1 ≤ x < 3 −(x − 31/10)3 − 23/250 : 3 ≤ x < 5 1083 200 (x − 6)2 − 6183/500 : x >= 5 4 2 0 2 4 6 8 10 0 20 40 60 As you can see, the function has a long nearly flat section (sometimes known as a plateau) just before the minimum.1 Plateaux can cause big problems during optimization. To show this: i. Derive by hand the (piece-wise) function ∇f and implement it in Python/numpy. ii. Use your implementation of ∇f to conduct gradient descent for T = 100 iterations. Always start from an initial x = −3. Try using various learning rates: 1e-3, 1e-2, 1e-1, 1e0, 1e1. Plot f, ∇f, as well as superimposed dots that show the sequence ((x (1), y(1)), . . . ,(x (T) , y(T) )) of gradient descent. Use plt.legend to indicate which scatter plot corresponds to which learning rate. iii. Describe in 1-2 sentences what you observe during gradient descent for the set of learning rates listed above. iv. Find a learning rate ϵ for which gradient descent successfully converges to min f(x), and report ϵ in the PDF file. (b) [8 pts]: Even a convex paraboloid – i.e., a parabola in multiple dimensions that has only one local minimum and no plateaux – can cause problems for “vanilla” SGD (i.e., the kind we’ve learned in 1The ugly constants in f were chosen to give rise to these characteristics while ensuring that it remains differentiable. 3 class so far). Examine the scatter-plot below which shows the sequence ((x (1), y(1)), . . . ,(x (T) , y(T) )) of gradient descent on a convex paraboloid f, starting at x (1) = [1, −3]⊤, where each x ∈ R 2 . The descent produces a zig-zag pattern that takes a long time to converge to the local minimum. 4 2 0 2 4 x1 3 2 1 0 1 2 3 4 x2 i. Speculate how the SGD trajectory would look if the learning rate were made to be very small (e.g., 100x smaller than in the figure above). ii. Let f(x1, x2) = a1(x1 − c1) 2 + a2(x2 − c2) 2 . Pick values for a1, a2, c1, c2 so that – when the scatter-plot shown above is superimposed onto it – the gradient descent is realistic for f. Rather than just guess randomly, consider: why would the zig-zag be stronger in one dimension than another, and how would this be reflected in the function’s constants? Plot a contour graph using plt.contour and superimpose the scatter-plot using plt.scatter. You can find the gradient descent sequence in gradient descent sequence.txt. Note that you are not required to find the exactly correct constants or even to estimate them algorithmically. Rather, you can combine mathematical intuition with some trial-and-error. Your solution should just look visually “pretty close” to get full credit. Note: to ensure proper rendering, use plt.axis(’equal’) right before calling plt.show(). (c) [6 pts]: This problem is inspired by this paper. Consider the function f(x) = 2 3 |x| 3/2 . Derive ∇f, implement gradient descent and plot the descent trajectories ((x (1), y(1)), . . . ,(x (T) , y(T) )) for a variety of learning rates 1e-3, 1e-2, 1e-1 and a variety of starting points. See what trend emerges, and report it in the PDF. (d) [6 pts]: While very (!) unlikely, it is theoretically possible for gradient descent to converge to a local maximum. Give the formula of a function (in my own solution, I used a degree-4 polynomial), a starting point, and a learning rate such that gradient descent converges to a local maximum after 1 descent iteration (i.e., after 1 iteration, it reaches the local maximum, and the gradient is exactly 0). Prove (by deriving the exact values for the descent trajectory) that this is true. You do not need to implement this in code (and, in fact, due to finite-precision floating-point arithmetic, it might not actually converge as intended). Submission: Create a Zip file containing both your Python and PDF files, and then submit on Canvas. If you are working as part of a group, then only one member of your group should submit (but make sure you have already signed up in a pre-allocated team for the homework on Canvas). 4