solvePosDef: Solve a System of Equations for Positive Definite Matrices

In RandomFieldsUtils: Utilities for the Simulation and Analysis of Random Fields and Genetic Data

Solve a System of Equations for Positive Definite Matrices

Description

This function solves the equality a x = b for x where a is a positive definite matrix and b is a vector or a matrix. It is slightly faster than the inversion by the cholesky decomposition and clearly faster than solve. It also returns the logarithm of the determinant at no additional computational costs.

Usage

solvex(a, b=NULL, logdeterminant=FALSE)

Arguments

a	a square real-valued matrix containing the coefficients of the linear system. Logical matrices are coerced to numeric.
b	a numeric or complex vector or matrix giving the right-hand side(s) of the linear system. If missing, b is taken to be an identity matrix and solvex will return the inverse of a.
logdeterminant	logical. whether the logarithm of the determinant should also be returned

Details

If the matrix is diagonal direct calculations are performed.

Else if the matrix is sparse the package spam is used.

Else the Cholesky decomposition is tried. Note that with RFoptions(pivot= ) pivoting can be enabled. Pivoting is about 30% slower.

If it fails, the eigen value decomposition is tried.

Value

If logdeterminant=FALSE the function returns a vector or a matrix, depending on b which is the solution to the linear equation. Else the function returns a list containing both the solution to the linear equation and the logarithm of the determinant of a.

References

See chol.spam of the package spam.

Examples

RFoptions(solve_method = "cholesky", printlevel=1)
set.seed(1)
n <- 1000
x <- 1:n
y <- runif(n)


## FIRST EXAMPLE: full rank matrix
M <- exp(-as.matrix(dist(x) / n)) 
b0 <- matrix(nr=n, runif(n * 5))
b <- M %*% b0 + runif(n)

## standard with 'solve'
print(system.time(z <- zR <- solve(M, b)))
print(range(b - M %*% z))
stopifnot(all(abs((b - M %*% z)) < 2e-11))

## using exactly the algorithm used in R
RFoptions(pivot=PIVOT_NONE, la_mode=LA_R) ## (default)
print(system.time(z <- solvex(M, b)))
print(range(b - M %*% z))
stopifnot(all(z == zR))

## Without pivoting, internal code:
RFoptions(pivot=PIVOT_NONE, la_mode=LA_INTERN) ## (default)
print(system.time(z <- solvex(M, b)))
print(range(b - M %*% z))
stopifnot(all(abs((b - M %*% z)) < 2e-11))

## Pivoting is slower here:
RFoptions(pivot=PIVOT_DO, la_mode=LA_INTERN) 
print(system.time(z <- solvex(M, b)))
print(range(b - M %*% z))
stopifnot(all(abs((b - M %*% z)) < 2e-11))


## SECOND EXAMPLE: low rank matrix
M <- x %*% t(x) + rev(x) %*% t(rev(x)) + y %*% t(y)
b1 <- M %*% b0

## Without pivoting, it does not work
RFoptions(pivot=PIVOT_NONE, la_mode=LA_R)
## Not run: try(solve(M, b1))
RFoptions(pivot=PIVOT_NONE, la_mode=LA_INTERN)
## Not run: try(solvex(M, b1))

## Pivoting works -- the precision however is reduced :
RFoptions(pivot=PIVOT_DO, la_mode=LA_INTERN) 
print(system.time(z1 <- solvex(M, b1)))
print(range(b1 - M %*% z1))
stopifnot(all(abs((b1 - M %*% z1)) < 2e-6))

## Pivoting fails, when the equation system is not solvable:
b2 <- M + runif(n)
## Not run: try(solvex(M, b2))


RFoptions(pivot = PIVOT_AUTO, la_mode = LA_AUTO)

RandomFieldsUtils documentation built on April 19, 2022, 5:09 p.m.