机器学习练习之正则化

 这个练习是来自http://openclassroom.stanford.edu/MainFolder/DocumentPage.php?course=MachineLearning&doc=exercises/ex5/ex5.html

Regularized linear regression

        之前已经做过线性回归， 我们知道线性回归就是要找到 一条直线去拟合训练数据，模型的损失函数为

        

        但是，这个模型可能会有些权重很大，有些权重很小，导致过拟合，就是过分拟合了训练数据，使得模型的复杂度提高，泛化能力较差，泛化能力就是对未知数据的预测能力。为了防止过拟合，通常会加入权重惩罚项，就是模型的正则项。正则项可以取不同的形式，在回归问题中取平方损失，就是参数的L2范数，也可以取L1范数。取平方损失时，模型的损失函数变为：

      

      lambda是正则项系数，如果它的值很大，说明对模型的复杂度惩罚大，对拟合数据的损失惩罚小，这样它就不会过分拟合数据，在训练数据上的偏差较大，在未知数据上的variance较小，但是可能出现欠拟合的现象；相反，如果lambda值很小，说明比较注重对训练数据的拟合，在训练数据上的bias会小，但是可能会导致过拟合。

      

     要求出theta使得损失J最小，可以将J对theta求导，得到结果如下：

    

    所以直接使用这个方程解就可以求出最佳的theta，不过也可以使得梯度下降算法，但是要折腾一下调步长参数。

  

    程序如下：

[plain]  view plain copy

