Anche se OP non è attivo da più di un anno, decido ancora di inviare una risposta. Vorrei usare X invece di S, come nelle statistiche, spesso vogliamo matrice di proiezione nel contesto di regressione lineare, dove X è la matrice del modello, è il vettore di risposta, mentre H = X(X'X)^{-1}X' è matrice di cappi/proiezioni in modo che Hy fornisca valori predittivi.
Questa risposta presuppone il contesto dei minimi quadrati ordinari. Per i minimi quadrati pesati, vedere Get hat matrix from QR decomposition for weighted least square regression.
Una panoramica
solve si basa sulla fattorizzazione LU di una matrice quadrata generale. Per X'X (deve essere calcolato da crossprod(X) anziché da t(X) %*% X in R, leggi ?crossprod per ulteriori informazioni), che è simmetrico, possiamo usare chol2inv basato sulla fattorizzazione di Choleksy.
Tuttavia, la fattorizzazione triangolare è meno stabile del fattorizzazione QR. Questo non è difficile da capire. Se X ha il numero condizionale kappa, X'X avrà il numero condizionale kappa^2. Questo può causare grosse difficoltà numeriche. Il messaggio di errore si ottiene:
# system is computationally singular: reciprocal condition number = 2.26005e-28
sta dicendo questo.kappa^2 è circa e-28, molto più piccolo della precisione della macchina intorno a e-16. Con la tolleranza tol = .Machine$double.eps, il valore X'X verrà considerato insufficiente, pertanto la fattorizzazione LU e Cholesky si interromperanno.
In genere, in questa situazione passiamo a SVD o QR, ma imperniato a La fattorizzazione di Cholesky è un'altra scelta.
- SVD è il metodo più stabile, ma troppo costoso;
- QR è soddisfacentemente stabile, a costi computazionali moderati ed è comunemente utilizzato nella pratica;
- Pivoted Cholesky è veloce, con stabilità accettabile. Per la matrice grande questo è preferito.
Di seguito, spiegherò tutti e tre i metodi.
Uso fattorizzazione QR
noti che la matrice di proiezione è permutazione indipendente, cioè, non importa se eseguiamo fattorizzazione QR con o senza basculante.
In R, qr.default è possibile chiamare la routine LINPACK DQRDC per la scomposizione QR non imperniata e la routine LAPACK DGEQP3 per la fattorizzazione QR imperniata a blocchi. Facciamo generare una matrice giocattolo e testare entrambe le opzioni:
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"
Annotare il $pivot è diverso per i due oggetti.
Ora, definiamo una funzione wrapper per calcolare QQ':
f <- function (QR) {
## thin Q-factor
Q <- qr.qy(QR, diag(1, nrow = nrow(QR$qr), ncol = QR$rank))
## QQ'
tcrossprod(Q)
}
vedremo che qr_linpack e qr_lapack invia la stessa matrice di proiezione:
H1 <- f(qr_linpack)
H2 <- f(qr_lapack)
mean(abs(H1 - H2))
# [1] 9.530571e-17
Uso decomposizione in valori singolari
In R, svd calcola la decomposizione del valore singolare. Usiamo ancora l'esempio matrice sopra 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
Nuovamente, si ottiene la stessa matrice di proiezione.
Uso pivot Cholesky
Per dimostrazione, usiamo ancora l'esempio X sopra.
## 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
Ancora una volta, otteniamo la stessa matrice di proiezione.
C'è un motivo per cui non si desidera utilizzare princomp o prcomp? Il calcolo dei componenti principali non deve essere eseguito a mano. –
Paura che non ci sia una soluzione generale, se questo è il numero della condizione si ha un problema. –
Ciao, non sto provando a fare PCA tanto quanto implementare uno stimatore di cui ho letto. Trovo strano che non riesca a farlo funzionare anche per una semplice matrice di strumenti fittizi, quando sembra garantito che non sia collineare. – bikeclub