Julia中的线性回归和矩阵除法答案

【问题标题】：Linear regression and matrix division in JuliaJulia中的线性回归和矩阵除法
【发布时间】：2018-07-25 10:52:50
【问题描述】：

众所周知的OLS公式是(X'X)^(-1)X'y，其中X是nxK，y是nx1。

在 Julia 中实现此功能的一种方法是 (X'*X)\X'*y。

但我发现X\y 给出了几乎相同的输出，直到出现微小的计算错误。

他们是否总是计算相同的东西（只要n>k）？如果是这样，我应该使用哪一个？

【问题讨论】：

标签： julia linear-algebra

【解决方案1】：

当X 为正方形时，存在唯一解，LU-factorization (with pivoting) 是一种数值稳定的计算方法。这就是反斜杠在这种情况下使用的算法。

当X 不是平方时（大多数回归问题都是这种情况），则不存在唯一解，但存在唯一最小二乘解。求解Xβ = y 的QR factorization 方法是一种生成最小二乘解的数值稳定方法，在这种情况下X\y 使用QR 分解，从而得到OLS 解。

注意数字稳定的词。虽然(X'*X)\X'*y 理论上总是给出与反斜杠相同的结果，但实际上反斜杠（使用正确的因式分解选择）会更精确。这是因为分解算法被实现为numerically stable。由于在执行(X'*X)\X'*y 时会累积浮点错误，因此不建议您将此表单用于任何实际的数值工作。

相反，(X'*X)\X'*y 有点等价于 SVD 分解，这是数值上最稳定的算法，但也是最昂贵的（事实上，它基本上是写出 Moore-Penrose pseudoinverse，这就是 SVD factorization 用于求解线性系统）。要使用旋转 SVD 直接进行 SVD 分解，请在 v0.6 上执行 svdfact(X) \ y 或在 v0.7 上执行 svd(X) \ y。直接这样做比(X'*X)\X'*y 更稳定。请注意，qrfact(X) \ y 或 qr(X) \ y (v0.7) 用于 QR。 See the factorizations portion of the documentation 了解有关所有选择的更多详细信息。

【讨论】：

【解决方案2】：

按照文档，X\y 的结果是（使用符号 \(A, B) 而不是 X 和 y）：

对于矩形 A，结果是最小范数最小二乘解

这是你的情况，我猜你假设n>k（所以你的矩阵不是正方形的）。所以你可以放心地使用X\y。实际上，使用它比使用标准公式更好，因为即使X 的等级小于min(n,k)，您也会得到结果，而标准公式(X'*X)^(-1)*X'*y 将失败或在X'*X 接近时产生数值不稳定的结果单数。

如果X 是方形的（这不是你的情况），那么我们在文档中有一些不同的规则：

对于输入矩阵 A 和 B，当 A 为正方形时，结果 X 是 A*X == B

这意味着如果您的矩阵是奇异的，\ 算法会产生错误，或者如果矩阵几乎是奇异的，则会产生数值不稳定的结果（实际上最常见的 lu 函数在内部为一般密集矩阵调用可能会抛出SingularException)。

如果您想要一个包罗万象的解决方案（对于方阵和非方阵），那么可以使用 qr(X, Val(true)) \ y。

【讨论】：

一个更好的包罗万象的是 SVD，但它比 QR 贵。有些矩阵可以在 SVD 下稳定求解，但在 QR 下则不稳定。 This has an extended discussion on the numerical stability of the two algorithms.
同意，我使用 QR，因为这是 Julia 默认使用的 AFAIK。

【解决方案3】：

简答：不，使用第一个（众所周知的）。

长答案：

线性回归模型是Xβ = y，很容易推导出β = X \ y，这是你的第二种方法。但是，在大多数情况下（当X 不可逆时），这是错误的，因为您不能简单地左乘X^-1。正确的做法是改为解决β = argmin{‖y - Xβ‖^2}，这就引出了第一种方法。

为了表明它们并不总是相同的，简单地构造一个 X 不可逆的情况：

julia> X = rand(10, 10)
10×10 Array{Float64,2}:
 0.938995  0.32773   0.740556  0.300323   0.98479    0.48808    0.748006   0.798089  0.864154  0.869864
 0.973832  0.99791   0.271083  0.841392   0.743448   0.0951434  0.0144092  0.785267  0.690008  0.494994
 0.356408  0.312696  0.543927  0.951817   0.720187   0.434455   0.684884   0.72397   0.855516  0.120853
 0.849494  0.989129  0.165215  0.76009    0.0206378  0.259737   0.967129   0.733793  0.798215  0.252723
 0.364955  0.466796  0.227699  0.662857   0.259522   0.288773   0.691278   0.421251  0.593215  0.542583
 0.126439  0.574307  0.577152  0.664301   0.60941    0.742335   0.459951   0.516649  0.732796  0.990509
 0.430213  0.763126  0.737171  0.433884   0.85549    0.163837   0.997908   0.586575  0.257428  0.33239
 0.28398   0.162054  0.481452  0.903363   0.780502   0.994575   0.131594   0.191499  0.702596  0.0967979
 0.42463   0.142     0.705176  0.0481886  0.728082   0.709598   0.630134   0.139151  0.423227  0.942262
 0.197805  0.526095  0.562136  0.648896   0.805806   0.168869   0.200355   0.557305  0.69514   0.227137

julia> y = rand(10, 1)
10×1 Array{Float64,2}:
 0.7751785556478308
 0.24185992335144801
 0.5681904264574333
 0.9134364924569847
 0.20167825754443536
 0.5776727022413637
 0.05289808385359085
 0.5841180308242171
 0.2862768657856478
 0.45152080383822746

julia> ((X' * X) ^ -1) * X' * y
10×1 Array{Float64,2}:
 -0.3768345891121706
  0.5900885565174501
 -0.6326640292669291
 -1.3922334538787071
  0.06182039005215956
  1.0342060710792016
  0.045791973670925995
  0.7237081408801955
  1.4256831037950832
 -0.6750765481219443

julia> X \ y
10×1 Array{Float64,2}:
 -0.37683458911228906
  0.5900885565176254
 -0.6326640292676649
 -1.3922334538790346
  0.061820390052523294
  1.0342060710793235
  0.0457919736711274
  0.7237081408802206
  1.4256831037952566
 -0.6750765481220102

julia> X[2, :] = X[1, :]
10-element Array{Float64,1}:
 0.9389947787349187
 0.3277301697101178
 0.7405555185711721
 0.30032257202572477
 0.9847899425069042
 0.48807977638742295
 0.7480061513093117
 0.79808859136911
 0.8641540973071822
 0.8698636291189576

julia> ((X' * X) ^ -1) * X' * y
10×1 Array{Float64,2}:
  0.7456524759867015
  0.06233042922132548
  2.5600126098899256
  0.3182206475232786
 -2.003080524452619
  0.272673133766017
 -0.8550165639656011
  0.40827327221785403
  0.2994419115664999
 -0.37876151249955264

julia> X \ y
10×1 Array{Float64,2}:
  3.852193379477664e15
 -2.097948470376586e15
  9.077766998701864e15
  5.112094484728637e15
 -5.798433818338726e15
 -2.0446050874148052e15
 -3.300267174800096e15
  2.990882423309131e14
 -4.214829360472345e15
  1.60123572911982e15

【讨论】：

正如我在回答 X\y 中所指出的，如果 X 不是正方形，则实际上给出了 argmin{‖y - Xβ‖^2} 的解决方案，并且比标准公式更有效。如果X 是正方形（但这不是问题，因为假设n>k）并且排名不足，如您的示例所示，实际上理论上这两种方法都应该引发错误。您的示例仅适用于数字舍入错误。原则上(X'*X)^(-1) 应该为秩低矩阵X 抛出错误。
这个答案不正确。正如@BogumiłKamiński 所说，当矩阵不是方阵时，使用 QR，它给出最小二乘解，因此X\y 是 OLS 解。并且通过 QR 的反斜杠在数值上比“传统公式”更稳定，因此应该始终使用反斜杠。