LDL: LDL' Decomposition of a Symmetric Positive-Definite Matrix

Description

Performs an LDL' factorization of a symmetric positive-definite matrix $X$, such that $$X = L D L^\prime,$$ where $L$ is unit lower-triangular (ones on the diagonal) and $D$ is diagonal.

Usage

LDL(x, epsilon = 1e-10)
InvLDL(s_l, uc_d)

Value

LDL(): a list with components:
- l: a unit lower-triangular matrix $L$
- s_l: a strictly lower-triangular part of $L$
- d: a vector of diagonal entries of $D$
- uc_d: unconstrained vector with $\mathrm{softplus}(uc\_d) = d$
- x: input matrix (with diagonal zeros possibly replaced by epsilon)
- epsilon: the epsilon value used
InvLDL(): a symmetric positive definite matrix

Arguments

x: Numeric matrix. Assumed symmetric positive-definite (not checked). Note: LDL() may error if the implied diagonal entries of $D$ are not strictly positive.
epsilon: Numeric. Small positive value used to replace exactly zero diagonal entries of x prior to factorization.
s_l: Matrix. Strictly lower-triangular part of $L$. In InvLDL(), only the strictly lower triangle is used (upper triangle and diagonal are ignored).
uc_d: Vector. Unconstrained vector such that Softplus(uc_d) = d, where d are the diagonal entries of $D$.

Details

LDL() returns both the unit lower-triangular factor $L$ and the diagonal factor $D$. The strictly lower-triangular part of $L$ is also provided for convenience. The function additionally computes an unconstrained vector uc_d such that Softplus(uc_d) = d. This uses a numerically stable inverse softplus implementation based on log(expm1(d)) (and a large-d rewrite), rather than the unstable expression $\log(\exp(d) - 1)$.

InvLDL() returns a symmetric positive definite matrix from the strictly lower-triangular part of $L$ and the unconstrained vector uc_d. The reconstructed matrix is symmetrized as $(\Sigma + \Sigma^\prime)/2$ to reduce numerical asymmetry.

Examples

Run this code

set.seed(123)
x <- crossprod(matrix(rnorm(16), 4, 4)) + diag(1e-6, 4)
ldl <- LDL(x = x)
ldl
inv_ldl <- InvLDL(s_l = ldl$s_l, uc_d = ldl$uc_d)
inv_ldl
max(abs(x - inv_ldl))