用技术发现美 | R 解方程、做矩阵分解

解方程

解矩阵方程 Ax = b 时，使用 solve(A, b) 可以解出 x。例如求解《孙子算经》中著名的鸡兔同笼问题：“今有雉（鸡）兔同笼，上有三十五头，下有九十四足。问雉兔各几何。” 列方程:

1 鸡头 * 鸡只数 + 1 兔头 * 兔只数 = 35
2 鸡脚 * 鸡只数 + 4 兔脚 * 兔只数 = 94

即

$ \begin{bmatrix} 1 & 1 \\ 2 & 4 \end{bmatrix} \begin{bmatrix} 鸡 \\ 兔 \end{bmatrix} = \begin{bmatrix} 35 \\ 94 \end{bmatrix} $

使用 R 来求解:

> A <- matrix(c(1,2,1,4),2)
> A
     [,1] [,2]
[1,]    1    1
[2,]    2    4
> b <- c(35, 94)
> solve(A, b)
[1] 23 12

即得到答案为 23 只鸡，12 只兔。

可以使用同一方法来求矩阵 A 的逆: solve(A) 。不过，在很多情况下，并不需要真正去计算一个矩阵的逆。

计算二次型 $ x A^{-1} x^T $ 时，可以无需求 A 的逆，使用如下方法会更快:

x %*% solve(A,x)

使用 eigen(Sm) 可以计算对称阵 Sm 的特征值，并得到列表，其中含有 val 和 vec 。可以指定只计算特征值而不计算特征向量:

eigen(Sm, only.values=TRUE)$values

使用 svd(M) 可以计算任意矩阵的奇异值分解，得到 u,v,d。 u 是和 M 具有相同列空间的正交列矩阵U， v 是和 M 具有相同行空间的正交列矩阵V， d 是一个以向量形式返回的正数值对角阵，是 M 的奇异值。且有 $ M = U \times D \times V^T $ 。

如果 M 是方阵，则它的行列式的绝对值为 prod(svd(M)$d) 。它的迹为 tr()，即 sum(svd(M)$d) 。

使用 det() 可以计算方阵的行列式，并且包含了符号。使用 determinant() 可以计算行列式的模数（即绝对值的对数）和符号。不过，一般情况下为了解决实际问题，通常是不需要真正计算矩阵的行列式的。

最小二乘法为了拟合 y = X b + e 中的 b。 lsfit() 返回一个列表包含了最小二乘拟合的结果:

ans <- lsfit(X, y)

参见手册。并且需要注意的是，在解决实际问题时，更应该使用 lm() 来代替 lsfit() 进行回归建模。对于非线性回归，使用 nls() ，参见R - 数据拟合。

qr() 也是与最小二乘紧密相关的另一个函数。

 Xplus <- qr(X)
 b <- qr.coef(Xplus, y)
 fit <- qr.fitted(Xplus, y)
 res <- qr.resid(Xplus, y)

如果一个方程组无法用 solve() 来求解，也就是 A 不是满秩矩阵时，无法求得 A 的逆阵，这时可以使用 qr.solve() 来求最小二乘解:

qr.solve(A, B)

也可以使用 svd 来求最小二乘解，比 qr() 更加稳定:

s <- svd(A)
x <- s$v %*% diag(1/s$d) %*% t(s$u) %*% B

需要注意的是：在做 svd 时，如果有特别接近 0 的奇异值，则说明 A 中的某些列很有可能线性相关。这时可以简单地把相应的奇异值设为 0，再求伪逆 a，可得拟合解 x。

例如 s$d[4] 之后接近 0

s$d[4:5] <- 0
s$u[,4:5] <- 0
s$v[4:5,] <- 0
x <- s$v %*% diag(c(1/s$d[1:3],0,0)) %*% t(s$u) %*% B