用于求解线性回归的梯度下降和正规方程方法给出了不同的解

Question

我正在研究机器学习问题,并希望使用线性回归作为学习算法.我已经实现了两种不同的方法来查找theta线性回归模型的参数:梯度(最陡)下降和法线方程.在相同的数据上,它们都应该给出近似相等的theta向量.但他们没有.

两个theta向量在所有元素上非常相似,但第一个.这是用于将添加到数据中的所有1的向量相乘的那个.

这是thetas的样子(第一列是Gradient下降的输出,是Normal方程的第二输出):

Grad desc Norm eq
-237.7752 -4.6736
-5.8471   -5.8467
9.9174    9.9178
2.1135    2.1134
-1.5001   -1.5003
-37.8558  -37.8505
-1.1024   -1.1116
-19.2969  -19.2956
66.6423   66.6447
297.3666  296.7604
-741.9281 -744.1541
296.4649  296.3494
146.0304  144.4158
-2.9978   -2.9976
-0.8190   -0.8189

有什么可以导致差异theta(1, 1)的梯度下降返回相比,theta(1, 1)通过正常的公式返回？我的代码中有错误吗？

这是我在Matlab中实现的正规方程:

function theta = normalEque(X, y)
    [m, n] = size(X);
    X = [ones(m, 1), X];
    theta = pinv(X'*X)*X'*y;
end

这是梯度下降的代码:

function theta = gradientDesc(X, y)
    options = optimset('GradObj', 'on', 'MaxIter',  9999);
    [theta, ~, ~] = fminunc(@(t)(cost(t, X, y)),...
                    zeros(size(X, 2), 1), options);
end

function [J, grad] = cost(theta, X, y)
    m = size(X, 1);
    X = [ones(m, 1), X];
    J = sum((X * theta - y) .^ 2) ./ (2*m);
    for i = 1:size(theta, 1)
        grad(i, 1) = sum((X * theta - y) .* X(:, i)) ./ m;
    end
end

我传递完全相同的数据X和y两个函数(我没有规范化X).

编辑1:

基于答案和评论,我检查了几个代码并运行了一些测试.

首先,我想根据@ user1489497的回答,检查问题是否可能是由于单独的X beeing引起的.所以我用inv替换了pinv - 当它运行时我真的得到了警告Matrix is close to singular or badly scaled..为了确保这不是问题,我获得了更大的数据集并使用这个新数据集运行测试.这次inv(X)没有显示警告并使用pinv并inv给出相同的结果.所以我希望X不再接近单数.

然后,我改变了normalEque代码,建议由木片所以现在它看起来像:

function theta = normalEque(X, y)
    X = [ones(size(X, 1), 1), X];
    theta = pinv(X)*y;
end

但问题仍然存在.normalEque不接近单数的新数据的新函数给出不同theta的gradientDesc.

为了找出哪个算法是错误的,我在相同的数据上运行了数据挖掘软件Weka的线性回归算法.Weka计算theta非常类似于输出normalEque但与输出不同gradientDesc.所以我想这normalEque是正确的,并且有一个错误gradientDesc.

这是thetaWeka计算的s的比较,normalEque并且GradientDesc:

Weka(correct) normalEque    gradientDesc
779.8229      779.8163      302.7994
  1.6571        1.6571        1.7064
  1.8430        1.8431        2.3809
 -1.5945       -1.5945       -1.5964
  3.8190        3.8195        5.7486
 -4.8265       -4.8284      -11.1071
 -6.9000       -6.9006      -11.8924
-15.6956      -15.6958      -13.5411
 43.5561       43.5571       31.5036
-44.5380      -44.5386      -26.5137
  0.9935        0.9926        1.2153
 -3.1556       -3.1576       -1.8517
 -0.1927       -0.1919       -0.6583
  2.9207        2.9227        1.5632
  1.1713        1.1710        1.1622
  0.1091        0.1093        0.0084
  1.5768        1.5762        1.6318
 -1.3968       -1.3958       -2.1131
  0.6966        0.6963        0.5630
  0.1990        0.1990       -0.2521
  0.4624        0.4624        0.2921
-12.6013      -12.6014      -12.2014
 -0.1328       -0.1328       -0.1359

我还根据Justin Peel的回答计算出错误.输出normalEque给出略小的平方误差,但差异很小.更重要的是,当我计算theta使用函数的成本梯度(与使用的函数cost相同gradientDesc)时,我得到了接近零的梯度.在输出上完成同样gradientDesc不会给出接近零的梯度.这就是我的意思:

>> [J_gd, grad_gd] = cost(theta_gd, X, y, size(X, 1));
>> [J_ne, grad_ne] = cost(theta_ne, X, y, size(X, 1));
>> disp([J_gd, J_ne])
  120.9932  119.1469
>> disp([grad_gd, grad_ne])
  -0.005172856743846  -0.000000000908598
  -0.026126463200876  -0.000000135414602
  -0.008365136595272  -0.000000140327001
  -0.094516503056041  -0.000000169627717
  -0.028805977931093  -0.000000045136985
  -0.004761477661464  -0.000000005065103
  -0.007389474786628  -0.000000005010731
   0.065544198835505  -0.000000046847073
   0.044205371015018  -0.000000046169012
   0.089237705611538  -0.000000046081288
  -0.042549228192766  -0.000000051458654
   0.016339232547159  -0.000000037654965
  -0.043200042729041  -0.000000051748545
   0.013669010209370  -0.000000037399261
  -0.036586854750176  -0.000000027931617
  -0.004761447097231  -0.000000027168798
   0.017311225027280  -0.000000039099380
   0.005650124339593  -0.000000037005759
   0.016225097484138  -0.000000039060168
  -0.009176443862037  -0.000000012831350
   0.055653840638386  -0.000000020855391
  -0.002834810081935  -0.000000006540702
   0.002794661393905  -0.000000032878097

这表明梯度下降根本没有收敛到全局最小值......但是,当我运行数千次迭代时,情况就不是这样了.那么bug在哪里？

Answer 1

我终于有时间回到这里了.没有"错误".

如果矩阵是单数的,那么存在无限多的解.您可以从该组中选择任何解决方案,并获得同样好的答案.pinv(X)*y解是一个很好的解决方案,因为它是最小的范数解决方案.

绝对没有充分的理由使用inv(X)*y.更糟糕的是,对正规方程使用逆,因此inv(X'*X)*X'*y只是数值废话.我不在乎是谁告诉你使用它,它们会引导你到错误的地方.(是的,对于条件良好的问题,它可以接受,但大多数时候你不知道什么时候会给你废话.所以为什么要使用它？)

正常方程通常是一件坏事,即使您正在解决正则化问题.有办法做到这一点,避免平方系统的条件数,虽然我不会解释它们,除非被问及这个答案已经足够长.

X\y也会产生合理的结果.

绝对没有理由在这个问题上抛出一个无约束的优化器,因为这会产生不稳定的结果,完全取决于你的起始值.

举个例子,我将从一个单一的问题开始.

X = repmat([1 2],5,1);
y = rand(5,1);

>> X\y
Warning: Rank deficient, rank = 1, tol =  2.220446e-15. 
ans =
                         0
         0.258777984694222

>> pinv(X)*y
ans =
         0.103511193877689
         0.207022387755377

pinv和反斜杠返回略有不同的解决方案.事实证明,有一个基本的解决方案,我们可以为X的行空间添加任意数量的零空间向量.

null(X)
ans =
         0.894427190999916
        -0.447213595499958

pinv生成最小范数解.在可能产生的所有解决方案中,这个解决方案至少有2个标准.

相反,反斜杠生成的解决方案将一个或多个变量设置为零.

但是如果使用无约束优化器,它将生成一个完全取决于起始值的解决方案.同样,可以将任意数量的空向量添加到您的解决方案中,并且您仍然拥有完全有效的解决方案,并且具有相同的误差平方值.

请注意,即使没有返回奇点,也不一定意味着你的矩阵不接近单数.你对这个问题几乎没有什么改变,所以它仍然很接近,只是不足以触发警告.