m x n 矩阵的叉积的倒数

Question

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

我从以下开始，给出了结果：

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'

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

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

Answer 1

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

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

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

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

Answer 2

如果 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

m x n 矩阵的叉积的倒数

The inverse of the cross product of a m x n matrix

r

matrix

inverse

cross-product