docs: def null_space(M, k, k_skip=1, eigen_solver='arpack',
                 random_state=None, solver_kwds=None):
      """
      Find the null space of a matrix M: eigenvectors associated with 0 eigenvalues
  
      Parameters
      ----------
      M : {array, matrix, sparse matrix, LinearOperator}
          Input covariance matrix: should be symmetric positive semi-definite
      k : integer
          Number of eigenvalues/vectors to return
      k_skip : integer, optional
          Number of low eigenvalues to skip.
      eigen_solver : {'auto', 'dense', 'arpack', 'lobpcg', or 'amg'}
          'auto' :
              algorithm will attempt to choose the best method for input data
          'dense' :
              use standard dense matrix operations for the eigenvalue decomposition.
              For this method, M must be an array or matrix type.  This method should be avoided for large problems.
          'arpack' :
              use arnoldi iteration in shift-invert mode. For this method,
              M may be a dense matrix, sparse matrix, or general linear operator.
              Warning: ARPACK can be unstable for some problems.  It is best to
              try several random seeds in order to check results.
          'lobpcg' :
              Locally Optimal Block Preconditioned Conjugate Gradient Method.
              A preconditioned eigensolver for large symmetric positive definite
              (SPD) generalized eigenproblems.
          'amg' :
              AMG requires pyamg to be installed. It can be faster on very large,
              sparse problems, but may also lead to instabilities.
      random_state: numpy.RandomState or int, optional
          The generator or seed used to determine the starting vector for arpack
          iterations.  Defaults to numpy.random.
      solver_kwds : any additional keyword arguments to pass to the selected eigen_solver
  
      Returns
      -------
      null_space : estimated k vectors of the null space
      error : estimated error (sum of eigenvalues)
  
      Notes
      -----
      dense solver key words: see
          http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.eigh.html
          for symmetric problems and
          http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.eig.html#scipy.linalg.eig
          for non symmetric problems.
      arpack sovler key words: see
          http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.eigsh.html
          for symmetric problems and http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.eigs.html#scipy.sparse.linalg.eigs
          for non symmetric problems.
      lobpcg solver keywords: see
          http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lobpcg.html
      amg solver keywords: see
          http://pyamg.googlecode.com/svn/branches/1.0.x/Docs/html/pyamg.aggregation.html#module-pyamg.aggregation.aggregation
          (Note amg solver uses lobpcg and also accepts lobpcg keywords)
      """
      eigen_solver, solver_kwds = check_eigen_solver(eigen_solver, solver_kwds,
                                                     size=M.shape[0],
                                                     nvec=k + k_skip)
      random_state = check_random_state(random_state)
  
      if eigen_solver == 'arpack':
          # This matches the internal initial state used by ARPACK
          v0 = random_state.uniform(-1, 1, M.shape[0])
          try:
              eigen_values, eigen_vectors = eigsh(M, k + k_skip, sigma=0.0,
                                                  v0=v0,**(solver_kwds or {}))
          except RuntimeError as msg:
              raise ValueError("Error in determining null-space with ARPACK. "
                               "Error message: '%s'. "
                               "Note that method='arpack' can fail when the "
                               "weight matrix is singular or otherwise "
                               "ill-behaved.  method='dense' is recommended. "
                               "See online documentation for more information."
                               % msg)
  
          return eigen_vectors[:, k_skip:], np.sum(eigen_values[k_skip:])
      elif eigen_solver == 'dense':
          if hasattr(M, 'toarray'):
              M = M.toarray()
          eigen_values, eigen_vectors = eigh(M, eigvals=(0, k+k_skip),overwrite_a=True,
                                             **(solver_kwds or {}))
          index = np.argsort(np.abs(eigen_values))
          eigen_vectors = eigen_vectors[:, index]
          eigen_values = eigen_values[index]
          return eigen_vectors[:, k_skip:k+1], np.sum(eigen_values[k_skip:k+1])
          # eigen_values, eigen_vectors = eigh(
              # M, eigvals=(k_skip, k + k_skip - 1), overwrite_a=True)
          # index = np.argsort(np.abs(eigen_values))
          # return eigen_vectors[:, index], np.sum(eigen_values)
      elif (eigen_solver == 'amg' or eigen_solver == 'lobpcg'):
          # M should be positive semi-definite. Add 1 to make it pos. def.
          try:
              M = sparse.identity(M.shape[0]) + M
              n_components = min(k + k_skip + 10, M.shape[0])
              eigen_values, eigen_vectors = eigen_decomposition(M, n_components,
                                                                eigen_solver = eigen_solver,
                                                                drop_first = False,
                                                                largest = False,
                                                                random_state=random_state,
                                                                solver_kwds=solver_kwds)
              eigen_values = eigen_values -1
              index = np.argsort(np.abs(eigen_values))
              eigen_values = eigen_values[index]
              eigen_vectors = eigen_vectors[:, index]
              return eigen_vectors[:, k_skip:k+1], np.sum(eigen_values[k_skip:k+1])
          except np.linalg.LinAlgError: # try again with bigger increase
              warnings.warn("LOBPCG failed the first time. Increasing Pos Def adjustment.")
              M = 2.0*sparse.identity(M.shape[0]) + M
              n_components = min(k + k_skip + 10, M.shape[0])
              eigen_values, eigen_vectors = eigen_decomposition(M, n_components,
                                                                eigen_solver = eigen_solver,
                                                                drop_first = False,
                                                                largest = False,
                                                                random_state=random_state,
                                                                solver_kwds=solver_kwds)
              eigen_values = eigen_values - 2
              index = np.argsort(np.abs(eigen_values))
              eigen_values = eigen_values[index]
              eigen_vectors = eigen_vectors[:, index]
              return eigen_vectors[:, k_skip:k+1], np.sum(eigen_values[k_skip:k+1])
      else:
          raise ValueError("Unrecognized eigen_solver '%s'" % eigen_solver)