Mimo że OP nie był aktywny przez ponad rok, nadal decyduję się na opublikowanie odpowiedzi. Używałbym X zamiast S, ponieważ w statystykach często potrzebujemy matrycy projekcyjnej w kontekście regresji liniowej, gdzie X jest matrycą modelu, y jest wektorem odpowiedzi, podczas gdy H = X(X'X)^{-1}X' jest macierzą kapelusz/rzutowanie, dzięki czemu Hy podaje wartości predykcyjne.
Ta odpowiedź przyjmuje kontekst zwykłych najmniejszych kwadratów. W przypadku ważonych najmniejszych kwadratów patrz Get hat matrix from QR decomposition for weighted least square regression.
Przegląd
solve opiera się na LU faktoryzacji ogólny kwadratowy podłoża. Dla X'X (powinno być obliczone przez crossprod(X) zamiast t(X) %*% X w R, przeczytaj ?crossprod po więcej), co jest symetryczne, możemy użyć chol2inv, który jest oparty na faktoryzacji Choleksy'ego.
Jednak rozkład trójkątny jest mniej stabilny niż faktoryzacja QR. To nie jest trudne do zrozumienia. Jeśli X ma numer warunkowy kappa, X'X będzie mieć numer warunkowy kappa^2. Może to spowodować duże trudności liczbowe. Komunikat o błędzie:
# system is computationally singular: reciprocal condition number = 2.26005e-28
mówi właśnie to.kappa^2 to około e-28, o wiele mniejsza niż precyzja maszyny około e-16. Z tolerancją tol = .Machine$double.eps, X'X będzie postrzegany jako niedobór rangi, a więc rozkład LU i Cholesky ulegnie rozpadowi.
Zasadniczo w tej sytuacji przełączamy się na SVD lub QR, ale obrócona Cholesky factorization to kolejny wybór.
- SVD jest najbardziej stabilną metodą, ale za drogie;
- QR jest satysfakcjonująco stabilny, przy umiarkowanych kosztach obliczeniowych i jest powszechnie stosowany w praktyce;
- Pivoted Cholesky jest szybka, z akceptowalną stabilnością. W przypadku dużej matrycy jest to preferowane.
Poniżej wyjaśnię wszystkie trzy metody.
Stosując Qr faktoryzacji
Należy zauważyć, że macierz projekcji jest permutacją niezależnie, to znaczy, że nie ma znaczenia, czy wykonać QR faktoryzacji z lub bez obracania.
W języku R, qr.default można wywołać procedurę LINPACK DQRDC dla nieindukowanej faktoryzacji QR i procedurę RAPORTU LAPACK DGEQP3 dla blokowej analizy współczynnika QR. Załóżmy wygenerować macierz zabawki i przetestować obie opcje:
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"
Zanotuj $pivot jest różna dla dwóch obiektów.
Teraz zdefiniowania funkcji opakowania w celu obliczenia QQ':
f <- function (QR) {
## thin Q-factor
Q <- qr.qy(QR, diag(1, nrow = nrow(QR$qr), ncol = QR$rank))
## QQ'
tcrossprod(Q)
}
zobaczymy że qr_linpack i qr_lapack dają taką samą matrycę do projekcji
H1 <- f(qr_linpack)
H2 <- f(qr_lapack)
mean(abs(H1 - H2))
# [1] 9.530571e-17
za pomocą rozkładu wartości pojedynczą
W R, svd oblicza dekompozycję wartości osobliwych. My nadal korzystać z powyższego przykładu macierz 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
Znów mamy tę samą macierz projekcji.
Korzystanie obracany Choleskiego faktoryzacji
Do demonstracji, my nadal korzystać z przykładu X powyżej.
## 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
Znów otrzymujemy tę samą matrycę projekcji.
Czy istnieje powód, dla którego nie chcesz używać polecenia princomp lub prcomp? Obliczanie głównych składników nie musi być wykonywane ręcznie. –
Obawiam się, że nie ma ogólnego rozwiązania, jeśli jest to numer warunku, w którym masz problem. –
Witam, nie próbuję zrobić PCA tak bardzo jak zaimplementować estymator, o którym czytałem. Uważam za dziwne, że nie mogę tego zrobić nawet dla prostej matrycy instrumentów dummy, gdy wydaje się, że to nie jest współliniowość. – bikeclub