docs: def compute_Lk(laplacian_matrix,n_samples,subset):
      """
      Compute sparse L matrix, neighbors and subset to L matrix index map.
  
      Returns
      -------
      Lk_tensor : array-like. Length = n
          each component correspond to the sparse matrix of Lk, which is
          generated by extracting the kth row of laplacian and removing zeros.
      nbk : array-like. Length = n
          each component correspond to the neighbor index of point k, which is
          used in slicing the gradient, Y or S arrays.
      si_map : dictionary.
          subset index to Lk_tensor (or nbk) index mapping.
      """
      Lk_tensor = []
      nbk = []
      row,column = laplacian_matrix.T.nonzero()
      nnz_val = np.squeeze(np.asarray(laplacian_matrix.T[(row,column)]))
      sorted_col_args = np.argsort(column)
      sorted_col_vals = column[sorted_col_args]
  
      breaks_row = np.diff(row).nonzero()[0]
      breaks_col = np.diff(sorted_col_vals).nonzero()[0]
  
      si_map = {}
  
      for idx,k in enumerate(subset):
          if k == 0:
              nbk.append( column[:breaks_row[k]+1].T )
              lk = nnz_val[np.sort(sorted_col_args[:breaks_col[k]+1])]
          elif k == n_samples-1:
              nbk.append( column[breaks_row[k-1]+1:].T )
              lk = nnz_val[np.sort(sorted_col_args[breaks_col[k-1]+1:])]
          else:
              nbk.append( column[breaks_row[k-1]+1:breaks_row[k]+1].T )
              lk = nnz_val[np.sort(
                  sorted_col_args[breaks_col[k-1]+1:breaks_col[k]+1])]
  
          npair = nbk[idx].shape[0]
          rk = (nbk[idx] == k).nonzero()[0]
          Lk = sp.sparse.lil_matrix((npair,npair))
          Lk.setdiag(lk)
          Lk[:,rk] = -(lk.reshape(-1,1))
          Lk[rk,:] = -(lk.reshape(1,-1))
          Lk_tensor.append(sp.sparse.csr_matrix(Lk))
          si_map[k] = idx
  
      assert len(Lk_tensor) == subset.shape[0], \
          'Size of Lk_tensor should be the same as subset.'
      return Lk_tensor, nbk, si_map