R에서 삼각 행렬을 반전시키는 방법은 무엇입니까?

Question

필자는 위 삼각 행렬을 가지고 있으며 그 역수를 빠르게 계산하려고합니다. 나는 qr.solve()을 시도했지만 solve()과 같으며 입력 행렬의 삼각형 특성을 이용하지 않는다고 생각합니다. 가장 좋은 방법은 무엇입니까?

Answer 1

solve은 qr.solve보다 다소 빠른 성능을 나타내지 만 qr.solve은 더욱 견고합니다. 매트릭스 세포를 편집 할 때 우리는 +3를 제거하면, solve 거의 항상 실패하고 qr.solve 여전히 일반적으로 대답을 줄 것이다

는

n <- 50 
mymatrix <- matrix(0, nrow=n, ncol=n) 

fun1 <- function() { 
    for (i in 1:n) { 
    mymatrix[i, i:n] <- rnorm(n-i+1)+3 
    } 
    solve(mymatrix) 
} 

fun2 <- function() { 
    for (i in 1:n) { 
    mymatrix[i, i:n] <- rnorm(n-i+1)+3 
    } 
    qr.solve(mymatrix) 
} 

> system.time(for (i in 1:1000) fun1()) 
    user system elapsed 
    1.92 0.03 1.95 
> system.time(for (i in 1:1000) fun2()) 
    user system elapsed 
    2.92 0.00 2.92 

참고.

> set.seed(0) 
> fun1() 
Error in solve.default(mymatrix) : 
    system is computationally singular: reciprocal condition number = 3.3223e-22 
> set.seed(0) 
> fun2() 
[returns the inverted matrix] 

Answer 2

backsolve()을 시도하고 적절한 크기의 단위 행렬을 오른쪽 값으로 사용하십시오. 당신이 당신의 매트릭스 M의 역을 얻고 싶다면

library(microbenchmark) 
n <- 2000 
n.cov <- 1000 
X <- matrix(sample(c(0L, 1L), size = n * n.cov, replace = TRUE), 
      nrow = n, ncol = n.cov) 
R <- chol(crossprod(X)) 
rm(X) 
microbenchmark(
    backsolve = backsolve(r = R, x = diag(ncol(R))), 
    solve = solve(R), 
times = 10) 

Unit: milliseconds 
     expr  min  lq  mean median  uq  max neval 
backsolve 467.2802 467.5142 469.4457 468.1578 468.6501 482.2431 10 
    solve 748.2351 748.8318 750.0764 750.3319 750.9583 751.5005 10 

Answer 3

, 당신은 단순히 backsolve(M, diag(dim(M)[1]))를 사용할 수 있습니다. 예를 들어 :

M <- matrix(rnorm(100), 10, 10) # random matrix 
M[lower.tri(M)] <- 0 # make it triangular 
Minv <- backsolve(M, diag(dim(M)[1])) 
M%*%Minv 

이 코드를 실행하면 인해 수치 근사 매우 낮은 계수 (~ 10^-15)으로 M%*%Minv를 볼 수 있습니다.

Answer 4

R베이스 함수 chol2inv은 삼각 행렬을 반전시키기위한 트릭을 수행 중일 수 있습니다. 대칭의 양의 정사각형 정사각형 행렬을 Choleski decomposition에서 반전합니다. 동등하게, QR decomposition of X의 (R 부분)에서 (X'X)^(-1)을 계산하십시오. 더욱 일반적으로, 위 삼각 행렬 R이 주어진다면, (R'R)^(-1)을 계산하십시오. 나는 마이크로 벤치의 방법으로 자사의 속도를 테스트 https://stat.ethz.ch/R-manual/R-devel/library/base/html/chol2inv.html

참조 : qr.solve가 느린이다, 속도, 등급 3를 해결 backsolve는 chol2inv가 가장 빠른 속도에서 2 위를 기록하고 있습니다.

CMA < - CHOL (MA < - cbind (1, 1 : 3, C (1,3,7)))

마이크로 벤치 (qr.solve (MA), 해결 (MA) 당신은 몇 가지 샘플 입력을 제공하는 경우 CMA < -chol (MA), chol2inv (CMA), invcma < -backsolve (CMA, DIAG (nrow (CMA))), backinvma < -tcrossprod (invcma는))