ISDF decomposition using pivoted Cholesky, generalized for tensors.
:param rho: Input tensor (N_b, N_grid, *trailing_dims) or (N_b, N_grid)
:param n_rank: Number of interpolation points
Returns:
Selected columns from rho (N_b, n_rank, *trailing_dims) or (N_b, n_rank)
xi: Interpolation coefficients (n_rank, N_grid, *trailing_dims) or (n_rank, N_grid)