m x n 矩阵的叉积的逆答案

【问题标题】：The inverse of the cross product of a m x n matrixm x n 矩阵的叉积的逆
【发布时间】：2019-02-16 09:16:03
【问题描述】：

我想用随机矩阵重新创建以下计算：

我从以下开始，结果如下：

kmin1 <- cbind(1:10,1:10,6:15,1:10,1:10,6:15,1:10,1:10,6:15)
C <- cbind(1, kmin1)                # Column of 1s
diag(C) <- 1
Ccrosprod <-crossprod(C)            # C'C
Ctranspose <- t(C)                  # C'
CCtransposeinv <- solve(Ccrosprod)  # (C'C)^-1
W <- Ctranspose %*% CCtransposeinv  # W=(C'C)^-1*C'

然而，我的假设是 C 应该能够成为 m x n 矩阵，因为没有充分的理由假设因子等于观察值。

编辑：根据 Hong Ooi 的评论，我将 kmin1 <- matrix(rexp(200, rate=.1), ncol=20) 更改为 kmin1 <- matrix(rexp(200, rate=.1), nrow=20)

我检查了wikipedia 并了解到m x n 可能有左或右反转。为了将其付诸实践，我尝试了以下方法：

kmin1 <- matrix(rexp(200, rate=.1), nrow=20)
C <- cbind(1, kmin1)                # Column of 1s
Ccrosprod <-crossprod(C)            # C'C
Ctranspose <- t(C)                  # C'
CCtransposeinv <- solve(Ccrosprod)  # (C'C)^-1
W <- Ctranspose %*% CCtransposeinv  # W=(C'C)^-1*C'

编辑：根据下面的问题，一切正常。

如果我确定这与语法没有任何关系，我会将其发布在 stackexchange 上，但由于我对矩阵没有经验，所以我不确定。

【问题讨论】：

要求n > m是普通最小二乘回归的基本假设之一/
感谢您的评论。我不认为这是矩阵代数的问题。我已经编辑了我的问题。错误现在已向下移动一个命令：Error in Ctranspose %*% CCtransposeinv : non-conformable arguments。有什么想法吗？
应该是W <- CCtransposeinv %*% Ctranspose。
谢谢，我意识到，应该对此发表评论！

标签： r matrix inverse cross-product

【解决方案1】：

如果 C 的列是线性独立的，则 C'C 是可逆的并且 (C'C)^-1C' 等于以下任何一个：

set.seed(123)
kmin1 <- matrix(rexp(200, rate=.1), nrow=20)
C <- cbind(1, kmin1)

r1 <- solve(crossprod(C), t(C))
r2 <- qr.solve(crossprod(C), t(C))

r3 <- chol2inv(chol(crossprod(C))) %*% t(C)

r4 <- with(svd(C), v %*% diag(1/d) %*% t(u))
r5 <- with(eigen(crossprod(C)), vectors %*% diag(1/values) %*% t(vectors)) %*% t(C)

r6 <- coef(lm.fit(C, diag(nrow(C))))

# check

all.equal(r1, r2)
## [1] TRUE

all.equal(r1, r3)
## [1] TRUE

all.equal(r1, r4)
## [1] TRUE

all.equal(r1, r5)
## [1] TRUE

dimnames(r6) <- NULL
all.equal(r1, r6)
## [1] TRUE

如果 C'C 不一定是可逆的，那么答案不一定是唯一的（尽管如果我们对 C(C'C)^-C' 感兴趣，那么即使伪逆也是唯一的C'C 的可能不是）。无论如何，我们可以通过取奇异值分解（或特征值分解）并使用奇异值（或特征值）的倒数并对接近 0 的值使用 0 来形成一个伪逆。这相当于使用 Moore彭罗斯伪逆。（上面显示的 lm.fit 方法也可以，但会在结果中产生一些 NA。）

set.seed(123)
kmin1 <- matrix(rexp(200, rate=.1), nrow=20)
C <- cbind(1, kmin1)
C[, 11] <- C[, 2] + C[, 3] # force singularity

eps <- 1.e-5 
s1 <- with(svd(C), v %*% diag(ifelse(abs(d) < eps, 0, 1/(d))) %*% t(u))
s2 <- with(eigen(crossprod(C)), 
  vectors %*% diag(ifelse(abs(values) < eps, 0, 1/values)) %*% t(vectors)) %*% t(C)

# check
all.equal(s1, s2)
## [1] TRUE

【讨论】：

【解决方案2】：

首先，我不熟悉您的研究/工作领域（计量经济学？），所以我不确定从特定领域知识的角度来看以下内容是否合理。

除此之外，库MASS 允许计算非方阵的Moore-Penrose generalised inverse。

因此，将您的计算推广到非方阵可能看起来像

library(MASS)
W <- ginv(t(C) %*% C) %*% t(C)

【讨论】：

尊敬的 Maurits，非常感谢您的评论，我会记住的！我自己实际上并不熟悉确切的主题，但它确实与计量经济学有关。