MakePositiveDefinite {fake}

Making positive definite matrix

Description

Determines the diagonal entries of a symmetric matrix to make it is positive definite.

Usage

MakePositiveDefinite(
  omega,
  pd_strategy = "diagonally_dominant",
  ev_xx = NULL,
  scale = TRUE,
  u_list = c(1e-10, 1),
  tol = .Machine$double.eps^0.25
)

Arguments

omega

input matrix.

pd_strategy

method to ensure that the generated precision matrix is positive definite (and hence can be a covariance matrix). If pd_strategy="diagonally_dominant", the precision matrix is made diagonally dominant by setting the diagonal entries to the sum of absolute values on the corresponding row and a constant u. If pd_strategy="min_eigenvalue", diagonal entries are set to the sum of the absolute value of the smallest eigenvalue of the precision matrix with zeros on the diagonal and a constant u.

ev_xx

expected proportion of explained variance by the first Principal Component (PC1) of a Principal Component Analysis. This is the largest eigenvalue of the correlation (if scale_ev=TRUE) or covariance (if scale_ev=FALSE) matrix divided by the sum of eigenvalues. If ev_xx=NULL (the default), the constant u is chosen by maximising the contrast of the correlation matrix.

scale

logical indicating if the proportion of explained variance by PC1 should be computed from the correlation (scale=TRUE) or covariance (scale=FALSE) matrix.

u_list

vector with two numeric values defining the range of values to explore for constant u.

tol

accuracy for the search of parameter u as defined in optimise.

Details

Two strategies are implemented to ensure positive definiteness: by diagonally dominance or using eigendecomposition.

A diagonally dominant symmetric matrix with positive diagonal entries is positive definite. With pd_strategy="diagonally_dominant", the diagonal entries of the matrix are defined to be strictly higher than the sum of entries on the corresponding row in absolute value, which ensures diagonally dominance. Let \Omega* denote the input matrix with zeros on the diagonal and \Omega be the output positive definite matrix. We have:

\Omega_{ii} = \sum_{j = 1}^p | \Omega_{ij}* | + u, where u > 0 is a parameter.

A matrix is positive definite if all its eigenvalues are positive. With pd_strategy="diagonally_dominant", diagonal entries of the matrix are defined to be higher than the absolute value of the smallest eigenvalue of the same matrix with a diagonal of zeros. Let \lambda_1 denote the smallest eigenvvalue of the input matrix \Omega* with a diagonal of zeros, and v_1 be the corresponding eigenvector. Diagonal entries in the output matrix \Omega are defined as:

\Omega_{ii} = | \lambda_1 | + u, where u > 0 is a parameter.

It can be showed that \Omega has stricly positive eigenvalues. Let \lambda and v denote any eigenpair of \Omega*:

\Omega* v = \lambda v

\Omega* v + (| \lambda_1 | + u) v = \lambda v + (| \lambda_1 | + u) v

(\Omega* + (| \lambda_1 | + u) I) v = (\lambda + | \lambda_1 | + u) v

\Omega v = (\lambda + | \lambda_1 | + u) v

The eigenvalues of \Omega are equal to the eigenvalues of \Omega* plus | \lambda_1 |. The smallest eigenvalue of \Omega is (\lambda_1 + | \lambda_1 | + u) > 0.

Considering the matrix to make positive definite is a precision matrix, its standardised inverse matrix is the correlation matrix. In both cases, the magnitude of correlations is controlled by the constant u.

If ev_xx=NULL, the constant u is chosen to maximise the Contrast of the corresponding correlation matrix.

If ev_xx is provided, the constant u is chosen to generate a correlation matrix with required proportion of explained variance by the first Principal Component, if possible. This proportion of explained variance is equal to the largest eigenvalue of the correlation matrix divided by the sum of its eigenvalues. If scale=FALSE, the covariance matrix is used instead of the correlation matrix for faster computations.

Value

A list with:

omega

positive definite matrix.

u

value of the constant u.

References

Bodinier B, Filippi S, Nost TH, Chiquet J, Chadeau-Hyam M (2021). “Automated calibration for stability selection in penalised regression and graphical models: a multi-OMICs network application exploring the molecular response to tobacco smoking.” https://arxiv.org/abs/2106.02521.

Examples

# Simulation of a symmetric matrix
p <- 5
set.seed(1)
omega <- matrix(rnorm(p * p), ncol = p)
omega <- omega + t(omega)
diag(omega) <- 0

# Diagonal dominance maximising contrast
omega_pd <- MakePositiveDefinite(omega,
  pd_strategy = "diagonally_dominant"
)
eigen(omega_pd$omega)$values # positive eigenvalues

# Diagonal dominance with specific proportion of explained variance by PC1
omega_pd <- MakePositiveDefinite(omega,
  pd_strategy = "diagonally_dominant",
  ev_xx = 0.55
)
lambda_inv <- eigen(cov2cor(solve(omega_pd$omega)))$values
max(lambda_inv) / sum(lambda_inv) # expected ev

# Version not scaled (using eigenvalues from the covariance)
omega_pd <- MakePositiveDefinite(omega,
  pd_strategy = "diagonally_dominant",
  ev_xx = 0.55, scale = FALSE
)
lambda_inv <- 1 / eigen(omega_pd$omega)$values
max(lambda_inv) / sum(lambda_inv) # expected ev

# Non-negative eigenvalues maximising contrast
omega_pd <- MakePositiveDefinite(omega,
  pd_strategy = "min_eigenvalue"
)
eigen(omega_pd$omega)$values # positive eigenvalues

# Non-negative eigenvalues with specific proportion of explained variance by PC1
omega_pd <- MakePositiveDefinite(omega,
  pd_strategy = "min_eigenvalue",
  ev_xx = 0.7
)
lambda_inv <- eigen(cov2cor(solve(omega_pd$omega)))$values
max(lambda_inv) / sum(lambda_inv)

# Version not scaled (using eigenvalues from the covariance)
omega_pd <- MakePositiveDefinite(omega,
  pd_strategy = "min_eigenvalue",
  ev_xx = 0.7, scale = FALSE
)
lambda_inv <- 1 / eigen(omega_pd$omega)$values
max(lambda_inv) / sum(lambda_inv)