def compute_Lk(laplacian_matrix,n_samples,subset):
Lk_tensor = []<EOL>nbk = []<EOL>row,column = laplacian_matrix.T.nonzero()<EOL>nnz_val = np.squeeze(np.asarray(laplacian_matrix.T[(row,column)]))<EOL>sorted_col_args = np.argsort(column)<EOL>sorted_col_vals = column[sorted_col_args]<EOL>breaks_row = np.diff(row).nonzero()[<NUM_LIT:0>]<EOL>breaks_col = np.diff(sorted_col_vals).nonzero()[<NUM_LIT:0>]<EOL>si_map = {}<EOL>for idx,k in enumerate(subset):<EOL><INDENT>if k == <NUM_LIT:0>:<EOL><INDENT>nbk.append( column[:breaks_row[k]+<NUM_LIT:1>].T )<EOL>lk = nnz_val[np.sort(sorted_col_args[:breaks_col[k]+<NUM_LIT:1>])]<EOL><DEDENT>elif k == n_samples-<NUM_LIT:1>:<EOL><INDENT>nbk.append( column[breaks_row[k-<NUM_LIT:1>]+<NUM_LIT:1>:].T )<EOL>lk = nnz_val[np.sort(sorted_col_args[breaks_col[k-<NUM_LIT:1>]+<NUM_LIT:1>:])]<EOL><DEDENT>else:<EOL><INDENT>nbk.append( column[breaks_row[k-<NUM_LIT:1>]+<NUM_LIT:1>:breaks_row[k]+<NUM_LIT:1>].T )<EOL>lk = nnz_val[np.sort(<EOL>sorted_col_args[breaks_col[k-<NUM_LIT:1>]+<NUM_LIT:1>:breaks_col[k]+<NUM_LIT:1>])]<EOL><DEDENT>npair = nbk[idx].shape[<NUM_LIT:0>]<EOL>rk = (nbk[idx] == k).nonzero()[<NUM_LIT:0>]<EOL>Lk = sp.sparse.lil_matrix((npair,npair))<EOL>Lk.setdiag(lk)<EOL>Lk[:,rk] = -(lk.reshape(-<NUM_LIT:1>,<NUM_LIT:1>))<EOL>Lk[rk,:] = -(lk.reshape(<NUM_LIT:1>,-<NUM_LIT:1>))<EOL>Lk_tensor.append(sp.sparse.csr_matrix(Lk))<EOL>si_map[k] = idx<EOL><DEDENT>assert len(Lk_tensor) == subset.shape[<NUM_LIT:0>],'<STR_LIT>'<EOL>return Lk_tensor, nbk, si_map<EOL>
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.