Solve a System of Equations for Positive Definite Matrices
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.
solvex(a, b=NULL, logdeterminant=FALSE)
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, |
logdeterminant |
logical. whether the logarithm of the determinant should also be returned |
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.
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.
Martin Schlather, schlather@math.uni-mannheim.de, http://ms.math.uni-mannheim.de
See chol.spam of the package spam
chol.spam in the package spam
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 <- solve(M, b))) print(range(b - M %*% z)) stopifnot(all(abs((b - M %*% z)) < 2e-11)) ## Without pivoting: RFoptions(pivot=PIVOT_NONE) ## (default) print(system.time(z <- solvex(M, b))) print(range(b - M %*% z)) stopifnot(all(abs((b - M %*% z)) < 2e-11)) ## Pivoting is 35% slower here: RFoptions(pivot=PIVOT_DO) 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) #try(solve(M, b1)) #try(solvex(M, b1)) ## Pivoting works -- the precision however is reduced : RFoptions(pivot=PIVOT_DO) 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) #try(solvex(M, b2))
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.