准备
我们首先以矩阵形式A * x = b 表达您的线性系统。如果您不清楚如何执行此操作，请阅读 General forms。对于您的示例，您可以将其表示为：
## x = r_Aus, r_Chi, r_Fra, r_Ser, r_USA, r_Ven
r_Aus - r_Fra - r_Ser - r_USA = 8.7
- r_Aus - r_Chi + r_Fra - r_Ser = 2.7
- r_Aus - r_Chi + r_USA - r_Ven = 37
+ r_Chi - r_Fra - r_USA - r_Ven = -29.7
- r_Aus - r_Fra + r_Ser - r_Ven = 2.7
- r_Chi - r_Ser - r_USA + r_Ven = -21.3
然后定义系数矩阵A和RHS向量b：
A <- matrix(c( 1, 0, -1, -1, -1, 0,
-1, -1, 1, -1, 0, 0,
-1, -1, 0, 0, 1, -1,
0, 1, -1, 0, -1, -1,
-1, 0, -1, 1, 0, -1,
0, -1, 0, -1, -1, 1),
nrow = 6, ncol = 6, byrow = TRUE)
b <- as.matrix(c(8.7, 2.7, 37, -29.7, 2.7, -21.3))
正在尝试solve()
几乎总是，我们首先考虑solve。但是solve()是基于LU分解的，需要满秩系数矩阵A；当发现A 秩不足时，LU 分解满足 0 对角元素并失败。让我们试试你的A 和b：
solve(A, b)
#Error in solve.default(A, b) :
# Lapack routine dgesv: system is exactly singular: U[6,6] = 0
U[0,0] = 0 告诉你，你的A 的排名只有 5。
一种稳定的方法：QR 分解
众所周知，QR 分解是一种非常稳定的方法。我们可以使用.lm.fit() 来做到这一点：
x <- .lm.fit(A, b)
x$coef
# [1] 4.783333 -5.600000 -21.450000 -18.650000 40.866667 0.000000
x$rank
# [1] 5
您的系统等级为 5，因此执行最小二乘拟合。第 6 个值是 r_Ven 被限制为 0，并且没有一个方程完全满足。 x$resi 为您提供残差，即b - A %*% x$beta。
高斯消元法
为了完成这张图，我不得不提到高斯消元法。 理论上这是最好的方法，因为您可以确定：
- 有一个独特的解决方案；
- 没有解决办法；
- 有无数种解决方案
以及求解线性系统。
周围有一个小的 R 包optR，但我发现它做得并不完美。
#install.packages("optR")
library(optR)
?optR 给出了一个满秩线性系统作为例子，它工作得很好（简单地使用solve(A, b) 也可以工作！）。但是对于您的 rank-5 系统，它给出：
optR(A, b, method="gauss")
call:
optR.default(x = A, y = b, method = "gauss")
Coefficients:
[,1]
[1,] 9.466667
[2,] -24.333333
[3,] -16.766667
[4,] -4.600000
[5,] 22.133333
[6,] 0.000000
Warning messages:
1: In opt.matrix.reorder(A, tol) : Singular Matrix
2: In opt.matrix.reorder(A, tol) : Singular Matrix
请注意您的线性系统秩不足的警告消息。要了解optR 在这种情况下的作用，请将b 与
进行比较
A %*% x$beta
# [,1]
#[1,] 8.7
#[2,] 2.7
#[3,] 37.0
#[4,] -29.7
#[5,] 2.7
#[6,] 6.8
前 5 个方程得到满足，第 6 个除外。所以，optR 放弃了最后一个方程来解决秩不足问题，而不是进行最小二乘拟合。