1. Introduction

In this exercise, you will implement regularized linear regression and use it to study models with different bias-variance properties. To get started with the exercise, you will need to download the starter code and unzip its contents to the directory where you wish to complete the exercise. You will be required to modify several parts of the code, by following the instructions in this assignment. Files needed for this practical include:


2. Regularized linear regression

In the first half of the exercise, you will implement regularized linear regression to predict the amount of product (monoclonal antibodies) flowing out of a biorreactor using the change of culture medium level in it. In the next half, you will go through some diagnostics of debugging learning algorithms and examine the effects of bias vs variance.


2.1. Visualizing the data

We will begin by visualizing the dataset containing historical records on the change in the culture medium level, \(x\), and the amount of product flowing out of the bioreactor, \(y\).

This dataset is divided into three parts: - A training set that your model will learn on: \(X\), \(y\) - A cross validation set for determining the regularization parameter: \(Xval\), \(yval\) - A test set for evaluating performance. These are “unseen” examples which your model did not use during training: \(Xtest\), \(ytest\)

The next piece of code will plot the training set. In the following parts, you will implement linear regression and use that to fit a straight line to the data and plot learning curves. Following that, you will implement polynomial regression to find a better fit to the data.

load("data1.Rdata")  # load data into R
# You will now have X, y, Xval, yval, Xtest, ytest in your environment

m <- length(y)  # Number of training examples

plot(X, y, xlab = "Change in culture medium level", ylab = "Product flowing out of the bioreactor",
     col = "blue", pch = 19, ylim = c(-5, 40))


2.2. Regularized linear regression cost function

Recall that regularized linear regression has the following cost function:

\[ J(\theta) = \frac{1}{2m}\Big(\sum\limits_{i=1}^m(h_\theta(x^{(i)})-y^{(i)})^2\Big) + \frac{\lambda}{2m}\Big(\sum\limits_{j=1}^n\theta^2_j\Big) \]

where \(\theta\) is a regularization parameter which controls the degree of regu- larization (thus, help preventing overfitting). The regularization term puts a penalty on the overall cost \(J\). As the magnitudes of the model parameters \(\theta_j\) increase, the penalty increases as well. Note that you should not regularize the \(\theta_0\) term. (Remember that in our script, the \(\theta_0\) term is represented as theta[1] since indexing in R starts from 1). You should now complete the code in the file LinearRCostFunctionReg.R Your task is to write a function to calculate the regularized linear regression cost function. If possible, try to vectorize your code and avoid writing loops. When you are finished, the following code will run your cost function using theta initialized at [1, 1]. You should expect to see an output of approximately 303.993.

source("LinearRCostFunctionReg.R")
theta <- c(1, 1) 
J <- LinearRCostFunctionReg(theta, cbind(1, X), y, 1)


2.3. Regularized linear regression gradient

Correspondingly, the partial derivative of regularized linear regression’s cost for \(\theta_j\) is defined as:

  • for \(j = 0\) \[ \frac{\partial J(\theta)}{\partial \theta_0} = \frac{1}{m}\sum\limits_{i=1}^m(h_{\theta}(x^{(i)}) - y^{(i)})x_j^{(i)} \]

  • for \(j \ge 1\)

\[ \frac{\partial J(\theta)}{\partial \theta_j} = \bigg( \frac{1}{m}\sum\limits_{i=1}^m(h_{\theta}(x^{(i)}) - y^{(i)})x_j^{(i)}\bigg) + \frac{\lambda}{m}\theta_j \]

In LinearRGradientReg.R, add code to calculate the gradient, returning it in the variable grad. When you are finished, the following code will run your gradient function using \(\theta\) initialized at [1, 1]. You should expect to see a gradient of [-15.30, 598.250].

source("LinearRGradientReg.R")
theta <- c(1, 1) 
grad <- LinearRGradientReg(theta, cbind(1, X), y, 1)


2.4. Fitting linear regression

Once your cost function and gradient are working correctly, the next code will run the function trainLinearReg.m to compute the optimal values of \(\theta\). This training function uses optim to optimize the cost function. In this part, we set the regularization parameter \(\lambda\) to zero. Because our current implementation of linear regression is trying to fit a 2-dimensional \(\theta\), regularization will not be incredibly helpful for a \(\theta\) of such low dimension. In the later parts of the exercise, you will be using polynomial regression with regularization. The code below should also plot the best fitted line, which tells us that the model does not provide a good fit as the data has a non-linear pattern.

