R - Standardize matrix to have unit diagonals

Admin 1
Viewed 77

I am seeking to generate the below matrix:

Θ = B + δIp ∈ Rp×p, where Ip is the identity matrix, each off-diagonal entry in B (symmetric matrix) is generated independently and equals 0.5 with probability 0.1 or 0 with probability 0.9. The parameter δ > 0 is chosen such that Θ is positive definite. The matrix is standardized to have unit diagonals (transforming from covariance matrix to correlation matrix).

I think that I have most of the code, but I'm unsure of how to standardize the matrix to have unit diagonals syntactically in R (and theoretically, why that is a useful feature of a matrix).

# set number of cols/rows
p <- 5
set.seed(123)
# generate matrix B with values of 0.5 given probabilities
B <- matrix(sample(c(0,0.5), p^2, replace=TRUE, prob=c(0.9,0.1)), p)
# call the matrix lower triangle, need a symmetric matrix
i <- lower.tri(B)
B[i] <- t(B)[i]
diag(B) <- rep(0, p)
# finding parameter delta, such that Θ is positive definite.
(delta <- -min(eigen(B, symmetric=TRUE, only.values=TRUE)$values))
# set theta (delta is 2.8802)
theta <- B + 2.89*(diag(p))
# now to standardize the matrix to have unit diagonals ?
r matrix
1 Answers

There are many ways to do this, but the following is very fast in timing experiments:

v <- 1/sqrt(diag(theta))
B <- theta * outer(v, v)

This divides all rows and columns by their standard deviations, which are the square roots of the diagonal elements.

It will fail whenever any diagonal is zero (or negative): but in that case such standardization isn't possible. Computing the square roots and their reciprocals first allows you to learn as soon as possible--with minimal computation--whether the procedure will succeed.

BTW, a direct way to compute B in the first part of your code (which has a zero diagonal) is

B <- as.matrix(structure(sample(c(0,1/2), p*(p-1)/2, replace=TRUE, prob=c(.9,.1),
               Size=p, Diag=TRUE, class="dist"))

This eliminates the superfluous sampling.

user_8447230
Related