Инверс обратного произведения матрицы a 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'

Однако я предполагаю, что С должен быть в состоянии m x n матрица, так как нет веских оснований предполагать, что факторы равны наблюдениям.

РЕДАКТИРОВАТЬ: На основании комментария Hong Ooi я изменил kmin1 <- matrix(rexp(200, rate=.1), ncol=20) в kmin1 <- matrix(rexp(200, rate=.1), nrow=20)

Я проверил википедию и узнал, что 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, если бы был уверен, что это не имеет ничего общего с синтаксисом, но так как у меня нет опыта работы с матрицами, я не уверен.

Answer 1

Прежде всего, я не знаком с вашей областью исследований / работ (эконометрика?), Поэтому я не уверен, имеет ли смысл следующее с точки зрения предметной области знаний.

Что в стороне, библиотека MASS позволяет вычислить обобщенную инверсию Мура-Пенроуза неквадратной матрицы.

Таким образом, потенциальное обобщение ваших расчетов на неквадратные матрицы может выглядеть следующим образом

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

Answer 2

Если столбцы C линейно независимы, то C'C обратимо и (C'C)^-1 C 'равно любому из них:

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. Это эквивалентно использованию Мура Псевдообратный Пенроуз. (Подход 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