code: 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