code: def spd_inv_split(W, epsilon=1e-10, method='QR', canonical_signs=False):
      
      if (_np.shape(W)[0] == 1):
          if W[0,0] < epsilon:
              raise _ZeroRankError(
                  'All eigenvalues are smaller than %g, rank reduction would discard all dimensions.' % epsilon)
          L = 1./_np.sqrt(W[0,0])
      else:
          sm, Vm = spd_eig(W, epsilon=epsilon, method=method, canonical_signs=canonical_signs)
          L = _np.dot(Vm, _np.diag(1.0/_np.sqrt(sm)))
  
      # return split
      return L
docs: Compute :math:`W^{-1} = L L^T` of the symmetric positive-definite matrix :math:`W`.
  
      by first reducing W to a low-rank approximation that is truly spd.
  
      Parameters
      ----------
      W : ndarray((m,m), 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 L, s. t. the largest element of in every row of L is positive.
  
      Returns
      -------
      L : ndarray((n, r))
          Matrix :math:`L` from the decomposition :math:`W^{-1} = L L^T`.