source("trainLinearReg.R")
lambda <- 0
theta <- trainLinearReg(cbind(1, X), y, lambda)

# Plot fit over the data
plot(X, y, xlab = "Change in culture medium level", ylab = "Product flowing out of the bioreactor",
     col = "blue", pch = 19, ylim = c(-5, 40))
lines(X, cbind(1, X)%*%theta, col = "red")

While visualizing the best fit is one possible way to debug your learning algorithm, it is not always easy to visualize the data and model. In the next section, you will implement a function to generate learning curves that can help you debug your learning algorithm even when it is not easy to visualize the data.


3. Bias vs variance

An important concept in machine learning is the bias-variance tradeoff. Models with high bias are not complex enough for the data and tend to underfit, while models with high variance overfit to the training data. In this part of the exercise, you will plot training and test errors on a learning curve to diagnose bias-variance problems.


3.1. Learning curves

You will now implement code to generate the learning curves that will be useful in debugging your learning algorithms. Recall that a learning curve plots training and cross validation error as a function of training set size. Your job is to fill in learningCurve.R so that it returns a vector of errors for the training and cross validation sets.

To plot the learning curve, we need a training and cross validation set error for different training set sizes. To obtain different training set sizes, you should use different subsets of the original training set \(X\). Specifically, for a training set size of \(i\), you should use the first \(i\) examples (i.e., X[1:i, ] and y[1:i]).

You can use the trainLinearReg function to find the \(\theta\) parameters. Note that \(\lambda\) is passed as a parameter to the function. After learning the \(\theta\) parameters, you should compute the error on the training and cross validation sets. Recall that the training error for a dataset is defined as:

\[ J_{train}(\theta) = \frac{1}{2m}\Big[\sum\limits_{i=1}^m(h_\theta(x^{(i)})-y^{(i)})^2\Big] \]

In particular, note that the training error does not include the regularization term. One way to compute the training error is to use your existing cost function and set \(\lambda\) to 0 only when using it to compute the training error and cross validation error. When you are computing the training set error, make sure you compute it on the training subset (i.e., X[1:i, ] and y[1:i]) (instead of the entire training set). However, for the cross validation error, you should compute it over the entire cross validation set. You should store the computed errors in the vectors error_train and error_val. Once you are done, the following code will print the learning curves.

source("learningCurve.R")
lambda <- 0
errors <- learningCurve(cbind(1, X), y, cbind(1, Xval), yval, lambda)
error_train <- errors$error_train
error_val <- errors$error_val

plot(1:m, error_train, main = "Learning curve", 
     xlab = "Number of training examples",
     ylab = "Error", col = "blue", type = "l", 
     xlim = c(0, 13), ylim = c(0, 200))
lines(1:m, error_val, col = "darkgreen")
legend("topright", c("Training", "Cross validation"),
           col = c("blue", "darkgreen"), lty = 1)

As you can observe, both the train error and cross validation error are high when the number of training examples is increased. This reflects a high bias problem in the model: the linear regression model is too simple and thus unable to fit well our dataset. In the next section, you will implement polynomial regression to fit a better model for this dataset.


4. Polynomial regression

The problem with our linear model was that it was too simple for the data and resulted in underfitting (high bias). In this part of the exercise, you will address this problem by adding more features. For polynomial regression, our hypothesis has the form:

\[ h\theta(x) = \theta_0 + \theta_1 * (mediumLevel) + \theta_2 * (mediumLevel)^2 + ... + \theta_p * (mediumLevel)^p \\ = \theta_0 + \theta_1x_1 + \theta_2x_2 + ... + \theta_px_p \]

Notice that by defining \(x_1\) = (culture medium level), \(x_2\) = (culture medium level)\(^2\), … , \(x_p\) = (culture medium level)\(^p\), we obtain a linear regression model where the features are the various powers of the original value (culture medium level). Now, you will add more features using the higher powers of the existing feature \(x\) in the dataset. Your task in this part is to complete the code in polyFeatures.R so that the function maps the original training set \(X\) of size \(m \times 1\) into its higher powers. Specifically, when a training set \(X\) of size \(m \times 1\) is passed into the function, the function should return a \(m \times p\) matrix \(X_{poly}\) where column 1 holds the original values of \(X\), column 2 holds the values of \(X^2\), column 3 holds the values of \(X^3\), and so on. Note that you don’t have to account for the zero-eth power in this function. Once that you have a function that will map features to a higher dimension, you can apply it to the training set, the cross validation set and the test set (which you haven’t used yet).

