<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def interpolate(G, f_subsampled, keep_inds, order=100, reg_eps=0.005, **kwargs):
 r"""Interpolate a graph signal. Parameters G : Graph f_subsampled : ndarray A graph signal on the graph G. keep_inds : ndarray List of indices on which the signal is sampled. order : int Degree of the Chebyshev approximation (default = 100). reg_eps : float The regularized graph Laplacian is $\bar{L}=L+\epsilon I$. A smaller epsilon may lead to better regularization, but will also require a higher order Chebyshev approximation. Returns ------- f_interpolated : ndarray Interpolated graph signal on the full vertex set of G. References See :cite:`pesenson2009variational` """

    L_reg = G.L + reg_eps * sparse.eye(G.N)
    K_reg = getattr(G.mr, 'K_reg', kron_reduction(L_reg, keep_inds))
    green_kernel = getattr(G.mr, 'green_kernel',
                           filters.Filter(G, lambda x: 1. / (reg_eps + x)))

    alpha = K_reg.dot(f_subsampled)

    try:
        Nv = np.shape(f_subsampled)[1]
        f_interpolated = np.zeros((G.N, Nv))
    except IndexError:
        f_interpolated = np.zeros((G.N))

    f_interpolated[keep_inds] = alpha

    return _analysis(green_kernel, f_interpolated, order=order, **kwargs)