Function kernel_nystroem

ot/lowrank.py:530–587

r""" Compute left and right factors corresponding to the Nystroem method on the Gaussian kernel :math:`K(x^s_i, x^t_j) = \exp(-\|x^s_i-x^t_j\|^2/2\sigma^2)`. The Nystroem approximation is computed by sampling :math:`\min(n, \lceil(c / 2))\rceil)` components in each distribution, where :m

(X_s, X_t, anchors=50, sigma=1.0, random_state=None)

def kernel_nystroem(X_s, X_t, anchors=50, sigma=1.0, random_state=None):
    r"""
    Compute left and right factors corresponding to the Nystroem method on the Gaussian kernel :math:`K(x^s_i, x^t_j) = \exp(-\|x^s_i-x^t_j\|^2/2\sigma^2)`.
    The Nystroem approximation is computed by sampling :math:`\min(n, \lceil(c / 2))\rceil)` components in each distribution,
    where :math:`n` is the number of samples in the distribution and :math:`c` the total number of anchor points.

    Parameters
    ----------
    X_s : array-like, shape (n_samples_a, dim)
        samples in the source domain
    X_t : array-like, shape (n_samples_b, dim)
        samples in the target domain
    anchors : int, optional
        The total number of anchors sampled for the Nystroem approximation (anchors/2 in each distribution), default 50.
    sigma : float, optional
        The standard deviation parameter for the Gaussian kernel.
    random_state : int, optional
        The random state for sampling the components in each distribution.

    Returns
    -------
    left_factor : array-like, shape (n_samples_a, dim_r)
        Left factor of Nystroem
    right_factor : array-like, shape (n_samples_b, dim_r)
        Right factor of Nystroem
    """
    nx = get_backend(X_s, X_t)
    if random_state is not None:
        nx.seed(random_state)

    n, m = X_s.shape[0], X_t.shape[0]

    if not isinstance(anchors, int) or anchors < 2:
        raise ValueError("anchors must be an integer >= 2")

    n_components_source = min(n, anchors // 2)
    n_components_target = min(m, anchors // 2)
    # draw n_components/2 points in each distribution
    inds_source = nx.randperm(n)  # sample n_components_source uniformly
    basis_source = X_s[inds_source[:n_components_source]]

    inds_target = nx.randperm(m)
    basis_target = X_t[inds_target[:n_components_target]]

    basis = nx.concatenate((basis_source, basis_target))

    Mzz = dist(basis, metric="sqeuclidean")  # compute \|z_i - z_j\|_2^2
    basis_kernel = nx.exp(-Mzz / (2.0 * sigma**2))

    normalization = nx.pinv(basis_kernel, hermitian=True)

    Mxz = dist(X_s, basis, metric="sqeuclidean")
    Myz = dist(X_t, basis, metric="sqeuclidean")

    left_factor = nx.exp(-Mxz / (2.0 * sigma**2)) @ normalization
    right_factor = nx.exp(-Myz / (2.0 * sigma**2))

    return left_factor, right_factor  # left_factor @ right_factor.T approx K