Open shz9 opened 2 years ago
Another option that seems to work well in some settings is to use the shrinkage algorithm of Higham et al. (2014). Here's a rough implementation of the idea:
from scipy.sparse import identity
from scipy.sparse.linalg import eigsh, ArpackNoConvergence
def isPSD_sparse(B):
"""
Check for positive semi-definitness using the ARPACK library by
examining the smallest eigenvalues.
"""
try:
return eigsh(B, 1, which='SA', return_eigenvectors=False)[0] >= 0.
except ArpackNoConvergence:
return eigsh(B, 1, which='SA', sigma=0, return_eigenvectors=False)[0] >= 0.
def nearcorr_bisection(m0, tol=1e-5):
"""
Find the nearest positive semi-definite correlation matrix
using the bisection method (Algorithm 3.1) of Higham et al. (2014).
The method finds an intermediate matrix between the original `m0`
and the identity matrix of the same shape. The benefit of this is that
it retains the sparsity structure in the original matrix, which is
important for large covariance matrices in genomics.
One major difference from Higham et al. is that we use the scipy's
and ARPACK's sparse matrix operations to find the smallest eigen
values efficiently, instead of performing cholesky decomposition to check for
positive semi-definitness.
The method specifically returns the `alpha` value from the shrinkage
procedure, instead of the updated matrix. This is because under this
formulation, we'd only need to multiply the off-diagonal elements by
(1 - a_*) to obtain the shrunk correlation matrix.
"""
a_l = 0.
a_r = 1.
m1 = identity(m0.shape[0])
while a_r - a_l > tol:
print(a_l, a_r)
a_m = .5*(a_l + a_r)
if not isPSD_sparse(a_m*m1 + (1. - a_m)*m0):
a_l = a_m
else:
a_r = a_m
return a_r
Needs more testing, but looks like a promising approach that can be very fast.
In some applications, we may require that the LD matrix (SNP-by-SNP correlation matrix) be positive definite (PD) or positive semi-definite (PSD). Very often, this is not the case, which may result in instabilities for some downstream models. Therefore, it would be great to augment our LD functionalities to check for PSD and to also find the nearest PSD matrix to the sample LD matrix.
A good starting point may be the code snippet in this StackOverflow thread. I modified this function to find the nearest PSD while retaining the sparsity structure in the original matrix:
The main bottleneck here is that we need to perform SVD on the LD matrix, which may be computationally expensive, even with sparse arrays. One way to get around this difficulty is to deploy this within LD blocks, which may be more manageable.
To experiment with this, we can try to tackle the following tasks:
LDMatrix
corresponding to the LD blocks defined by the block estimator.viprs
).