# solve

##### 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
`chol`

esky decomposition and clearly faster than
`solve`

.
It also returns the logarithm of the
determinant at no additional computational costs.

- Keywords
- math

##### 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

##### See Also

chol.spam in the package spam

##### Examples

```
# NOT RUN {
<!-- % library(RandomFieldsUtils) -->
# }
# NOT RUN {
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))
# }
```

*Documentation reproduced from package RandomFieldsUtils, version 0.5.3, License: GPL (>= 3)*