def spd_eig(W, epsilon=1e-10, method='QR', canonical_signs=False):
    assert _np.allclose(W.T, W), 'W is not a symmetric matrix'
    if method.lower() == 'qr':
        from .eig_qr.eig_qr import eig_qr
        s, V = eig_qr(W)
    elif method.lower() == 'schur':
        from scipy.linalg import schur
        S, V = schur(W)
        s = _np.diag(S)
    else:
        raise ValueError('method not implemented: ' + method)
    s, V = sort_by_norm(s, V)
    evmin = _np.min(s)
    if evmin < 0:
        epsilon = max(epsilon, -evmin + 1e-16)
    evnorms = _np.abs(s)
    n = _np.shape(evnorms)[0]
    m = n - _np.searchsorted(evnorms[::-1], epsilon)
    if m == 0:
        raise _ZeroRankError('All eigenvalues are smaller than %g, rank reduction would discard all dimensions.'%epsilon)
    Vm = V[:, 0:m]
    sm = s[0:m]
    if canonical_signs:
        for j in range(m):
            jj = _np.argmax(_np.abs(Vm[:, j]))
            Vm[:, j] *= _np.sign(Vm[jj, j])
    return sm, Vm