Aunque OP no ha estado activo durante más de un año, todavía decido publicar una respuesta. Me gustaría utilizar en lugar de XS, al igual que en las estadísticas, muchas veces queremos matriz de proyección en el contexto de regresión lineal, donde X es la matriz del modelo, y es el vector de respuesta, mientras que H = X(X'X)^{-1}X' es la matriz sombrero/proyección de manera que Hy da valores predictivos.
Esta respuesta asume el contexto de mínimos cuadrados ordinarios. Para mínimos cuadrados ponderados, vea Get hat matrix from QR decomposition for weighted least square regression.
Una visión general
solve se basa en factorización LU de una matriz cuadrada general. Para X'X (debe calcularse por crossprod(X) en lugar de t(X) %*% X en R, lea ?crossprod para obtener más información) que es simétrico, podemos usar chol2inv que se basa en la factorización de Choleksy.
Sin embargo, la factorización triangular es menos estable que la factorización QR. Esto no es difícil de entender. Si X tiene el número condicional kappa, X'X tendrá el número condicional kappa^2. Esto puede causar una gran dificultad numérica. El mensaje de error que aparece:
# system is computationally singular: reciprocal condition number = 2.26005e-28
es solo decir esto.kappa^2 es aproximadamente e-28, mucho más pequeño que la precisión de la máquina en alrededor de e-16. Con la tolerancia tol = .Machine$double.eps, X'X se considerará como rango deficiente, por lo que las factorías LU y Cholesky se descompondrán.
Generalmente, cambiamos a SVD o QR en esta situación, pero pivotó La factorización de Cholesky es otra opción.
- SVD es el método más estable, pero demasiado caro;
- QR es satisfactoriamente estable, a costos computacionales moderados, y se usa comúnmente en la práctica;
- Pivoted Cholesky es rápido, con una estabilidad aceptable. Para una matriz grande, esta es la preferida.
A continuación, explicaré los tres métodos.
Uso de factorización QR
Tenga en cuenta que la matriz de proyección es de permutación independiente, es decir, no importa si realizamos factorización QR con o sin pivotar.
En R, qr.default puede llamar a la rutina LINPACK DQRDC para la factorización QR no pivotada, y la rutina LAPACK DGEQP3 para la factorización QR pivotante de bloque. Vamos a generar una matriz de juguete y probar las dos opciones:
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"
Nota del $pivot es diferente para los dos objetos.
Ahora, se define una función de contenedor para calcular QQ':
f <- function (QR) {
## thin Q-factor
Q <- qr.qy(QR, diag(1, nrow = nrow(QR$qr), ncol = QR$rank))
## QQ'
tcrossprod(Q)
}
veremos que qr_linpack y qr_lapack dan la misma matriz de proyección:
H1 <- f(qr_linpack)
H2 <- f(qr_lapack)
mean(abs(H1 - H2))
# [1] 9.530571e-17
El uso de la descomposición en valores singulares
En R, svd calcula la descomposición en valores singulares. Todavía usamos el ejemplo de matriz anterior 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
Una vez más, se obtiene la misma matriz de proyección.
Usando Pivotado factorización de Cholesky
Para la demostración, todavía usamos el ejemplo anterior 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
De nuevo, obtenemos la misma matriz de proyección.
¿Hay algún motivo por el que no desee utilizar princomp o prcomp? El cálculo de los componentes principales no necesita hacerse a mano. –
Miedo a que no haya una solución general, si ese es el número de condición, tiene un problema. –
Hola, no estoy tratando de hacer PCA tanto como implementar un estimador sobre el que he leído. Me parece extraño que no pueda hacer que esto funcione, incluso para una matriz simple de instrumentos ficticios, cuando parece que eso no es colineal. – bikeclub