source("polyFeatures.R")
source("featureScaling.R") # Reuse it from Practical I
p <- 8

# Map X onto polynomial features
X_poly <- polyFeatures(X, p)

# Center and scale features
fS <- featureScaling(X_poly)         
X_poly <- fS$Xnorm
mu <- fS$mu
sigma <- fS$sigma

# Add column of ones
X_poly <- cbind(1, X_poly) 

# Map X_poly_test, center and scale (using mu and sigma), add column of ones
X_poly_test <- polyFeatures(Xtest, p)
X_poly_test <- t(apply(X_poly_test, 1, function(x){(x-mu)/sigma}))
X_poly_test <- cbind(1, X_poly_test)

# Map X_poly_val, center and scale (using mu and sigma), add column of ones
X_poly_val <- polyFeatures(Xval, p)
X_poly_val <- t(apply(X_poly_val, 1, function(x){(x-mu)/sigma}))
X_poly_val <- cbind(1, X_poly_val)


4.1. Learning polynomial regression

After you have completed polyFeatures.R, the following code will proceed to train polynomial regression using your linear regression cost function. Keep in mind that even though we have polynomial terms in our feature vector, we are still solving a linear regression optimization problem. The polynomial terms have simply turned into features that we can use for linear regression. We are using the same cost function and gradient that you wrote for the earlier part of this exercise.

Here you will be using a polynomial of degree 8. It turns out that if we run the training directly on the projected data, will not work well as the features would be badly scaled (e.g., an example with \(x\) = 40 will now have a feature \(x_8 = 40^8 = 6.5 \times 10^{12}\)). That is why we need to use feature scaling. Therefore, before learning the parameters \(\theta\) for the polynomial regression, you will first call featureScaling to center and scale the features of the training set, storing the mu, sigma parameters separately. You can reuse the function you wrote for practical I.

lambda <- 0
theta <- trainLinearReg(X_poly, y, lambda)

After learning the \(\theta\) parameters, the code below will generate two plots for polynomial regression with \(\lambda\) = 0. In the first one, you should see that the polynomial fit is able to follow the datapoints almost perfectly, obtaining a low training error. However, the fit is very complex and even drops off at the extremes. This is an indicator that the polynomial regression model is overfitting the training data and will not generalize well. To better understand the problems with the unregularized (\(\lambda\) = 0) model, you can see have a look at the learning curve (second plot): the training error is low, but the cross validation error is high. There is a gap between the training and cross validation errors, indicating a high variance problem.

par(mfrow=c(1, 2))
source("plotFit.R")
# Plot fit over the data
plot(X, y, xlab = "Change in culture medium level", 
     ylab = "Product flowing out of the bioreactor",
     main = sprintf("Polynomial Regression Fit (lambda = %0.2f)", lambda),
     col = "blue", pch = 19, ylim = c(-20, 60), xlim = c(-60, 60))
plotFit(min(X), max(X), mu, sigma, theta, p)

# Plot learning curve
errors <- learningCurve(cbind(1, X_poly), y, cbind(1, X_poly_val), yval, lambda)
error_train <- errors$error_train
error_val <- errors$error_val

plot(1:m, error_train, 
     main = sprintf("Polynomial Regression Learning Curve (lambda = %0.1f)", lambda),
     xlab = "Number of training examples",
     ylab = "Error", col = "blue", type = "l", 
     xlim = c(0, 13), ylim = c(0, 200))
lines(1:m, error_val, col = "darkgreen")
legend("topright", c("Training", "Cross validation"),
           col = c("blue", "darkgreen"), lty = 1)

As you already know, one way to combat the overfitting (high-variance) problem is to add regularization to the model. In the next section, you will get to try different \(\lambda\) parameters to see how regularization can lead to a better model.


4.2. Studying the effect of the regularization parameter

