CS5487 Programming Assignment 1 Regression Antoni Chan Department of Computer Science City University of Hong Kong In this programming assignment you will implement and test some of the regression methods that were introduced in the lecture and in the problem sets. Let f(x, θ) be a function with input x ∈ Rd and parameters θ ∈ RD, and f(x, θ) = φ(x)T θ, (1) where φ(x) : Rd → RD is a feature transformation of x. For example, the K-th order polynomial with input x ∈ R can be expressed as f(x, θ) = K∑ k=0 xkθk = φ(x) T θ, (2) where the feature transformation and parameters are φ(x) = [ 1, x, x2, · · · , xK]T ∈ RK+1, θ = [θ0, · · · , θK]T ∈ RK+1. (3) Our goal is to obtain the best estimate of the function given iid samples D = {(x1, y1), . . . , (xn, yn)}, where yi are the noisy observations of f(xi, θ). We have seen several ways of doing the regression, using different assumed noise models and formulations. For convenience, define the following quan- tities, y = [ y1, · · · , yn ]T , Φ = [ φ(x1), · · · , φ(xn) ] , X = [ x1, · · · , xn ] . (4) Here is a summary of the various regression algorithms we have seen so far: method objective function parameter estimate θˆ prediction f∗ for input x∗ least-squares (LS) ∥∥y − ΦT θ∥∥2 θˆLS = (ΦΦT )−1Φy f∗ = φ(x∗)T θˆ regularized LS (RLS) ∥∥y − ΦT θ∥∥2 + λ ‖θ‖2 θˆRLS = (ΦΦT + λI)−1Φy f∗ = φ(x∗)T θˆ L1-regularized LS (LASSO) ∥∥y − ΦT θ∥∥2 + λ ‖θ‖1 QP solver (see Prob. 3.12) f∗ = φ(x∗)T θˆ robust regression (RR) ∥∥y − ΦT θ∥∥ 1 LP solver (see Prob. 2.10) f∗ = φ(x∗)T θˆ distributions posterior predictive Bayesian regression (BR) θ ∼ N (0, αI) θ|X, y ∼ N (µˆθ, σˆ2θ) f∗|X, y, x∗ ∼ N (µˆ∗, σˆ2∗) y|x, θ ∼ N (f(x, θ), σ2) µˆθ = 1σ2 ΣˆθΦy µˆ∗ = φ(x∗)T µˆθ Σˆθ = ( 1 αI + 1 σ2 ΦΦ T )−1 σˆ2∗ = φ(x∗) T Σˆθφ(x∗) 1 Part 1 Polynomial function In the first problem, you will test these regression methods on a simple dataset, and look at the effects of using different objective functions, regularization terms, and formulations. The MATLAB data file poly data.mat (or poly data *.txt if you are not using MATLAB) contains a set of samples generated from a 5th order polynomial function, with observation noise σ2 = 5. The file contains the following variables: • sampx – sample input values (xi), each entry is an input value. • sampy – sample output values (yi), each entry is an output. • polyx – input values for the true function (also these are the test inputs x∗). • polyy – output values for the true function. • thtrue – the true value of the θ used to generate the data. A plot of the samples (sampx,sampy) and the true polynomial function (polyx, polyy) are shown below. −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −20 −10 0 10 20 30 40 x y samples function The goal is to try to estimate the polynomial function from the samples. (a) Implement the above 5 regression algorithms for the K-th order polynomial given in (2). In the next problem, you will use these regression methods with a different feature transforma- tion φ(x). Hence, it would be better to separate the regression algorithm and the feature transformation in your implementation. (b) For each regression method, use the sample data (sampx, sampy) to estimate the parameters of a 5th order polynomial function. Make a plot of the estimated function using polyx as inputs, along with the sample data. For BR, also plot the standard deviation around the mean. What is the mean-squared error between the learned function outputs and the true function outputs (polyy), averaged over all input values in polyx? For algorithms with hyperparameters, select some values that tend to work well. 2 (c) Repeat (b), but reduce the amount of training data available, by selecting a subset of the samples (e.g., 10%, 25%, 50%, 75%). Plot the estimated functions. Which models are more robust with less data, and which tend to overfit? Make a plot of error versus training size, and comment on any important trends and findings. (You should run multiple trials with different random subsets, and take the average error). (d) Add some outliers output values (e.g., add large numbers to a few values in sampy), and repeat (b). Which methods are robust to the presence of outliers, and which are most sensitive? Why? (e) Repeat (b) but estimate a higher-order polynomial (e.g., 10th order). Which models tend to overfit the data when learning a more complex model, and which models do not? Verify your observations by examining the estimated parameter values. . . . . . . . . . Part 2 A real world regression problem – counting people Now we will consider a real world regression problem – estimating the number of people in an image. Below is an example image of a sidewalk with people. From each image, we extract a feature vector xi ∈ R9, which is based on some properties of the crowd (e.g., the area covered by the crowd, the perimeter of the crowd, etc). We also have the number of people in the image, yi ∈ R. The goal is to learn a regression function between this feature vector and the number of people in the image. (Let’s ignore the fact that the count is actually a non-negative integer, and just apply our standard regression methods). The data can be found in the count data.mat MATLAB file (or count data *.txt if you are not using MATLAB). The mat file contains the following variables: • trainx – training inputs. Each column is a 9-dimensional feature vector. • trainy – training outputs (mean-subtracted counts for each input vector). • testx – test inputs. • testy – true outputs for the test inputs. Note that the output values are actually the mean-subtracted counts, so you will see some negative outputs. 3 (a) Let’s first look at using the features directly, i.e., set φ(x) = x. Use the training set (trainx, trainy), estimate a function with some of the regression algorithms above. Use the test set inputs testx to predict the output, and compare the predictions to the true output testy. (You can round the function outputs so they are counting numbers). Calculate the mean-absolute error and mean-squared error. Which method works the best? Plot the test predictions and the true counts, and discuss any interesting findings. (b) Now try some other feature transformations. For example, you can create a simple 2nd order polynomial as φ(x) = [ x1, · · · , x9, x21, · · · , x29 ]T . This could be expanded to include more cross- terms xixj . Try some other non-linear transformation of the features. Can you improve your results from (a)? . . . . . . . . . Part 3 Estimating hyperparameters (optional) So far you probably have been setting the hyperparameters of each method (e.g., λ) by hand to get the best performance on the test set. This is not a fair way to do an experiment, since it amounts to optimizing to the test set, and doesn’t reflect whether the model can generalize well to unseen data. In a real scenario, you wouldn’t have the true outputs for the test set! Here are two common ways to set the hyperparameters using only the training set: 1. N-fold cross-validation: With this method, the original training set is split into a new training set and a validation set. The regression function is estimated from the new (smaller) training set, and tested on the validation set, for different hyperparameter values. This is repeated for different partitions of the original training set. Formally, separate the training set D into N equal partitions, {D1, · · · ,DN}. For each partition i and candidate hyperparameter value λj : 1) train a regression function using λj on D\Di (all the data except for Di); 2) test the function on Di; 3) record the error ei,j . Finally, sum the error over the N partitions, Ej = ∑N i=1 ei,j . Select the λj with the lowest overall error, and train the final regression function on the full training set D. (usually people use between N = 3 to N = 10 folds, depending on the size of the
training set). 2. maximum marginal likelihood: For Bayesian regression, another way to select the hyper- parameters is to maximize the marginal likelihood of the training data (see Problem 3.11 for details), αˆ, σˆ2 = argmax α,σ2 p(y|X,σ2, α) = argmax α,σ2 log p(y|X,σ2, α) (5) where the marginal log-likelihood is log p(y|X,σ2, α) = log ∫ p(y|θ,X, σ2)p(θ|α)dθ (6) = −D 2 logα− n 2 log σ2 − 1 2σ2 ∥∥y − ΦT µˆθ∥∥2 − 1 2α ‖µˆθ‖2 − 1 2 log ∣∣∣Σˆ−1θ ∣∣∣− n2 log(2pi). (7) The procedure is as follows: for each candidate hyperparameter pair {αj , σ2j }, calculate the marginal likelihood in (7). Select the hyperparameter pair that results in the largest marginal likelihood. 4 (a) For Problem 2, use cross-validation or maximum marginal likelihood to select the hyperparam- eters using only the training data. Does the accuracy increase or decrease? How close were the automatically selected parameters to your hand-selected ones? . . . . . . . . . • What to hand in – You need to turn in the following things: 1. Answers, plots, analysis, discussion, etc. for the above problems. 2. Source code files. You must submit your programming assignment using the Canvas website. Go to “Assign- ments” ⇒ “Programming Assignment 1”. • Plagiarism – You should implement each regression algorithm using your own code. Do not use someone else’s implementation of the regression algorithms (including MATLAB’s implementation). It is okay to use common library components and functions, e.g. standard matrix operations (e.g., inv, eig), general purpose optimization toolboxes (e.g., quadprog, linprog), etc. If you use any special toolboxes or libraries, make sure to note them in your report. • Grading – The marks for this assignment will be distributed as follows: – 10% – Implementation of the regression algorithms (1a). – 40% – Results testing the regression functions (1b, 1c, 1d, 1e). – 20% – Results on counting people (2a). – 10% – Trying out other features for counting people (2b). – 20% – Quality of the written report. More points given to insightful observations and analysis. Note: if you really cannot implement the algorithms correctly, it is okay to use 3rd party software. You will not get marks for implementation, but you can still get marks for presenting the results. • Optional Problems – The optional problem 3 will not be graded. 5
