CODEWITHSUNDEEP
rmatrix-inversetriangular
Mark Morrisson
Mark Morrisson

Reputation: 2703

What is the way to invert a triangular matrix in R?

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?

Upvotes: 3

Views: 4025

Answers (5)

NdNg
NdNg

Reputation: 11

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).

Upvotes: 1

J.G.
J.G.

Reputation: 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))

Upvotes: -1

ClementWalter
ClementWalter

Reputation: 5314

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.

Upvotes: 2

NoBackingDown
NoBackingDown

Reputation: 2194

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

Upvotes: 4

Will Beason
Will Beason

Reputation: 3561

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]

Upvotes: 0

Related Questions