In this section, you will get to observe how the regularization parameter affects the bias-variance of regularized polynomial regression. You should now modify the the lambda parameter in and try \(\lambda\) = 1, 100. For each of these values, the script should generate a polynomial fit to the data and a learning curve.

For \(\lambda\) = 1, you should see a polynomial fit that follows the data trend well and a learning curve showing that both the cross validation and training error converge to a relatively low value. This shows the \(\lambda\) = 1 regularized polynomial regression model does not have high- bias or high-variance problems. In effect, it achieves a good trade-off between bias and variance.

For \(\lambda\) = 100, you should see a polynomial fit that does not follow the data well. In this case, there is too much regularization and the model is unable to fit the training data.

par(mfrow=c(2, 2))

for (lambda in c(1, 100)){
  
  # Train the model
  theta <- trainLinearReg(X_poly, y, lambda)

  # Plot fit over the data
  source("plotFit.R")
  plot(X, y, xlab = "Change in culture medium level", 
       ylab = "Product flowing out of the bioreactor",
       main = sprintf("Polynomial Regression Fit (lambda = %0.1f)", lambda),
       col = "blue", pch = 19, ylim = c(-20, 60), xlim = c(-60, 60))
  plotFit(min(X), max(X), mu, sigma, theta, p)
  
  # Plot learning curve
  errors <- learningCurve(cbind(1, X_poly), y, cbind(1, X_poly_val), yval, lambda)
  error_train <- errors$error_train
  error_val <- errors$error_val
  
  plot(1:m, error_train, 
       main = sprintf("Polynomial Regression Learning Curve (lambda = %0.1f)", lambda),
       xlab = "Number of training examples",
       ylab = "Error", col = "blue", type = "l", 
       xlim = c(0, 13), ylim = c(0, 200))
  lines(1:m, error_val, col = "darkgreen")
  legend("topright", c("Training", "Cross validation"),
             col = c("blue", "darkgreen"), lty = 1) 
}


4.3. Selecting \(\lambda\) using a cross-validation set

From the previous parts of the exercise, you observed that the value of \(\lambda\) can significantly affect the results of regularized polynomial regression on the training and cross validation set. In particular, a model without regular- ization (\(\lambda\) = 0) fits the training set well, but does not generalize. Conversely, a model with too much regularization (\(\lambda\) = 100) does not fit the training set and testing set well. A good choice of \(\lambda\) (e.g., \(\lambda\) = 1) can provide a good fit to the data.

In this section, you will implement an automated method to select the \(\lambda\) parameter. Concretely, you will use the cross validation set to evaluate how good each \(\lambda\) value is. After selecting the best \(\lambda\) value using the cross validation set, we can then evaluate the model on the test set to estimate how well the model will perform on actual unseen data. Your task is to complete the code in validationCurve.R Specifically, you should use the trainLinearReg function to train the model using different values of \(\lambda\) and compute the training error and cross validation error. You should try \(\lambda\) in the range 0-10.

After you have completed the code, the next code will run your function to plot a cross validation curve of error vs \(\lambda\) that allows you select which \(\lambda\) parameter to use. In this figure, you will see that the best value of \(\lambda\) is around 3. Due to randomness in the training and validation splits of the dataset, the cross validation error can sometimes be lower than the training error.

source("validationCurve.R")
val_out <- validationCurve(X_poly, y, X_poly_val, yval)
lambda_vec <- val_out$lambda_vec
error_train <- val_out$error_train
error_val <- val_out$error_val

plot(lambda_vec, error_train, 
     xlab = "Number of training examples",
     ylab = "Error", col = "blue", type = "l", 
     ylim = c(0, 20))
lines(lambda_vec, error_val, col = "darkgreen")
legend("topright", c("Training", "Cross validation"),
           col = c("blue", "darkgreen"), lty = 1)


4.4. Compute the test set error

In the previous part of the exercise, you implemented code to compute the cross validation error for various values of the regularization parameter \(\lambda\). However, to get a better indication of the model’s performance in the real world, it is important to evaluate the “final” model on a test set that was not used in any part of training (that is, it was neither used to select the \(\lambda\) parameters, nor to learn the model parameters \(\theta\)). Here, you should compute the test error using the best value of \(\lambda\) you found. In our cross validation, we obtained an error of 3.8599 for \(\lambda\) = 3.

# Compute the the set error