Инверс обратного произведения матрицы amxn

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)

Я проверил Википедию и узнал, что 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

Если столбцы 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

If 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

Answer 2

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

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

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

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