matrix.py 文件源码

def check_covmat(C,N=None,eps=1e-6):
    '''
    Verify that matrix M is a size NxN precision or covariance matirx
    '''
    if not type(C)==np.ndarray:
        raise ValueError("Covariance matrix should be a 2D numpy array")
    if not len(C.shape)==2:
        raise ValueError("Covariance matrix should be a 2D numpy array")
    if N is None: 
        N = C.shape[0]
    if not C.shape==(N,N):
        raise ValueError("Expected size %d x %d matrix"%(N,N))
    if np.any(~np.isreal(C)):
        raise ValueError("Covariance matrices should not contain complex numbers")
    C = np.real(C)
    if np.any(~np.isfinite(C)):
        raise ValueError("Covariance matrix contains NaN or ±inf!")
    if not np.all(np.abs(C-C.T)<eps):
        raise ValueError("Covariance matrix is not symmetric up to precision %0.1e"%eps)

    # Get just highest eigenvalue
    maxe = np.real(scipy.linalg.decomp.eigh(C,eigvals=(N-1,N-1))[0][0])

    # Get just lowest eigenvalue
    mine = np.real(scipy.linalg.decomp.eigh(C,eigvals=(0,0))[0][0])

    #if np.any(w<-eps):
    #    raise ValueError('Covariance matrix contains eigenvalue %0.3e<%0.3e'%(np.min(w),-eps)) 
    if mine<0:
        raise ValueError('Covariance matrix contains negative eigenvalue %0.3e'%mine) 
    if (mine<eps):
        C = C + eye(N)*(eps-mine)
    # trucate spectrum at some small value
    # w[w<eps]=eps
    # Very large eigenvalues can also cause numeric problems
    # w[w>1./eps]=1./eps;
    # maxe = np.max(np.abs(w))
    # if maxe>10./eps:
    #     raise ValueError('Covariance matrix eigenvalues %0.2e is larger than %0.2e'%(maxe,10./eps))
    # Rebuild matrix
    # C = v.dot(np.diag(w)).dot(v.T)
    # Ensure symmetry (only occurs as a numerical error for very large matrices?)
    C = 0.5*(C+C.T)
    return C

# need a faster covariance matrix checker