हालांकि ओ पी एक वर्ष से अधिक के लिए सक्रिय नहीं किया गया है, मैं अभी भी एक जवाब पोस्ट करने के लिए निर्णय लेते हैं। मैं S के बजाय X का प्रयोग करेंगे, सांख्यिकी में के रूप में, हम अक्सर रेखीय प्रतीपगमन संदर्भ में, जहां X मॉडल मैट्रिक्स है में प्रक्षेपण मैट्रिक्स चाहते हैं, y प्रतिक्रिया वेक्टर है, जबकि H = X(X'X)^{-1}X' टोपी/प्रक्षेपण मैट्रिक्स है ताकि Hy भविष्य कहनेवाला मूल्यों देता है।
यह उत्तर साधारण कम से कम वर्गों के संदर्भ को मानता है। भारित कम से कम वर्गों के लिए, Get hat matrix from QR decomposition for weighted least square regression देखें।
अवलोकन
solve एक सामान्य वर्ग मैट्रिक्स की LU गुणन पर आधारित है। X'X के लिए (crossprod(X) द्वारा t(X) %*% X के बजाय आर में ?crossprod के लिए गणना की जानी चाहिए), जो सममित है, हम chol2inv का उपयोग कर सकते हैं जो Choleksy कारककरण पर आधारित है।
हालांकि, त्रिभुज कारक 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 और Cholesky कारक टूट जाएगा।
आम तौर पर, हम इस स्थिति में एसवीडी या क्यूआर पर स्विच करते हैं, लेकिन पिवोट चोलस्की कारककरण एक और विकल्प है।
- एसवीडी सबसे स्थिर विधि है, लेकिन बहुत महंगा है;
- क्यूआर मध्यम कम्प्यूटेशनल लागत पर संतोषजनक रूप से स्थिर है, और आमतौर पर अभ्यास में उपयोग किया जाता है;
- स्वीकार्य स्थिरता के साथ पिवोटेड चॉस्की तेजी से है। बड़े मैट्रिक्स के लिए इसे प्राथमिकता दी जाती है।
निम्नलिखित में, मैं सभी तीन विधियों को समझाऊंगा।
क्यूआर गुणन का उपयोग
ध्यान दें कि प्रक्षेपण मैट्रिक्स क्रमचय स्वतंत्र है, जैसे कि, यह कोई बात नहीं हम साथ या पिवट बिना क्यूआर गुणन प्रदर्शन या नहीं।
आर में, qr.default क्यूआर गुणन पिवट किए Linpack दिनचर्या DQRDC गैर पिवट किए क्यूआर गुणन के लिए, और LAPACK दिनचर्या DGEQP3 ब्लॉक के लिए कॉल कर सकते हैं। के एक खिलौना मैट्रिक्स दोनों विकल्प पैदा करते हैं और जांचते हैं:
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' एक आवरण समारोह को परिभाषित:
H1 <- f(qr_linpack)
H2 <- f(qr_lapack)
mean(abs(H1 - H2))
# [1] 9.530571e-17
विलक्षण मूल्य अपघटन का उपयोग करना:
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 एक ही प्रक्षेपण मैट्रिक्स दे देखेंगे
आर में, svd एकवचन मूल्य अपघटन की गणना करता है। फिर
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
, हम एक ही प्रक्षेपण मैट्रिक्स मिलती है: हम अभी भी ऊपर के उदाहरण मैट्रिक्स X का उपयोग करें।
पाइवट किया गया Cholesky गुणन
का उपयोग करते हुए प्रदर्शन के लिए, हम अभी भी ऊपर के उदाहरण 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
फिर, हमें एक ही प्रक्षेपण मैट्रिक्स मिलता है।
क्या कोई कारण है कि आप प्रिंसकंप या प्रकोप का उपयोग नहीं करना चाहते हैं? प्रमुख घटकों की गणना हाथ से करने की आवश्यकता नहीं है। –
डर है कि कोई सामान्य समाधान नहीं है, अगर यह स्थिति संख्या है तो आपको कोई समस्या है। –
हाय, मैं पीसीए को जितना अनुमान लगाता हूं उतना ही लागू करने की कोशिश नहीं कर रहा हूं। मुझे यह अजीब लगता है कि मैं इसे डमी उपकरणों के एक साधारण मैट्रिक्स के लिए भी काम नहीं कर सकता, जब ऐसा लगता है कि कॉललाइनर नहीं है। – bikeclub