# Something wrong with my implementation of the bias/variance diagnostic in polynomial regression

I'm trying to diagnosing bias/variance so I have the below Octavecode:

%========== Begin - constants declaration ==========%
x_training_percent = 0.7;
cv_set_percent = 0.3;
%========== End - constants declaration ==========%

data = [X y];
data = data(:, randperm(size(data,2)));

[m, n] = size(X);% m: Number of examples, n: Number of features.

%========== Begin - Getting traingin and CV sets ==========%
training_set_size = round(m * x_training_percent);
cv_set_size = round(m * cv_set_percent);
%test_set_size = round(m * test_set_percent);

x_training_o = data(1:training_set_size, 1:(end - 1));
y_training = data(1:training_set_size, end);

x_cv_o = data((training_set_size + 1):(training_set_size + cv_set_size), 1:(end - 1));
y_cv = data((training_set_size + 1):(training_set_size + cv_set_size), end);

%========== End - Getting traingin and CV sets ==========%

max_p = 20; % Max degree polynomial

cv_error = zeros(max_p, 1);
training_error = zeros(max_p, 1);

for p = 1:max_p

% Processing training set
x_training = x_training_o;
x_training = polyFeatures(x_training, p); % Adding polynomial terms from 1 to p
x_training = [ones(training_set_size, 1) x_training];

% Processing cross validation set
x_cv = x_cv_o;
x_cv = polyFeatures(x_cv, p); % Adding polynomial terms from 1 to p
x_cv = [ones(cv_set_size, 1) x_cv];

%========== Begin - Training ==========%
lambda = 0
theta = trainLinearReg(x_training, y_training, lambda);
%========== End - Training ==========%

%========== Begin - Computing prediction errors with polinomial degree p ==========%
predictions = x_training * theta; % Predictions with training set
training_error(p, :) = (1 / (2 * training_set_size)) * sum((predictions - y_training) .^ 2);

cv_predictions = x_cv * theta; % Predictions with cross validation set
cv_error(p, :) = (1 / (2 * cv_set_size)) * sum((cv_predictions - y_cv) .^ 2);
%========== End - Computing prediction errors ==========%

end
plot(1:max_p, training_error, 1:max_p, cv_error);

legend('Train', 'Cross Validation')
xlabel('Degree of Polynomial')
ylabel('Error')


My outputs is:

But it's very different from the following result shown by Andrew Ng:

So what am I doing wrong here? P.D: You can see the complete source code here

• Given that higher-order polynomials must fit the training data better than the lower-order polynomials, your growing blue curve clearly indicates that you have a bug. It is still fine to ask this such a question, but I have edited the title accordingly (and added [matlab] tag). – amoeba says Reinstate Monica May 18 '16 at 16:19
• I have checked exhaustively but I found't bugs. 40% of features are booleans (0 and 1) so maybe adding polynomial degrees doesn't help to fit correctly (because 0^n is always 0 and 1^n is always 1). – Overflow012 May 18 '16 at 22:35
• It might not help, but it cannot hurt (on the training set). – amoeba says Reinstate Monica May 18 '16 at 22:39
• It looks okay, so maybe the bug is in one of the functions you use. Can you show the code of trainLinearReg? – tomka May 24 '16 at 21:10

Here is how you should debug your code:

1. For now, forget about cross-validation. Your most serious problem is the training error increasing with the polynomial degree, while it should decrease (and become zero for big enough polynomial degree). Make the code without any traces of cross-validation. It will be much easier to check.

2. If the error will not be reproduced (when you remove the cross-validation from your code, the training error suddenly decreases), then add the cross-validation part gradually (small chunks of code) until you find which part spoils the correct work, and understand why. Otherwise, go to step 3.

3. Now you have a code without cross-validation, which does not work properly. Make the simulated data which ideally fits, say, quadratic polynomial, and check if it works properly. If it doesn't, try the linear polynomial. At the end you will have something very simple which does not work. You will probably be able to figure out why.

4. Posting your whole code which does not work and asking people "what did I do wrong" is a bad idea. It means that you want people to debug your code instead of you. Why should they do it? You should first try to make your code as small as possible while still retaining the error. If you can remove something from your code, and the error is still here, then the code is NOT as simple as possible. If there is a part of your code which you did not try to remove, the code is NOT as simple as possible.

• "Why should they do it?" --> because they get +100? – tomka May 24 '16 at 21:12
• @ThomasKlausch it is absurd to do debugging of someone's else code for scores. – user31264 May 25 '16 at 0:12
• stackoverflow.com/help/mcve applies here too. – Nick Cox May 25 '16 at 10:55

I had a very similar problem before, and it was because the step size in my gradient descent implementation was too large and the tolerance too forgiving.

As you add higher order features their magnitude may be rapidly increasing, making deviations from the true coefficient values cause greater error. The fact that your training error increases as you add higher order features suggests this may be happening.

I have added feature normalization (mean = 0 and standard deviation = 1) and seems to be that I found the problem. See the new curve (Just training set). However, error is so high

And the code below...

addpath ./common;

data = [X y];
data = data(:, randperm(size(data,2)));

[m, n] = size(X);% m: Number of examples, n: Number of features.

%========== Begin - Getting traingin and CV sets ==========%
x_training_o = data(1:m, 1:(end - 1));
y_training = data(1:m, end);
%========== End - Getting traingin and CV sets ==========%

max_p = 10; % Max degree polynomial

training_error = zeros(max_p, 1);

csvwrite("log/Xy_tr.txt", [x_training_o y_training]);

for p = 1:max_p
x_training = x_training_o;
x_training = polyFeatures(x_training, p); % Adding polynomial terms from 1 to p
[x_training, mu, sigma] = featureNormalize(x_training);
x_training = [ones(m, 1) x_training];

%========== Begin - Training ==========%
lambda = 0
theta = trainLinearReg(x_training, y_training, lambda);
%========== End - Training ==========%

%========== Begin - Computing prediction errors with polinomial degree p ==========%
predictions = x_training * theta; % Predictions with training set

training_error(p, :) = (1 / (2 * m)) * sum((predictions - y_training) .^ 2);

%========== End - Computing prediction errors ==========%

end

• The training curve should decrease independent of feature normalization. – amoeba says Reinstate Monica May 25 '16 at 23:57
• @amoeba, it's the only change that I made. – Overflow012 May 26 '16 at 0:33
• I afraid you still have the bug, for the reason amoeba explained. Besides, why would your feature normalization increased the mean squared error from about 0.1 to about $3 \cdot 10^7$ (many orders of magnitude). Your MSE should remain the same under any linear transformations. – user31264 May 26 '16 at 11:04

I suspect you may find better results when you use the ordinary least squares estimator instead of gradient descent:

$$\hat{ \beta } = (X'X)^{-1} X'Y.$$