Matlab最小二乘估计通过\用于R中的欠定系统答案

【问题标题】：Matlabs least square estimate via \ for underdetermined system in RMatlab最小二乘估计通过\用于R中的欠定系统
【发布时间】：2019-12-16 13:22:42
【问题描述】：

我想找到一个欠定线性方程组的最小二乘估计，就像 Matlab 对 '\' 所做的那样。

我试图在 R 中重现的内容

% matlab code

X = [3.8642    9.6604;
   14.2000   35.5000;
   41.7832  104.4580;
    0.4084    1.0210];

y = [1.2300
    4.5200
   13.3000
    0.1300];

X\y %  => [0, 0.1273]

我尝试了 R 的 lsfit 方法，MASS 包中的广义逆 (ginv)，并使用 QR 组合 (R^{-1}Q'y)，但都返回不同的结果。

【问题讨论】：

很好地提供了您想要模拟的 MATLAB 代码，但就 minimal reproducible example 而言，最好还提供数据以供人们轻松剪切并粘贴到 R

标签： r matlab linear-equation

【解决方案1】：

为了将来参考，正如我在评论中所说，以预期回答者可以轻松使用的格式提供示例输入被认为是一种很好的做法（它可以帮助其他人帮助您！）。在这里，您需要一个 R 解决方案，因此示例输入也应该以在 R 中易于使用的方式提供。您可以提供以下方式：

x <- matrix(c(3.8642, 14.2, 41.7832, 0.4084, 9.6604, 35.5, 104.458, 1.021),
            ncol = 2)
y <- c(1.23, 4.52, 13.3, 0.13)

解决了这部分问题，我们可以继续这个问题：如何做相当于 MATLAB 的\ 的最小二乘估计？ The documentation for \ 说

x = A\B 求解线性方程组 A*x = B。

您在 R 中寻找的等价物是 solve() 或 qr.solve()，在这种情况下：

?solve

这个通用函数求解方程 a %*% x = b for x
. . .
qr.solve 可以处理非正方形系统。

所以，我们可以看到

qr.solve(x, y)
# [1] 0.003661243 0.125859408

这非常接近您的 MATLAB 解决方案。同样，lsfit() 或 lm()（当然）给你同样的答案：

coef(lm(y ~ x + 0))
#          x1          x2 
# 0.003661243 0.125859408 
lsfit(x, y, intercept = FALSE)$coef
#          X1          X2 
# 0.003661243 0.125859408

我们可以看到这个答案至少与您的 MATLAB 解决方案一样适合数据：

r_solution <- coef(lm(y ~ x + 0))
y - (x %*% r_solution)
              [,1]
[1,] -1.110223e-15
[2,]  1.366296e-06
[3,] -4.867456e-07
[4,]  2.292817e-06
matlab_solution <- c(0, 0.1273)
y - (x %*% matlab_solution)
           [,1]
[1,] 0.00023108
[2,] 0.00085000
[3,] 0.00249660
[4,] 0.00002670

【讨论】：

【解决方案2】：

当我使用ginv() 时，我能够达到与 MATLAB/OCTAVE 相同的结果：

在R中使用ginv()

library(MASS)
X <- matrix(c(3.8642, 14.2, 41.7832, 0.4084, 9.6604, 35.5, 104.458, 1.021),ncol = 2)
y <- matrix(c(1.23,4.52,13.3,0.13),ncol = 1)

> ginv(X)%*%y
            [,1]
[1,] 0.003661243
[2,] 0.125859408

在 MATLAB/OCTAVE 中使用 pinv() 或 \

X = [3.8642    9.6604;
   14.2000   35.5000;
   41.7832  104.4580;
    0.4084    1.0210];

y = [1.2300
    4.5200
   13.3000
    0.1300];

>> X\y
ans =

   0.0036612
   0.1258594

>> pinv(X)*y
ans =

   0.0036612
   0.1258594

【讨论】：