scipy.linalg.cho_solve counterpart in R?

I was wonderin开发者_JAVA百科g if there is a counterpart to scipy.linalg.cho_solve in R. What the function does is given the cholesky factor L of A (A = LL') and b, it solves the original problem, Ax = b. (not Lx = b)

(So it is different from backsolve/forwardsolve)

Thank you, Joon

I can't think of a function doing that for you automatically, but given you have the cholesky factor L, it's easily done in one line by reconstructing the A matrix as defined by the decomposition A=LL' :

 A=matrix(c(1,1,1,1,5,5,1,5,14),nrow=3)
 # Cholesky decomposition A = LL'
 L <- chol(A)

 # Make some b with known x
 x <- c(1,2,3)
 b <- A%*%x

 # Solve
 solve( t(L) %*% L, b)

edit: be aware that in R, the definition of the Cholesky factor is related to A = L'L, which is why you have to put the transposed first in the solve.

edit2 : After reading Bates article, I realized it should be:

> solve(crossprod(L),b)
     [,1]
[1,]    1
[2,]    2
[3,]    3

If I understand you correctly, then Doug Bates covered some of this in an article he wrote for R News in 2004 (see page 18 of the linK).

The relevant bit is:

ch <- chol(crossprod(X))
chol.sol <- backsolve(ch, forwardsolve(ch, crossprod(X, y),
                                       upper = TRUE, trans = TRUE))

where X is the matrix of predictor variables.

Doug's article goes on to show how functionality in the Matrix package (which comes with R) can be used solve the same system very quickly indeed.

I realise this question is a little old, but I see that the answer

forwardsolve(L, forwardsolve(L, b), transp=TRUE)

hasn't been given yet. This uses the triangular structure, while keeping to the original question. This should be faster and more accurate for larger matrices. It might also be worth noting that L <- t(chol(A)) since chol returns an upper triangular matrix.

A <- matrix(c(1,1,1,1,5,5,1,5,14), nrow=3)
# Cholesky decomposition A = LL'
L <- t(chol(A))

# Make some b with known x
x <- c(1, 2, 3)
b <- A %*% x

# Solve
forwardsolve(L, forwardsolve(L, b), transp=TRUE)

Giving the answer:

> forwardsolve(L, forwardsolve(L, b), transp=TRUE)
     [,1]
[1,]    1
[2,]    2
[3,]    3