1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
\name{solve}
\alias{solvePosDef}
\alias{solvex}
\alias{solve}
\title{Solve a System of Equations for Positive Definite Matrices}
\description{
This function solves the equality \eqn{a x = b} for \eqn{x}
where \eqn{a} is a \bold{positive definite}
matrix and \eqn{b} is a vector or a
matrix. It is slightly faster than the inversion by the
\code{\link[base]{chol}}esky decomposition and clearly faster than
\code{\link[base]{solve}}.
It also returns the logarithm of the
determinant at no additional computational costs.
}
\usage{
solvex(a, b=NULL, logdeterminant=FALSE)
%, sparse=NA, method=-1)
}
\arguments{
\item{a}{a square real-valued matrix containing the
coefficients of the linear system. Logical matrices are
coerced to numeric.
}
\item{b}{
a numeric or complex vector or matrix giving the right-hand
side(s) of the linear system. If missing, \code{b} is taken to be
an identity matrix and \code{solvex} will return the inverse of
\code{a}.
}
\item{logdeterminant}{logical.
whether the logarithm of the determinant should also be returned
}
}
\value{
If \code{logdeterminant=FALSE} the function returns a vector or a matrix,
depending on \code{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 \code{a}.
}
\details{
% The values of \code{method} could be:
% \itemize{
% \item \code{<0} :
If the matrix is diagonal direct calculations are performed.
Else if the matrix is sparse the package \pkg{spam} is used.
Else the Cholesky decomposition is tried. Note that with
\code{RFoptions(pivot= )} pivoting can be enabled. Pivoting is about
30\% slower.
If it fails, the eigen value decomposition is tried.
}
\references{
% See \link[spam]{chol.spam} of the package \pkg{spam}
See \code{chol.spam} of the package \pkg{spam}.
}
%\seealso{
% \link[spam]{chol.spam} in the package \pkg{spam}
%}
\me
\keyword{math}
\examples{
% library(RandomFieldsUtils)
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)
\dontrun{try(solve(M, b1))}
RFoptions(pivot=PIVOT_NONE, la_mode=LA_INTERN)
\dontrun{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)
\dontrun{try(solvex(M, b2))}
RFoptions(pivot = PIVOT_AUTO, la_mode = LA_AUTO)
}
|
