在R中反转三角矩阵的方法是什么？答案

【问题标题】：What is the way to invert a triangular matrix in R?在R中反转三角矩阵的方法是什么？
【发布时间】：2014-09-04 10:04:44
【问题描述】：

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

【问题讨论】：

最好提供一些示例输入和预期输出，以确保愿意提供可能答案的人回答正确的问题。
我认为这个问题已经很清楚了。我只是在寻找一个函数名称，或者找不到一个以有效方式计算三角矩阵逆的算法。
qr.solve 和 solve 在您似乎暗示的意义上并不等同。第一个使用 QR 分解，第二个使用 LU 分解。你为什么要逆？如果你想解决一个三角系统，看看backsolve和forwardsolve。
其实我已经通过QR分解得到了我的上三角矩阵R。我想要的只是高效地计算 R^-1。
这只是一个求解矩阵方程的例子，而不是矩阵求逆的例子。

标签： r matrix-inverse triangular

【解决方案1】：

尝试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

【讨论】：

【解决方案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)。

【讨论】：

【解决方案3】：

我最近也遇到了同样的问题。我的解决方案是加载 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 的倒数。不过，一个小问题可能是上述两个 R 函数中的每一个的结果都属于 Matrix 类。如果有人对此有疑问，请根据 r.h.s.是向量或矩阵（在代数意义上）。

【讨论】：

【解决方案4】：

似乎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]

【讨论】：

感谢您的测试，但我不确定
```
qrsolve
```
是否真的利用了输入矩阵的三角形形式。它必须尝试将其分解为 QR 形式。不过，有一些特定（和 FAST）算法旨在反转此类矩阵。

【解决方案5】：

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最慢，solve速度排名第三，backsolve速度排名第二，chol2inv最快。

cma

微基准（qr.solve（ma），solve（ma），cma

【讨论】：

chol2inv(M) 不返回 M 的倒数（即使 M 是三角形的）