Tính ma trận chiếu/hat qua QR thừa, SVD (và Cholesky thừa?)

Question

Tôi đang cố gắng để tính toán trong R một ma trận chiếu P của một ma trận J N x tùy ý S:

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

Tôi đã đã cố gắng để thực hiện điều này với các chức năng sau:

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

Nhưng khi tôi sử dụng này tôi nhận được lỗi mà trông như thế này:

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

Tôi nghĩ rằng đây là kết quả của sự cố và/hoặc sự bất ổn định số như được thảo luận ở nhiều nơi như r-help và here, nhưng tôi không đủ kinh nghiệm bằng cách sử dụng phân tích SVD hoặc QR để khắc phục sự cố hoặc đặt mã hiện tại này đi vao thực hiện. Tôi cũng đã thử mã được đề xuất, để viết giải quyết dưới dạng hệ thống:

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

Nhưng vẫn không hoạt động. Mọi lơi đê nghị đêu nên được đanh gia cao.

Tôi khá chắc chắn rằng ma trận của tôi sẽ không thể đảo ngược và không có bất kỳ tuyến tính nào, nếu chỉ vì tôi đã thử nghiệm điều này với ma trận của các biến giả trực giao và nó vẫn không hoạt động.

Ngoài ra, tôi muốn áp dụng điều này cho các ma trận khá lớn, vì vậy tôi đang tìm một giải pháp tổng quát gọn gàng.

Answer 1

Cách áp dụng chol đến S'S, sau đó chol2inv. Nên ổn định hơn:

chol2inv(chol(crossprod(S))) 

Answer 2

Mặc dù OP chưa hoạt động trong hơn một năm, tôi vẫn quyết định đăng câu trả lời. Tôi sẽ sử dụng X thay vì S, như trong thống kê, chúng tôi thường muốn ma trận chiếu trong bối cảnh hồi quy tuyến tính, trong đó X là ma trận mô hình, y là vector phản hồi, trong khi H = X(X'X)^{-1}X' là ma trận mũ/chiếu để Hy đưa ra giá trị tiên đoán.

Câu trả lời này giả định bối cảnh của các ô vuông nhỏ nhất thông thường. Đối với ô vuông có trọng số tối thiểu, hãy xem Get hat matrix from QR decomposition for weighted least square regression.

---

Tổng quan

solve được dựa trên LU thừa số của một ma trận vuông nói chung. Đối với X'X (nên được tính bằng crossprod(X) thay vì t(X) %*% X trong R, đọc ?crossprod để biết thêm) là đối xứng, chúng tôi có thể sử dụng chol2inv dựa trên hệ số Choleksy.

Tuy nhiên, hệ số hình tam giác kém ổn định hơn hệ số QR. Đây không phải là khó hiểu. Nếu X có số điều kiện kappa, X'X sẽ có số có điều kiện kappa^2. Điều này có thể gây ra khó khăn số lớn. Thông báo lỗi bạn nhận được:

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

chỉ nói điều này.kappa^2 là khoảng e-28, nhỏ hơn nhiều so với độ chính xác của máy tại khoảng e-16. Với dung sai tol = .Machine$double.eps, X'X sẽ bị coi là thiếu cấp, do đó hệ số LU và Cholesky sẽ bị hỏng.

Nói chung, chúng tôi chuyển sang SVD hoặc QR trong tình huống này, nhưng xoay vòng Hệ số yếu tố Cholesky là một lựa chọn khác.

• SVD là phương pháp ổn định nhất nhưng quá đắt;
• QR là đáp ứng ổn định, với chi phí tính toán vừa phải và thường được sử dụng trong thực tế;
• Xoay vòng Cholesky nhanh, có độ ổn định chấp nhận được. Đối với ma trận lớn, cái này được ưu tiên.

Sau đây, tôi sẽ giải thích cả ba phương pháp.

---

Sử dụng QR thừa

Lưu ý rằng ma trận chiếu là hoán vị độc lập, nghĩa là, nó không quan trọng cho dù chúng ta thực hiện QR thừa số có hoặc không xoay vòng.

Trong R, qr.default có thể gọi thường trình LINPACK DQRDC cho hệ số QR không được xoay vòng, và LAPACK thường quy DGEQP3 cho hệ số QR của khối được xoay vòng. Hãy tạo ma trận đồ chơi và kiểm tra cả hai tùy chọn:

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" 

Lưu ý rằng $pivot là khác nhau cho hai đối tượng.

Bây giờ, chúng ta định nghĩa một hàm wrapper để tính QQ':

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

Chúng ta sẽ thấy rằng qr_linpack và qr_lapack cung cấp cho các ma trận chiếu giống nhau:

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

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

---

Sử dụng phân hủy giá trị đặc biệt

Trong R, svd tính toán phân tách giá trị ít nhất. Chúng tôi vẫn sử dụng ma trận ví dụ trên 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 

Một lần nữa, chúng tôi có cùng ma trận chiếu.

---

Sử dụng xoay Cholesky thừa

Đối với cuộc biểu tình, chúng tôi vẫn sử dụng ví dụ X trên.

## 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 

Một lần nữa, chúng tôi có cùng ma trận chiếu.