blockChol {gek} | R Documentation |
Block Cholesky Decomposition
Description
blockChol
calculates the block Cholesky decomposition of a partitioned matrix of the form
A = \begin{pmatrix} K & R^{\rm T} \\ R & S \end{pmatrix}.
Usage
blockChol(K, R = NULL, S = NULL, tol = NULL)
Arguments
K |
a real symmetric positive-definite square submatrix. |
R |
an (optinal) rectangular submatrix. |
S |
an (optional) real symmetric positive-definite square submatrix. |
tol |
an (optional) numeric tolerance, see ‘Details’. |
Details
To obtain the block Cholesky decomposition
\begin{pmatrix} K & R^{\rm T} \\ R & S \end{pmatrix} = \begin{pmatrix} L^{\rm T} & 0 \\ Q^{\rm T} & M^{\rm T} \end{pmatrix} \begin{pmatrix} L & Q \\ 0 & M \end{pmatrix}
the following steps are performed:
Calculate
K = L^{\rm T}L
with upper triangular matrixL
.Solve
L^{\rm T}Q = R^{\rm T}
via forward substitution.Compute
N = S - Q^{\rm T}Q
the Schur complement of the blockK
of the matrixA
.Determine
N = M^{\rm T}M
with upper triangular matrixM
.
The upper triangular matrices L
and M
in step 1 and 4 are obtained by chol
.
Forward substitution in step 2 is carried out with backsolve
and the option transpose = TRUE
.
If tol
is specified a regularization of the form A_{\epsilon} = A + \epsilon I
is conducted.
Here, tol
is the upper bound for the logarithmic condition number of A_{\epsilon}
.
Then
\epsilon = \max\left\{ \frac{\lambda_{\max}(\kappa(A) - e^{\code{tol}})}{\kappa(A)(e^{\code{tol}} - 1)}, 0 \right\}
is chosen as the minimal "nugget" that is added to the diagonal of A
to ensure \log(\kappa(A_{\epsilon})) \leq
tol
.
Within gek this function is used to calculate the block Cholesky decomposition of the correlation matrix with derivatives.
Here K
is the Kriging correlation matrix.
R
is the matrix containing the first partial derivatives and S
consists of the second partial derivatives of the correlation matrix K
.
Value
blockChol
returns a list with the following components:
L |
the upper triangular factor of the Cholesky decomposition of |
Q |
the solution of the triangular system |
M |
the upper triangular factor of the Cholesky decomposition of the Schur complement |
If R
or S
are not specified, Q
and M
are set to NULL
,
i.e. only the Cholesky decomposition of K
is calculated.
The attribute "eps"
gives the minimum “nugget” that is added to the diagonal.
Warning
As in chol
there is no check for symmetry.
Author(s)
Carmen van Meegen
References
Chen, J., Jin, Z., Shi, Q., Qiu, J., and Liu, W. (2013). Block Algorithm and Its Implementation for Cholesky Factorization.
Gustavson, F. G. and Jonsson, I. (2000). Minimal-storage high-performance Cholesky factorization via blocking and recursion. IBM Journal of Research and Development, 44(6):823–850. doi:10.1147/rd.446.0823.
Ranjan, P., Haynes, R. and Karsten, R. (2011). A Computationally Stable Approach to Gaussian Process Interpolation of Deterministic Computer Simulation Data. Technometrics, 53:366–378. doi:10.1198/TECH.2011.09141.
See Also
chol
for the Cholesky decomposition.
backsolve
for backward substitution.
blockCor
for computing a correlation matrix with derivatives.
Examples
# Construct correlation matrix
x <- matrix(seq(0, 1, length = 5), ncol = 1)
res <- blockCor(x, theta = 1, covtype = "gaussian", derivatives = TRUE)
A <- cbind(rbind(res$K, res$R), rbind(t(res$R), res$S))
# Cholesky decomposition of correlation matix without derivatives
cholK <- blockChol(res$K)
cholK
cholK$L == chol(res$K)
# Cholesky decomposition of correlation matix with derivatives
cholA <- blockChol(res$K, res$R, res$S)
cholA <- cbind(rbind(cholA$L, matrix(0, ncol(cholA$Q), nrow(cholA$Q))),
rbind(cholA$Q, cholA$M))
cholA
cholA == chol(A)
# Cholesky decomposition of correlation matix with derivatives with regularization
res <- blockCor(x, theta = 2, covtype = "gaussian", derivatives = TRUE)
A <- cbind(rbind(res$K, res$R), rbind(t(res$R), res$S))
try(blockChol(res$K, res$R, res$S))
blockChol(res$K, res$R, res$S, tol = 35)