gek (version 1.2.0)

blockChol: 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)

Value

blockChol returns a list with the following components:

L

the upper triangular factor of the Cholesky decomposition of K.

Q

the solution of the triangular system t(L) %*% Q == t(R).

M

the upper triangular factor of the Cholesky decomposition of the Schur complement N.

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.

Arguments

K

a real symmetric positive-definite square submatrix.

R

an (optional) rectangular submatrix.

S

an (optional) real symmetric positive-definite square submatrix.

tol

an (optional) numeric tolerance, see ‘Details’.

Author

Carmen van Meegen

Warning

As in chol there is no check for symmetry.

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:

1. Calculate \(K = L^{\rm T}L\) with upper triangular matrix \(L\).

2. Solve \(L^{\rm T}Q = R^{\rm T}\) via forward substitution.

3. Compute \(N = S - Q^{\rm T}Q\) the Schur complement of the block \(K\) of the matrix \(A\).

4. Determine \(N = M^{\rm T}M\) with upper triangular matrix \(M\).

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.

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. tools:::Rd_expr_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. tools:::Rd_expr_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)