3

I have an upper triangular matrix and I'd like to compute its inverse in a fast way. I tried qr.solve() but I have the feeling that it's equivalent to solve(), and that it does not exploit the triangular nature of the input matrix. What is the best way to do it?

Mark Morrisson
  • 2,543
  • 4
  • 19
  • 25
  • 1
    It's best if you provide some sample input and expected output to make sure that those who are willing to provide possible answers are answering the right question. – A5C1D2H2I1M1N2O1R2T1 Sep 04 '14 at 10:06
  • 1
    I think the question is clear enough. I'm only looking for a function name, or failing that, an algorithm to compute the inverse of a triangular matrix in an efficient way. – Mark Morrisson Sep 04 '14 at 10:27
  • `qr.solve` and `solve` are not equivalent in the sense you seem to imply. The first uses a QR decomposition and the second uses an LU decomposition. Why do you want an inverse? If you want to solve a triangular system, have a look at `backsolve` and `forwardsolve`. – Bhas Sep 04 '14 at 10:40
  • In fact, I have already obtained my upper triangular matrix R with a QR decomposition. All I want is to compute efficiently R^-1. – Mark Morrisson Sep 04 '14 at 10:44
  • This is only an example of solving a matrix equation, not of a matrix inversion. – Mark Morrisson Sep 04 '14 at 11:29

5 Answers5

3

Try backsolve() and use the identity matrix with appropriate dimension as the right-hand value.

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
NoBackingDown
  • 2,144
  • 1
  • 19
  • 25
2

If you want to get the inverse of your matrix M, you can simply use backsolve(M, diag(dim(M)[1])). For example:

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

When running this code you may see M%*%Minv with very low coefficient (~10^-15) due to numerical approximation.

ClementWalter
  • 4,814
  • 1
  • 32
  • 54
1

I have, of late, confronted the same problem. My solution has been to load the Matrix package, and then define the following two R functions serving as wrappers for that task:

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

These need the Matrix package because triuand tril are defined in it. Then my.backsolve(A) (resp. my.forwardsolve(A)) computes the inverse of A for the upper (resp. lower) triangular case. A minor issue, though, might be that the result of each of the above two R functions is of class Matrix. If one has issues with that, then apply as.vector or as.matrixto the result according to whether the r.h.s. is a vector or a matrix (in the algebraic sense).

NdNg
  • 11
  • 1
0

It seems solve performs somewhat faster than qr.solve but qr.solve is more robust.

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

Note that if we remove the +3 when editing the matrix cells, solve will almost always fail and qr.solve will still usually give an answer.

> 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]
Will Beason
  • 3,417
  • 2
  • 28
  • 46
  • 1
    Thanks for the tests, but I'm not sure that 
    qrsolve
    actually exploits the triangular form of the input matrix. It must try to decompose it into a QR form. Though, there are specific (and FAST) algorithms designed to invert such matrices.     – Mark Morrisson Sep 04 '14 at 18:35
-1

The R base function chol2inv might be doing the trick for inverting triangular matrix. Inverts a symmetric, positive definite square matrix from its Choleski decomposition. Equivalently, compute (X'X)^(-1) from the (R part) of the QR decomposition of X. Even more generally, given an upper triangular matrix R, compute (R'R)^(-1). See https://stat.ethz.ch/R-manual/R-devel/library/base/html/chol2inv.html

I tested its speed by microbenchmark method: qr.solve is slowest, solve ranks 3rd in speed, backsolve ranks 2nd in speed, chol2inv is the fastest.

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

microbenchmark(qr.solve(ma), solve(ma), cma<-chol(ma), chol2inv(cma), invcma<-backsolve(cma, diag(nrow(cma))), backinvma<-tcrossprod(invcma))

Community
  • 1
  • 1
J.G.
  • 1
  • 1