用R解决非方程线性系统

Question

如何用R解决非方程线性系统：A X = B？

（在系统无解或有无限多个解的情况下）

示例：

A=matrix(c(0,1,-2,3,5,-3,1,-2,5,-2,-1,1),3,4,T)
B=matrix(c(-17,28,11),3,1,T)

A
     [,1] [,2] [,3] [,4]
[1,]    0    1   -2    3
[2,]    5   -3    1   -2
[3,]    5   -2   -1    1


B
     [,1]
[1,]  -17
[2,]   28
[3,]   11

- Basj

---

请查看此链接。一个好的可重现的例子将会帮助其他人更容易地解决你的问题。 - CHP

我添加了一个可重现的例子。 - Basj

在solve的文档文件中，他们提到“qr.solve可以处理非方系统。”所以当我使用这里给出的A和B，并尝试qr.solve(A,B)时，我会收到错误消息Error in qr.solve(A, B) : singular matrix 'a' in solve。有什么想法吗？ - Gimelist

Answer 1

如果矩阵A的行数多于列数，则应使用最小二乘法拟合。

如果矩阵A的行数少于列数，则应执行奇异值分解。每种算法都尽力使用假设来给出解决方案。

以下链接显示如何使用SVD作为求解器：

http://www.ecse.rpi.edu/~qji/CV/svd_review.pdf

我来将其应用到你的问题上，看看是否有效：

您的输入矩阵 A 和已知的右手边向量 B：

> A=matrix(c(0,1,-2,3,5,-3,1,-2,5,-2,-1,1),3,4,T)
> B=matrix(c(-17,28,11),3,1,T)
> A
     [,1] [,2] [,3] [,4]
[1,]    0    1   -2    3
[2,]    5   -3    1   -2
[3,]    5   -2   -1    1
> B
     [,1]
[1,]  -17
[2,]   28
[3,]   11

让我们分解您的A矩阵：

> asvd = svd(A)
> asvd
$d
[1] 8.007081e+00 4.459446e+00 4.022656e-16

$u
           [,1]       [,2]       [,3]
[1,] -0.1295469 -0.8061540  0.5773503
[2,]  0.7629233  0.2908861  0.5773503
[3,]  0.6333764 -0.5152679 -0.5773503

$v
            [,1]       [,2]       [,3]
[1,]  0.87191556 -0.2515803 -0.1764323
[2,] -0.46022634 -0.1453716 -0.4694190
[3,]  0.04853711  0.5423235  0.6394484
[4,] -0.15999723 -0.7883272  0.5827720

> adiag = diag(1/asvd$d)
> adiag
          [,1]      [,2]        [,3]
[1,] 0.1248895 0.0000000 0.00000e+00
[2,] 0.0000000 0.2242431 0.00000e+00
[3,] 0.0000000 0.0000000 2.48592e+15

这是关键：在d中的第三个特征值非常小；相反，在adiag中的对角元素非常大。 在解决之前，将其设置为零：

> adiag[3,3] = 0
> adiag
          [,1]      [,2] [,3]
[1,] 0.1248895 0.0000000    0
[2,] 0.0000000 0.2242431    0
[3,] 0.0000000 0.0000000    0

现在让我们计算解决方案（请参见我上面给你的链接中的第16张幻灯片）：

> solution = asvd$v %*% adiag %*% t(asvd$u) %*% B
> solution
          [,1]
[1,]  2.411765
[2,] -2.282353
[3,]  2.152941
[4,] -3.470588

现在我们有了一个解决方案，让我们将其替换回去，看看是否给我们相同的。

> check = A %*% solution
> check
     [,1]
[1,]  -17
[2,]   28
[3,]   11

这是你开始的B面，所以我认为我们做得很好。

以下是AMS提供的另一个不错的奇异值分解讨论：

http://www.ams.org/samplings/feature-column/fcarc-svd

- duffymo

---

既然我已经给你了相关的名称，也许你可以先进行一些研究并自己尝试一下。 - duffymo

我做了最小二乘法和奇异值分解。请查看R中的lm()和svd()方法。 - duffymo

我已经在我的答案中添加了SVD解决方案。 - duffymo

我认为这是因为线性求解器函数对于非方阵没有意义。就像我在原始答案中所说的，如果m>n，则使用最小二乘法；如果m<n，则使用SVD。就我所知，就是这样。 - duffymo

你正在构建的矩阵 asvd$v %*% adiag %*% t(asvd$u) 是 A 的伪逆，是吗？ - Basj

@Basj 是的，它就是摩尔-彭罗斯伪逆 - 你也可以使用 Mass::ginv 和 pracma::pinv 进行计算。 - Tom Wenseleers

Answer 2

旨在解决 Ax = b

其中 A 是 p 乘以 q，x 是 q 乘以 1，且 b 是 p 乘以 1，给定 A 和 b 求 x。

方法一：广义逆矩阵：Moore-Penrose https://en.wikipedia.org/wiki/Generalized_inverse

将等式两边乘以A'，我们得到：

A'Ax = A' b

这里的A'是A的转置。请注意，A'A现在是一个 q 乘以 q 的矩阵。 解决此方程的一种方法是将等式两边乘以 A'A 的逆矩阵。如下所示：

x = (A'A)^{-1} A' b

这就是广义逆矩阵的理论依据。 这里的 G = (A'A)^{-1} A' 是 A 的伪逆矩阵。

library(MASS)

ginv(A) %*% B

#          [,1]
#[1,]  2.411765
#[2,] -2.282353
#[3,]  2.152941
#[4,] -3.470588

方法2：使用SVD求广义逆矩阵。

@duffymo使用SVD获得A的伪逆矩阵。

然而，svd(A)$d的最后几个元素可能不像这个例子中那么小。因此，不能直接使用该方法。以下是一个示例，其中A的最后一个元素都不接近于零。

A <- as.matrix(iris[11:13, -5])    
A
#   Sepal.Length Sepal.Width Petal.Length Petal.Width
# 11          5.4         3.7          1.5         0.2
# 12          4.8         3.4          1.6         0.2
# 13          4.8         3.0          1.4         0.1

svd(A)$d
# [1] 10.7820526  0.2630862  0.1677126

有一个选项是查看cor（A）中的奇异值。

svd(cor(A))$d
# [1] 2.904194e+00 1.095806e+00 1.876146e-16 1.155796e-17

现在，很明显只有两个较大的奇异值存在。因此，现在可以对A应用svd来获取伪逆，如下所示。

svda <- svd(A)
G = svda$v[, 1:2] %*% diag(1/svda$d[1:2]) %*% t(svda$u[, 1:2])
# to get x
G %*% B

- Suren

---

它可以在CRAN上找到。 - kangaroo_cliff

我认为你的方法1和方法2是相同的，因为MASS::ginv（以及pracma::pinv）使用与你的方法2完全相同的代码来计算Moore-Penrose伪逆（即通过SVD）。 - Tom Wenseleers

这有点正确，因为我使用了 MASS::ginv 来演示这种方法。我会在以后编辑答案并提供不同的实现方式。干杯！ - kangaroo_cliff