10

如何用 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
4

2 回答 2

11

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

如果矩阵 A 的行数少于列数,则应执行奇异值分解。每个算法都尽其所能通过使用假设为您提供解决方案。

这是一个链接,显示如何使用 SVD 作为求解器:

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

让我们将其应用于您的问题,看看它是否有效:

您的输入矩阵A和已知的 RHS 向量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

这是关键:in的第三个特征值d非常小;相反,in的对角元素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

现在我们有了一个解决方案,让我们将其替换回去,看看它是否给我们相同的结果B

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

那是B你开始的那一面,所以我认为我们很好。

这是来自 AMS 的另一个很好的 SVD 讨论:

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

于 2013-11-04T12:13:14.680 回答
4

目的是解决Ax = b

其中Ap by qxq by 1并且bp by 1对于x给定Ab

方法 1:广义逆:Moore-Penrose https://en.wikipedia.org/wiki/Generalized_inverse

等式两边相乘,我们得到

A'Ax = A' b

其中A'是 A 的转。请注意,A'A现在是q by 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
于 2017-10-16T09:24:53.260 回答