Inverse from Cholesky Factor

Given formally upper and lower triangular matrices \(U\) and \(L\), compute \((U' U)^{-1}\) and \((L L')^{-1}\), respectively.

This function can be seen as way to compute the inverse of a symmetric positive definite matrix given its Cholesky factor. Equivalently, it can be seen as a way to compute \((X' X)^{-1}\) given the \(R\) part of the QR factorization of \(X\), if \(R\) is constrained to have positive diagonal entries.

chol2inv(x, ...)
# S4 method for class 'dtrMatrix'
chol2inv(x, ...)
# S4 method for class 'dtCMatrix'
chol2inv(x, ...)
# S4 method for class 'generalMatrix'
chol2inv(x, uplo = "U", ...)

Arguments

x

a square matrix or Matrix, typically the result of a call to chol. If x is square but not (formally) triangular, then only the upper or lower triangle is considered, depending on optional argument uplo if x is a Matrix.

uplo

a string, either "U" or "L", indicating which triangle of x contains the Cholesky factor. The default is "U", to be consistent with chol2inv from base.

...

further arguments passed to or from methods.

Value

A matrix, symmetricMatrix, or diagonalMatrix representing the inverse of the positive definite matrix whose Cholesky factor is x. The result is a traditional matrix if x is a traditional matrix, dense if x is dense, and sparse if x is sparse.

See also

The default method from base, chol2inv, called for traditional matrices x.

Generic function chol, for computing the upper triangular Cholesky factor \(L'\) of a symmetric positive semidefinite matrix.

Generic function solve, for solving linear systems and (as a corollary) for computing inverses more generally.

Examples

(A <- Matrix(cbind(c(1, 1, 1), c(1, 2, 4), c(1, 4, 16))))
#> 3 x 3 Matrix of class "dsyMatrix"
#>      [,1] [,2] [,3]
#> [1,]    1    1    1
#> [2,]    1    2    4
#> [3,]    1    4   16
(R <- chol(A))
#> 3 x 3 Matrix of class "dtrMatrix"
#>      [,1]    [,2]    [,3]   
#> [1,] 1.00000 1.00000 1.00000
#> [2,]       . 1.00000 3.00000
#> [3,]       .       . 2.44949
(L <- t(R))
#> 3 x 3 Matrix of class "dtrMatrix"
#>      [,1]    [,2]    [,3]   
#> [1,] 1.00000       .       .
#> [2,] 1.00000 1.00000       .
#> [3,] 1.00000 3.00000 2.44949
(R2i <- chol2inv(R))
#> 3 x 3 Matrix of class "dpoMatrix"
#>            [,1] [,2]       [,3]
#> [1,]  2.6666667 -2.0  0.3333333
#> [2,] -2.0000000  2.5 -0.5000000
#> [3,]  0.3333333 -0.5  0.1666667
(L2i <- chol2inv(R))
#> 3 x 3 Matrix of class "dpoMatrix"
#>            [,1] [,2]       [,3]
#> [1,]  2.6666667 -2.0  0.3333333
#> [2,] -2.0000000  2.5 -0.5000000
#> [3,]  0.3333333 -0.5  0.1666667
stopifnot(exprs = {
    all.equal(R2i, tcrossprod(solve(R)))
    all.equal(L2i,  crossprod(solve(L)))
    all.equal(as(R2i %*% A, "matrix"), diag(3L)) # the identity 
    all.equal(as(L2i %*% A, "matrix"), diag(3L)) # ditto
})