在R中计算非方阵的逆矩阵

Question

我对 R 语言还很陌生，想知道如何计算非方阵的矩阵逆（非方阵？非规则矩阵？我不确定正确术语）。从我的书和快速 Google 搜索中发现（参见来源），如果 a 是方阵，则可以使用 solve(a) 找到矩阵的逆。但我创建的矩阵不是方阵，据我的理解：

  > matX<-matrix(c(rep(1, 8),2,3,4,0,6,4,3,7,-2,-4,3,5,7,8,9,11), nrow=8, ncol=3);
  > matX

       [,1] [,2] [,3]
  [1,]    1    2   -2
  [2,]    1    3   -4
  [3,]    1    4    3
  [4,]    1    0    5
  [5,]    1    6    7
  [6,]    1    4    8
  [7,]    1    3    9
  [8,]    1    7   11
  > 

有没有一种解决这个大小的矩阵的函数，还是我需要对每个元素进行操作？因为solve()函数会报错：

  Error in solve.default(matX) : 'a' (8 x 3) must be square

我正在尝试从上述矩阵中获得的计算是：(matX'matX)^-1。提前致谢。

- Reanimation

---

关于你最后一行的代码，请尝试使用：solve(t(matX) %*% matX) - Mark Heckmann

Answer 1

ginv函数在MASS包中可以给出矩阵的广义逆矩阵。将原始矩阵与其进行左乘，将会得到单位矩阵：

library(MASS)
inv <- ginv(matX)

# test it out
inv %*% matX
##               [,1]         [,2]          [,3]
## [1,]  1.000000e+00 6.661338e-16  4.440892e-15
## [2,] -8.326673e-17 1.000000e+00 -1.110223e-15
## [3,]  6.938894e-17 8.326673e-17  1.000000e+00

正如评论中建议的那样，使用zapsmall可以以更好的方式显示此内容：

zapsmall(inv %*% matX)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

< p > matX'matX的逆矩阵如下：

tcrossprod(inv)
##              [,1]         [,2]         [,3]
## [1,]  0.513763935 -0.104219636 -0.002371406
## [2,] -0.104219636  0.038700372 -0.007798748
## [3,] -0.002371406 -0.007798748  0.006625269

解决方法，然而，如果你的目标是计算matX'matX的逆矩阵，那么你根本不需要它。这与MASS无关：

solve(crossprod(matX))
##              [,1]         [,2]         [,3]
## [1,]  0.513763935 -0.104219636 -0.002371406
## [2,] -0.104219636  0.038700372 -0.007798748
## [3,] -0.002371406 -0.007798748  0.006625269

奇异值分解（SVD）也可用于类似的操作，并且不需要使用MASS：

with(svd(matX), v %*% diag(1/d^2) %*% t(v))
##              [,1]         [,2]         [,3]
## [1,]  0.513763935 -0.104219636 -0.002371406
## [2,] -0.104219636  0.038700372 -0.007798748
## [3,] -0.002371406 -0.007798748  0.006625269

ADDED some additional info.

- G. Grothendieck

---

zapsmall 可能会很好地清晰地显示 inv %*% matX 的结果... - Ben Bolker

哦，太酷了:D感谢您抽出时间回答:D太棒了。 - Reanimation

对于SVD解法，为什么要将对角矩阵平方？ diag(1/d^2) - Katelynn ruan

@user1143669，请阅读问题的倒数第二行。 - G. Grothendieck

Answer 2

您可以进行所谓的“Moore-Penrose伪逆”。这里有一个名为exp.mat的函数可以为您完成此操作。还有一个例子在这里概述了它的用法。

exp.mat():

#The exp.mat function performs can calculate the pseudoinverse of a matrix (EXP=-1)
#and other exponents of matrices, such as square roots (EXP=0.5) or square root of 
#its inverse (EXP=-0.5). 
#The function arguments are a matrix (MAT), an exponent (EXP), and a tolerance
#level for non-zero singular values.
exp.mat<-function(MAT, EXP, tol=NULL){
    MAT <- as.matrix(MAT)
    matdim <- dim(MAT)
    if(is.null(tol)){
        tol=min(1e-7, .Machine$double.eps*max(matdim)*max(MAT))
    }
    if(matdim[1]>=matdim[2]){ 
        svd1 <- svd(MAT)
        keep <- which(svd1$d > tol)
        res <- t(svd1$u[,keep]%*%diag(svd1$d[keep]^EXP, nrow=length(keep))%*%t(svd1$v[,keep]))
    }
    if(matdim[1]<matdim[2]){ 
        svd1 <- svd(t(MAT))
        keep <- which(svd1$d > tol)
        res <- svd1$u[,keep]%*%diag(svd1$d[keep]^EXP, nrow=length(keep))%*%t(svd1$v[,keep])
    }
    return(res)
}

使用示例：

source("exp.mat.R")
X <- matrix(c(rep(1, 8),2,3,4,0,6,4,3,7,-2,-4,3,5,7,8,9,11), nrow=8, ncol=3)
iX <- exp.mat(X, -1)
zapsmall(iX %*% X) # results in identity matrix
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

- Marc in the box

---

你好。感谢您的时间。我需要使用特定的软件包才能使用上面提到的函数吗？ - Reanimation

哦，等等，我猜你需要在你的程序中自己创建那个函数？ - Reanimation

@Reanimation - 不用谢。你可以直接加载该函数（即将函数代码传递到你的R环境中）。如果你更喜欢加载一个包而不是使用source函数，那么在corpcor包中也可以找到类似的函数(pseudoinverse)。exp.mat函数的优点是它还可以接受其他指数值（例如EXP=-0.5）。 - Marc in the box

MASS::ginv 还实现了一个伪逆，可以在标准的 R 安装中使用... - Ben Bolker

@BenBolker - 谢谢。我简直不敢相信我之前不知道ginv，感谢你指出来。 - Marc in the box