A square matrix is PSD if all eigenvalues of its Hermitian part are
    non- negative. The Hermitian part is given by (self + M*)/2, where M* is
    the complex conjugate transpose of M

def check_psd(matrix, tolerance=1e-6):
    """ A square matrix is PSD if all eigenvalues of its Hermitian part are
    non- negative. The Hermitian part is given by (self + M*)/2, where M* is
    the complex conjugate transpose of M """

    hermitian = (matrix + matrix.T.conjugate()) / 2
    eigenvalues = np.linalg.eigh(hermitian)[0]
    return (eigenvalues > -tolerance).all()