sho*_*oli 5 linear-algebra julia
对于OLS众所周知的公式是(X'X)^(-1)X'y其中X是nxK和y是nx1.
在Julia中实现这一点的一种方法是(X'*X)\X'*y.
但我发现,X\y几乎相同的输出可以达到微小的计算误差.
他们总是计算相同的东西(只要n>k)？如果是这样,我应该使用哪一个？
什么X是正方形时,有一个独特的解决方案,LU分解(带有旋转)是一种数值稳定的计算方法.这是反斜杠在这种情况下使用的算法.
如果X不是正方形,这是大多数回归问题的情况,那么没有唯一的解决方案,但有一个独特的最小二乘解决方案.用于求解的QR分解方法X? = y是用于生成最小二乘解的数值稳定方法,并且在这种情况下X\y使用QR分解并因此给出OLS解.
注意这些词在数值上是稳定的.虽然(X'*X)\X'*y理论上总是会给出与反斜杠相同的结果,但实际上反斜杠(使用正确的分解选择)会更精确.这是因为分解算法被实现为在数值上稳定.由于浮点误差的变化在进行时会累积(X'*X)\X'*y,因此不建议您将此表单用于任何实际的数值工作.
相反,(X'*X)\X'*y它在某种程度上等同于SVD因子分解,这是最经济的稳定算法,但也是最昂贵的(事实上,它基本上写出了Moore-Penrose伪逆,这是SVD分解用于求解线性系统的方式).要使用旋转SVD直接进行SVD分解,请执行v0.6 svdfact(X) \ y或svd(X) \ yv0.7.直接这样做比比这更稳定(X'*X)\X'*y.请注意,qrfact(X) \ y或qr(X) \ y(v0.7)适用于QR.有关所有选项的更多详细信息,请参阅文档的分解部分.
继文档的结果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可以使用.