Вычислить матрицу проекции / шапки с помощью QR-факторизации, SVD (и факторизации Холецкого?) - PullRequest
Новая
       41

Вычислить матрицу проекции / шапки с помощью QR-факторизации, SVD (и факторизации Холецкого?)

9 голосов
/ 31 января 2012

Я пытаюсь вычислить в R матрицу проекции P произвольной матрицы N x J S: 

P = S (S'S) ^ -1 S'

Я пытался выполнить это с помощью следующей функции: 

P <- function(S){
  output <- S %*% solve(t(S) %*% S) %*% t(S)
  return(output)
}

Но когда я использую это, я получаю ошибки, которые выглядят так: 

# Error in solve.default(t(S) %*% S, t(S), tol = 1e-07) : 
#  system is computationally singular: reciprocal condition number = 2.26005e-28

Я думаю, что это результат численного занижения и / или нестабильности, о чем говорилось во многих местах, таких как r-help и здесь , но я не достаточно опытен в использовании SVD или QR-разложение для решения проблемы, или же приведите этот существующий код в действие Я также попробовал предложенный код, который должен написать решите как система: 

output <- S %*% solve (t(S) %*% S, t(S), tol=1e-7)

Но все равно это не работает. Любые предложения будут оценены.

Я почти уверен, что моя матрица должна быть обратимой и не иметь какой-либо коллинеарности хотя бы потому, что я пытался проверить это с помощью матрицы ортогональных фиктивных переменных, и она все еще не работает.

Кроме того, я хотел бы применить это к довольно большим матрицам, поэтому я ищу аккуратное общее решение.

1 Ответ

10 голосов
/ 02 сентября 2016

Хотя ОП не был активен более года, я все же решил опубликовать ответ.Я бы использовал X вместо S, как в статистике, мы часто хотим проекционную матрицу в контексте линейной регрессии, где X - матрица модели, y - вектор ответа, а H = X(X'X)^{-1}X' - шляпа /матрица проекции, так что Hy дает прогнозные значения.

Этот ответ принимает контекст обычных наименьших квадратов.Для взвешенных наименьших квадратов см. Получить матрицу шапки из разложения QR для взвешенной регрессии наименьших квадратов .

Обзор

solveоснован на факторизации LU общей квадратной матрицы.Для X'X (должно быть вычислено как crossprod(X), а не t(X) %*% X в R, читайте ?crossprod для более), который является симметричным, мы можем использовать chol2inv, который основан на факторизации Холекса.

Однако треугольная факторизация менее стабильна, чем QR факторизация.Это не сложно понять.Если X имеет условный номер kappa, X'X будет иметь условный номер kappa ^ 2.Это может вызвать большие численные трудности.Сообщение об ошибке: 

# system is computationally singular: reciprocal condition number = 2.26005e-28

просто говорит об этом.kappa ^ 2 составляет около e-28, намного меньше точности станка на отметке e-16.С допуском tol = .Machine$double.eps, X'X будет рассматриваться как недостаток ранга, поэтому LU и коэффициент Холецкого будут нарушены.

Обычно в этой ситуации мы переключаемся на SVD или QR, но pivoted Факторизация Холецкого - это еще один выбор.

  • SVD является наиболее стабильным методом, но слишком дорогим;
  • QR является достаточно стабильным, при умеренных вычислительных затратах и ​​обычно используется на практике;
  • Поворотный Холецкий быстр, с приемлемой стабильностью.Для большой матрицы этот вариант предпочтителен.

Далее я объясню все три метода.

Использование QR-факторизации

QR factorization

Обратите внимание, что матрица проекции не зависит от перестановки, т. Е. Не имеет значения, выполняем ли мы QR-факторизацию с поворотом или без него.

In R, qr.default может вызывать процедуру LINPACK DQRDC для неповоротной QR-факторизации, а подпрограмму LAPACK DGEQP3 для блочной развертки QR-факторизации.Давайте сгенерируем игрушечную матрицу и протестируем оба варианта: 

set.seed(0); X <- matrix(rnorm(50), 10, 5)
qr_linpack <- qr.default(X)
qr_lapack <- qr.default(X, LAPACK = TRUE)

str(qr_linpack)
# List of 4
# $ qr   : num [1:10, 1:5] -3.79 -0.0861 0.3509 0.3357 0.1094 ...
# $ rank : int 5
# $ qraux: num [1:5] 1.33 1.37 1.03 1.01 1.15
# $ pivot: int [1:5] 1 2 3 4 5
# - attr(*, "class")= chr "qr"

str(qr_lapack)
# List of 4
# $ qr   : num [1:10, 1:5] -3.79 -0.0646 0.2632 0.2518 0.0821 ...
# $ rank : int 5
# $ qraux: num [1:5] 1.33 1.21 1.56 1.36 1.09
# $ pivot: int [1:5] 1 5 2 4 3
# - attr(*, "useLAPACK")= logi TRUE
# - attr(*, "class")= chr "qr"

Обратите внимание, что $pivot отличается для двух объектов.

Теперь мы определим функцию-оболочку для вычисления QQ': 

f <- function (QR) {
  ## thin Q-factor
  Q <- qr.qy(QR, diag(1, nrow = nrow(QR$qr), ncol = QR$rank))
  ## QQ'
  tcrossprod(Q)
  }

Мы увидим, что qr_linpack и qr_lapack дают одинаковую матрицу проекции: 

H1 <- f(qr_linpack)
H2 <- f(qr_lapack)

mean(abs(H1 - H2))
# [1] 9.530571e-17

Использование разложения по сингулярным значениям

SVD

В R svd вычисляет разложение по сингулярным числам.Мы по-прежнему используем приведенный выше пример матрицы X: 

SVD <- svd(X)

str(SVD)
# List of 3
# $ d: num [1:5] 4.321 3.667 2.158 1.904 0.876
# $ u: num [1:10, 1:5] -0.4108 -0.0646 -0.2643 -0.1734 0.1007 ...
# $ v: num [1:5, 1:5] -0.766 0.164 0.176 0.383 -0.457 ...

H3 <- tcrossprod(SVD$u)

mean(abs(H1 - H3))
# [1] 1.311668e-16

Опять же, мы получаем ту же матрицу проекций.

Использование факторизации Pivoted-Cholesky

Pivoted Cholesky factorization

Для демонстрации мы все еще используем приведенный выше пример X. 

## pivoted Chol for `X'X`; we want lower triangular factor `L = R'`:
## we also suppress possible rank-deficient warnings (no harm at all!)
L <- t(suppressWarnings(chol(crossprod(X), pivot = TRUE)))

str(L)
# num [1:5, 1:5] 3.79 0.552 -0.82 -1.179 -0.182 ...
# - attr(*, "pivot")= int [1:5] 1 5 2 4 3
# - attr(*, "rank")= int 5

## compute `Q'`
r <- attr(L, "rank")
piv <- attr(L, "pivot")
Qt <- forwardsolve(L, t(X[, piv]), r)

## P = QQ'
H4 <- crossprod(Qt)

## compare
mean(abs(H1 - H4))
# [1] 6.983997e-17

Опять же, мы получаем ту же матрицу проекции.

...