信赖域法+狗腿法在回归分析中的应用

      第一节上节已经讲了信赖域法+狗腿法解决强凹凸二次函数，现在把这个方法引用于 Multivariance Linear Regression 中，以解决一些机器学习及其模式识别的问题。同时现在回归问题也大量运用于 Deep learning 之中。通常模式问题进行参数调整时都会涉及到学习率α的问题，而我们通常是自己给定的，比 如 0.01 或者 0.1 ，这样的结果可能会使迭代次数增加或者损失函数成震荡状况。下面我将介绍信赖域法+狗腿法在Multivariance Linear Regression中应用，以更快的速度使损失函数达到收敛。

    下面我们将介绍Multivariance Linear Regression回归的损失函数以及信赖域法+狗腿法的参数调节方式，并与梯度下降法方法进行比较。

   本文要解决的问题是给出了47个训练样本，训练样本的y值为房子的价格，x属性有2个，一个是房子的大小，另一个是房子卧室的个数。需要通过这些训练数据来学习系统的函数，从而预测房子大小为1650，且卧室有3个的房子的价格（需要的数据可以去我的资源找）。

                           

    此节的中的算法框架已经在上一节中给出，请参看。

结果图1

其中1,2，...,7 分别代表Trust Region with DogLeg，学习率为0.01, 0.03, 0.1, 0.3, 1, 1.3的最终损失函数的值。

                               
                                                

x = load('ex3x.dat');
y = load('ex3y.dat');

trustRegionBound = 1000;
x = [ones(size(x,1),1) x];
meanx = mean(x);%求均值
sigmax = std(x);%求标准偏差
x(:,2) = (x(:,2)-meanx(2))./sigmax(2);
x(:,3) = (x(:,3)-meanx(3))./sigmax(3);
itera_num = 1000; %尝试的迭代次数
sample_num = size(x,1); %训练样本的次数
figure

theta_grad_descent = zeros(size(x(1,:)));
theta = zeros(size(x,2),1); %theta的初始值赋值为0
Jtheta = zeros(itera_num, 1);
 for i = 1:itera_num %计算出某个学习速率alpha下迭代itera_num次数后的参数       
        Jtheta(i) = (1/(2*sample_num)).*(x*theta-y)'*(x*theta-y);%Jtheta是个行向量
        grad = (1/sample_num).*x'*(x*theta-y);
        B=x'*x;
        du = -grad' * grad * grad / (grad' * B * grad);
        dB = -B^-1 * grad;
        a = 2;
        if du'*du > trustRegionBound*trustRegionBound;
        a = trustRegionBound / sqrt((du'*du));
        else if dB'*dB > trustRegionBound*trustRegionBound
        a = sqrt((trustRegionBound*trustRegionBound - du'*du) / ((dB-du)'*(dB-du))) + 1;
            end    
        end
        if a < 1
        d = a * du;
        else
        d = du + (a - 1) * (dB - du);
        end
        Jtheta1(i)=(1/(2*sample_num)).*(x*(theta+d)-y)'*(x*(theta+d)-y);
        p = (Jtheta(i)-Jtheta1(i))/(-grad'*d-1/2*d'*B*d);
            if p > 0.75 && sqrt(abs(d'*d) - trustRegionBound) < 0.001
        trustRegionBound = min(2 * trustRegionBound, 10000);
           else if p < 0.25
            trustRegionBound = sqrt(abs(d'*d)) * 0.25;
               end
            end
         if p > 0%q(zeros(2,1),x) > q(d, x)
        theta = theta + d;
         end
 end
K(1)=Jtheta(1000)
    plot(0:350, Jtheta(1:351),'k--','LineWidth', 2)%此处一定要通过char函数来转换
    hold on
alpha = [0.01, 0.03, 0.1, 0.3, 1, 1.3];%因为差不多是选取每个3倍的学习率来测试，所以直接枚举出来
plotstyle = {'b', 'r', 'g', 'k', 'b--', 'r--'};

theta_grad_descent = zeros(size(x(1,:)));
for alpha_i = 1:length(alpha) %尝试看哪个学习速率最好
    theta = zeros(size(x,2),1); %theta的初始值赋值为0
    Jtheta = zeros(itera_num, 1);
    for i = 1:itera_num %计算出某个学习速率alpha下迭代itera_num次数后的参数       
        Jtheta(i) = (1/(2*sample_num)).*(x*theta-y)'*(x*theta-y);%Jtheta是个行向量
        grad = (1/sample_num).*x'*(x*theta-y);
        theta = theta - alpha(alpha_i).*grad;
    end
    K(alpha_i+1)=Jtheta(1000);
    plot(0:350, Jtheta(1:351),char(plotstyle(alpha_i)),'LineWidth', 2)%此处一定要通过char函数来转换
    hold on
end
legend('Trust Region with DogLeg','0.01','0.03','0.1','0.3','1','1.3');
xlabel('Number of iterations')
ylabel('Cost function')
figure
plot(1:7,K,'b-','LineWidth', 2);