docs: def compute_Rk(L,A,n_samples):
      # TODO: need to inspect more into compute Rk.
      """
      Compute sparse L matrix and neighbors.
  
      Returns
      -------
      Rk_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.
      """
      laplacian_matrix = L.copy()
      laplacian_matrix.setdiag(0)
      laplacian_matrix.eliminate_zeros()
  
      n = n_samples
      Rk_tensor = []
      nbk = []
      row_A,column_A = A.T.nonzero()
  
      row,column = laplacian_matrix.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_A = np.diff(row_A).nonzero()[0]
      breaks_col = np.diff(sorted_col_vals).nonzero()[0]
  
      for k in range(n_samples):
          if k == 0:
              nbk.append( column_A[:breaks_row_A[k]+1].T )
              Rk_tensor.append(
                  nnz_val[np.sort(sorted_col_args[:breaks_col[k]+1])])
          elif k == n_samples-1:
              nbk.append( column_A[breaks_row_A[k-1]+1:].T )
              Rk_tensor.append(
                  nnz_val[np.sort(sorted_col_args[breaks_col[k-1]+1:])])
          else:
              nbk.append( column_A[breaks_row_A[k-1]+1:breaks_row_A[k]+1].T )
              Rk_tensor.append(nnz_val[np.sort(
                  sorted_col_args[breaks_col[k-1]+1:breaks_col[k]+1])])
  
      return Rk_tensor, nbk