<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def spd_eig(W, epsilon=1e-10, method='QR', canonical_signs=False):
 """ Rank-reduced eigenvalue decomposition of symmetric positive definite matrix. Removes all negligible eigenvalues Parameters W : ndarray((n, n), dtype=float) Symmetric positive-definite (spd) matrix. epsilon : float Truncation parameter. Eigenvalues with norms smaller than this cutoff will be removed. method : str Method to perform the decomposition of :math:`W` before inverting. Options are: * 'QR': QR-based robust eigenvalue decomposition of W * 'schur': Schur decomposition of W canonical_signs : boolean, default = False Fix signs in V, s. t. the largest element of in every row of V is positive. Returns ------- s : ndarray(k) k non-negligible eigenvalues, sorted by descending norms V : ndarray(n, k) k leading eigenvectors """

    # check input
    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)
    # compute the Eigenvalues of C0 using Schur factorization
    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) # sort them

    # determine the cutoff. We know that C0 is an spd matrix,
    # so we select the truncation threshold such that everything that is negative vanishes
    evmin = _np.min(s)
    if evmin < 0:
        epsilon = max(epsilon, -evmin + 1e-16)

    # determine effective rank m and perform low-rank approximations.
    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:
        # enforce canonical eigenvector signs
        for j in range(m):
            jj = _np.argmax(_np.abs(Vm[:, j]))
            Vm[:, j] *= _np.sign(Vm[jj, j])

    return sm, Vm