1. %% Exercie: Regularized linear regression  
2. % 正则化的线性回归  
3. % 在简单的线性回归上加了正则项，用于惩罚权重，减少模型参数的搜索空间，  
4. % 使得权重不会偏向于某些特征，模型更加平滑  
5. %%  
6. %% 数据初始化  
7. x = load('ex5Linx.dat'); % x是只有一个特征，为一个列向量  
8. y = load('ex5Liny.dat');  
9. % 画图  
10. figure  
11. plot(x, y, 'o', 'MarkerFacecolor', 'r', 'MarkerSize', 8);  
12. m = length(x); % m样本数量  
13. x = [ones(m,1), x, x.^2, x.^3, x.^4, x.^5];  % 根据x构造特征向量，这里是多项式的特征，每个特征取(0-5)次方  
14. n = size(x, 2); % 特征数量  
15. %%  
16. %% 梯度下降算法  
17. % theta = zeros(n, 1); %初始化权重向量,一个列向量  
18. % alpha = 0.01; %学习率，步长  
19. % lambda = 0;  
20. % MAX_ITR = 500;  
21. % for i=1:MAX_ITR  
22. %     grad = (1/m).*(x'*(x*theta-y) + lambda*theta); %计算梯度  
23. %     theta = theta - alpha*grad; %更新权重  
24. % end  
25. % theta  
26. % J = (0.5/m) * (sum((x*theta-y).^2) + lambda*sum(theta.^2))   
27. % theta_norm = norm(theta)  
28. %%  
29. %% Normal equations  
30. lambda = [0 1 10];  
31. plotstyle = {'b', 'r', 'g'}  
32. for i = 1:3  
33.     L = lambda(i).*eye(n);  
34.     L(1,1) = 0;  
35.     theta = (x'*x+L)\x'*y;  
36.     theta  
37.     J = (0.5/m) * (sum((x*theta-y).^2) + lambda(i)*sum(theta.^2))   
38.     theta_norm = norm(theta)    %计算theta的L2范数  
39.     hold on  
40.     x_vals = (-1:0.05:1)'; %这个范围要根据x的范围确定  
41.     features = [ones(size(x_vals)), x_vals, x_vals.^2, x_vals.^3, ...  
42.             x_vals.^4, x_vals.^5];  
43.     plot(x_vals, features*theta, char(plotstyle(i)), 'LineWidth', 2)  
44.     legend(['lambda=', num2str(i)])  
45. end  
46. legend('Training data','\lambda=0', '\lambda=1', '\lambda=10');  
47. title('5th order fit')  
48. hold off  
49. %%  
50.       
51.       
52.       

结果：

Regularized Logistic regression

      正则化的逻辑斯蒂回归模型也是在之前的 logistic regression上增加了正则项。

   

     要使得上面的损失函数最小，同样这样使用牛顿方法，它的更新规则为：

   

      其中一阶导和海森矩阵的形式有所变化：

     

   

   程序如下：

[plain]  view plain copy

1. function out = map_feature(feat1, feat2)  
2. % MAP_FEATURE    Feature mapping function for Exercise 5  
3. %  
4. %   map_feature(feat1, feat2) maps the two input features  
5. %   to higher-order features as defined in Exercise 5.  
6. %  
7. %   Returns a new feature array with more features  
8. %  
9. %   Inputs feat1, feat2 must be the same size  
10. %  
11. % Note: this function is only valid for Ex 5, since the degree is  
12. % hard-coded in.  
13.     degree = 6;  
14.     out = ones(size(feat1(:,1)));  
15.     for i = 1:degree  
16.         for j = 0:i  
17.             out(:, end+1) = (feat1.^(i-j)).*(feat2.^j);  
18.         end  
19.     end  
20.       

[plain]  view plain copy

1. %% Regularized logistic regression  
2. % 加入正则项的逻辑斯蒂回归模型  
3. %%  
4. %% 数据初始化  
5. x = load('ex5Logx.dat');  % 这里x的每一行是一个二维向量，表示有两个特征  
6. y = load('ex5Logy.dat');  
7. % Find the indices for the 2 classes  
8. pos = find(y); neg = find(y==0);  
9. x_unfeatured = x; % 保存未构造特征前的x  
10. x = map_feature(x(:,1), x(:,2)); %由于原始数据下样本不是线性可分的，所以这里根据x的2个特征构造更多的特征  
11. [m n] = size(x);   
12. %%  
13. %% 定义sigmoid函数  
14. g = inline('1.0 ./ (1.0+exp(-z))', 'z');  
15. %%   
16. %% 牛顿算法  
17. lambda = [0 1 10];  %正则项系数  
18. MAX_ITR = 15;  
19. for i = 1:3  
20.     theta = zeros(n, 1); %权重初始化  
21.     J = zeros(MAX_ITR, 1); %保存损失  
22.     for j = 1:MAX_ITR  
23.         h = g(x*theta); %计算sigmoid值  
24.         L = lambda(i) .* eye(n);  
25.         L(1,1) = 0;  
26.         J(j) = -(1.0/m).*sum(y.*log(h) + (1-y).*log(1-h)) + ...  
27.             lambda(i)/(2*m).*norm(theta(2:end))^2;  
28.         grad = (1/m) .* (x'*(h-y) + L*theta); %计算梯度  
29.         H = (1/m) .* x'*diag(h)*diag(1-h)*x + lambda(i)/m.*L;  %计算海森矩阵  
30.           
31.         theta = theta - H\grad;  
32.     end  
33.     J  
34.     norm_theta = norm(theta)  
35.     figure  
36.     plot(x_unfeatured(pos, 1), x_unfeatured(pos, 2), 'k+', 'LineWidth', 2, 'MarkerSize', 7)  
37.     hold on  
38.     plot(x_unfeatured(neg, 1), x_unfeatured(neg, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7)  
39.     hold on  
40.     %显示结果  
41.     u = linspace(-1, 1.5, 200);  
42.     v = linspace(-1, 1.5, 200);  
43.     z = zeros(length(u), length(v));  
44.     for ii = 1:length(u)  
45.         for jj = 1:length(v)  
46.             z(ii, jj) = map_feature(u(ii), v(jj))*theta;  
47.         end  
48.     end  
49.     z = z';  
50.     % plot z = 0  
51.     contour(u, v, z, [0,0], 'LineWidth', 2)  
52.     legend('y=1', 'y=0', 'Decision boundary')  
53.     title(sprintf('\\lambda = %g', lambda(i)), 'FontSize', 14)  
54.     hold off  
55. end  
56.       

结果：