Denoise a Covariance Matrix using Constant Residual Eigen Value Method
Practical use:
The Eigen Values are naturally stored in decreasing order. Once we find the range of Eigen Values that correspond to signal (using fit to Marcenko Pastur Distribution), we can then reset all other eigen values below this threshold to a constant value. This basically vaccums all the eigen values containing random noise into a single Eigen value. These transformed Eigen values are used with the original Eigen Vectors, to recreate a denoised covariance matrix, that will be used for further processing.
Compute using Python:
# File: Denoise-Constant-Residual-Eigen.py
import numpy as np,pandas as pd
import matplotlib.pyplot as plt
#---------------------------------------------------
def mpPDF(var,q,pts):
# Marcenko-Pastur pdf
# q=T/N
eMin,eMax=var*(1-(1./q)**.5)**2,var*(1+(1./q)**.5)**2
eVal=np.linspace(eMin,eMax,pts)
pdf=q/(2*np.pi*var*eVal)*((eMax-eVal)*(eVal-eMin))**.5
pdf=pd.Series(pdf,index=eVal)
return pdf
from sklearn.neighbors.kde import KernelDensity
#---------------------------------------------------
def getPCA(matrix):
# Get eVal,eVec from a Hermitian matrix
eVal,eVec=np.linalg.eigh(matrix)
indices=eVal.argsort()[::-1] # arguments for sorting eVal desc
eVal,eVec=eVal[indices],eVec[:,indices]
eVal=np.diagflat(eVal)
return eVal,eVec
#---------------------------------------------------
def fitKDE(obs,bWidth=.25,kernel='gaussian',x=None):
# Fit kernel to a series of obs, and derive the prob of obs
# x is the array of values on which the fit KDE will be evaluated
if len(obs.shape)==1:obs=obs.reshape(-1,1)
kde=KernelDensity(kernel=kernel,bandwidth=bWidth).fit(obs)
if x is None:x=np.unique(obs).reshape(-1,1)
if len(x.shape)==1:x=x.reshape(-1,1)
logProb=kde.score_samples(x) # log(density)
pdf=pd.Series(np.exp(logProb),index=x.flatten())
return pdf
# Determine matrix with some signal
def getRndCov(nCols,nFacts):
w=np.random.normal(size=(nCols,nFacts)) # We start with a fully random matrix
cov=np.dot(w,w.T) # random cov matrix, however not full rank # We calculate the covariance matrix of this fully random matrix
cov+=np.diag(np.random.uniform(size=nCols)) # full rank cov # We add some constant terms to the above fully random matrix, to represent some signal
return cov
#---------------------------------------------------
def cov2corr(cov):
# Derive the correlation matrix from a covariance matrix
std=np.sqrt(np.diag(cov)) # Determine standard deviation as square root of diagnonal terms (variances)
corr=cov/np.outer(std,std) # Dividing covariance by product of standard deviations will give you correlation matrix
corr[corr<-1],corr[corr>1]=-1,1 # numerical error
return corr
#---------------------------------------------------
alpha,nCols,nFact,q=.995,1000,100,10
cov=np.cov(np.random.normal(size=(nCols*q,nCols)),rowvar=0) # Generate a random matrix and calculate its covariance matrix
cov=alpha*cov+(1-alpha)*getRndCov(nCols,nFact) # noise+signal # Add this random matrix to a matrix containing some signal, with mixing ratio alpha
corr0=cov2corr(cov) # Generate Correlation matrix from Covariance matrix so we can apply Marcenko Pastur Theorem
eVal0,eVec0=getPCA(corr0) # Determine Eigen Values & Eigen Vectors of the Correlation matrix
from scipy.optimize import minimize # Importing the required optimizers
#---------------------------------------------------
def errPDFs(var,eVal,q,bWidth,pts=1000): # This is the objective function we will minimise
# Fit error
pdf0=mpPDF(var[0],q,pts) # theoretical pdf. For each test variance, the optimiser passes, the PDF for that variance, will be calculated using Marcenko Pastur Theorem (MPT)
pdf1=fitKDE(eVal,bWidth,x=pdf0.index.values) # empirical pdf. The actual data distribution calculated using the Gaussian KDE
sse=np.sum((pdf1-pdf0)**2) # Calculate difference between ideal random distribution, and actual data distribution. Evaluation is only conducted within range spanned PDF from MPT
return sse
#---------------------------------------------------
def findMaxEval(eVal,q,bWidth):
# Find max random eVal by fitting Marcenko’s distp
out=minimize(lambda *x:errPDFs(*x),.5,args=(eVal,q,bWidth),bounds=((1E-5,1-1E-5),)) # The minimiser is optimising the function errPDF (its objective function) by Varying variance between bounds with .5 as initial estimate
# .. contd The other parameters required by the objective function are passed separately but remain constant, and are not varied within the optimiser
if out['success']:var=out['x'][0] # If the optimiser converges, there is an in-built flag called success, which is queried, and the minima is obtained
else:var=1
eMax=var*(1+(1./q)**.5)**2 # For the minima variance, calculate the maximum eigen value as per MPT
return eMax,var
#---------------------------------------------------
eMax0,var0=findMaxEval(np.diag(eVal0),q,bWidth=.01) # Calculate var0 (degree of randomness in actual data), and max eigen value explained by randomness, eigen values above this are due to signal
nFacts0=eVal0.shape[0]-np.diag(eVal0)[::-1].searchsorted(eMax0) # identifying the number of eigen values that contain signal
# Denoise correlation matrix using Constant Residual Eigen Method
def denoisedCorr(eVal,eVec,nFacts):
eVal_=np.diag(eVal).copy() # Copy the original eigen values into a list
eVal_[nFacts:] = (np.sum(eVal_[nFacts:]))/float(eVal_.shape[0]-nFacts) # Eigen Values beyond NFacts are modified. This is because by default Eigen values are stored in
# decreasing order, and the non-random eigen values are always greater in value. All non-random eigen values are re-set to the average of their joint sum to preserve the
# trace of the matrix. This is not critical, any identifiable value will do. Notice that the Eigen vectors are not modified
eVal_ = np.diag(eVal_) # Change the Eigen value list back into a diagonal matrix
corr1 = np.dot(eVec,eVal_).dot(eVec.T) # Use Linear algebra properties to convert the Eigen values, and Eigen Vectors into a Covariance matrix
corr1 = cov2corr(corr1) # Convert covariance matrix to correlation matrix using earlier function
return corr1
corr1 = denoisedCorr(eVal0,eVec0,nFacts0) # The denoised correlation matrix corr1 can now be used for further processing
eVal1,eVec1 = getPCA(corr1) # Use the new correlation matrix to show that the underlying Eigen values and Eigen vectors were indeed modified.
# The above getPCA is for demonstration of this technique and is not required for further processing
# Plot the results of Constant Residual Eigen Method
fig,ax= plt.subplots()
plt.plot(np.diag(eVal0),'r',label='Original')
plt.plot(np.diag(eVal1),'b', label = 'Denoised')
ax.set_xlabel('Eigen Value #')
ax.set_ylabel('Eigen Value')
ax.set_title('Denoise Covariance Matrix using Constant Residual Eigen Method')
ax.legend()
plt.show()
Comments
Post a Comment