3

我有一个上三角矩阵，我想快速计算它的逆矩阵。我试过qr.solve()了，但我觉得它等同于solve()，并且它没有利用输入矩阵的三角形性质。最好的方法是什么？

4

5 回答 5

3

尝试backsolve()使用具有适当维度的单位矩阵作为右手值。

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
于 2015-08-02T12:54:00.257 回答
2

如果你想得到你的矩阵的逆M，你可以简单地使用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

运行此代码时M%*%Minv，由于数值近似，您可能会看到系数非常低 (~10^-15)。

于 2016-04-13T08:53:44.027 回答
1

我最近也遇到了同样的问题。我的解决方案是加载 Matrix 包，然后定义以下两个 R 函数作为该任务的包装器：

my.backsolve <- function(A, ...) solve(triu(A), ...)  
my.forwardsolve <- function(A, ...) solve(tril(A), ...)

这些需要 Matrix 包，因为triu并tril在其中定义。然后my.backsolve(A)(resp. my.forwardsolve(A)) 计算A上 (resp. 下) 三角形情况的倒数。不过，一个小问题可能是上述两个 R 函数中的每一个的结果都是 class Matrix。如果有人对此有疑问，则根据 rhs 是向量还是矩阵（在代数意义上）应用as.vector或应用于结果。as.matrix

于 2019-07-10T09:18:00.213 回答
0

它似乎solve比执行速度更快，qr.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

请注意，如果我们+3在编辑矩阵单元格时删除 ，solve几乎总是会失败并且qr.solve通常仍然会给出答案。

> 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]
于 2014-09-04T16:06:56.007 回答
-1

R 基函数chol2inv可能正在做反转三角矩阵的技巧。从它的 反转一个对称的正定方阵Choleski decomposition。等效地，(X'X)^(-1)从 的 (R 部分)计算QR decomposition of X。更一般地，给定一个上三角矩阵 R，计算(R'R)^(-1)。见https://stat.ethz.ch/R-manual/R-devel/library/base/html/chol2inv.html

我用微基准法测试了它的速度：qr.solve最慢，solve速度排名第三，backsolve速度排名第二，chol2inv最快。

cma <- chol(ma <- cbind(1, 1:3, c(1,3,7)))

微基准（qr.solve（ma），solve（ma），cma<-chol（ma），chol2inv（cma），invcma<-backsolve（cma，diag（nrow（cma））），backinvma<-tcrossprod（invcma） )

于 2016-12-15T21:33:15.517 回答