Примечание

Перейти в конец, чтобы скачать полный код примера или запустить этот пример в браузере через JupyterLite или Binder.

Граница Джонсона — Линденштраусса для вложения с помощью случайных проекций

Лемма Джонсона — Линденштраусса утверждает, что любой набор данных большой размерности можно случайным образом спроецировать в евклидово пространство меньшей размерности, контролируя искажение попарных расстояний.

# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause

import sys
from time import time

import matplotlib.pyplot as plt
import numpy as np

from sklearn.datasets import fetch_20newsgroups_vectorized, load_digits
from sklearn.metrics.pairwise import euclidean_distances
from sklearn.random_projection import (
    SparseRandomProjection,
    johnson_lindenstrauss_min_dim,
)

Теоретические границы

Искажение, вносимое случайной проекцией p, ограничивается тем фактом, что p задаёт ε-вложение с высокой вероятностью, как определяется условием:

\[(1 - eps) \|u - v\|^2 < \|p(u) - p(v)\|^2 < (1 + eps) \|u - v\|^2\]

Здесь u и v — любые строки набора данных формы (n_samples, n_features), а p — проекция с помощью случайной гауссовской матрицы N(0, 1) формы (n_components, n_features) (или разреженной матрицы Ахлиоптаса).

Минимальное число компонент, необходимое для гарантии ε-вложения, задаётся формулой:

\[n\_components \geq 4 log(n\_samples) / (eps^2 / 2 - eps^3 / 3)\]

На первом графике показано, что при увеличении числа образцов n_samples минимальное число измерений n_components логарифмически возрастает, чтобы гарантировать ε-вложение eps.

# range of admissible distortions
eps_range = np.linspace(0.1, 0.99, 5)
colors = plt.cm.Blues(np.linspace(0.3, 1.0, len(eps_range)))

# range of number of samples (observation) to embed
n_samples_range = np.logspace(1, 9, 9)

plt.figure()
for eps, color in zip(eps_range, colors):
    min_n_components = johnson_lindenstrauss_min_dim(n_samples_range, eps=eps)
    plt.loglog(n_samples_range, min_n_components, color=color)

plt.legend([f"eps = {eps:0.1f}" for eps in eps_range], loc="lower right")
plt.xlabel("Number of observations to eps-embed")
plt.ylabel("Minimum number of dimensions")
plt.title("Johnson-Lindenstrauss bounds:\nn_samples vs n_components")
plt.show()

Johnson-Lindenstrauss bounds: n_samples vs n_components

На втором графике показано, что увеличение допустимого искажения eps позволяет значительно уменьшить минимальное число измерений n_components при заданном числе образцов n_samples

# range of admissible distortions
eps_range = np.linspace(0.01, 0.99, 100)

# range of number of samples (observation) to embed
n_samples_range = np.logspace(2, 6, 5)
colors = plt.cm.Blues(np.linspace(0.3, 1.0, len(n_samples_range)))

plt.figure()
for n_samples, color in zip(n_samples_range, colors):
    min_n_components = johnson_lindenstrauss_min_dim(n_samples, eps=eps_range)
    plt.semilogy(eps_range, min_n_components, color=color)

plt.legend([f"n_samples = {n}" for n in n_samples_range], loc="upper right")
plt.xlabel("Distortion eps")
plt.ylabel("Minimum number of dimensions")
plt.title("Johnson-Lindenstrauss bounds:\nn_components vs eps")
plt.show()

Johnson-Lindenstrauss bounds: n_components vs eps

Эмпирическая проверка

Проверим приведённые выше границы на текстовом наборе документов 20 newsgroups (частоты слов TF-IDF) или на наборе данных цифр:

для набора данных 20 newsgroups около 300 документов, содержащих в общей сложности 100 тыс. признаков, с помощью разреженной случайной матрицы проецируются в евклидовы пространства меньшей размерности с различными целевыми числами измерений n_components.
для набора данных цифр данные об уровнях серого для пикселей изображений 300 рукописных цифр размером 8×8 случайным образом проецируются в пространства с различными, более высокими значениями размерности n_components.

По умолчанию используется набор данных 20 newsgroups. Чтобы запустить пример на наборе данных цифр, передайте этому скрипту аргумент командной строки --use-digits-dataset.

if "--use-digits-dataset" in sys.argv:
    data = load_digits().data[:300]
else:
    data = fetch_20newsgroups_vectorized().data[:300]

Для каждого значения n_components строим графики:

двумерное распределение пар образцов, где попарные расстояния в исходном и спроецированном пространствах представлены соответственно по осям x и y;
одномерную гистограмму отношения этих расстояний (спроецированное / исходное).

n_samples, n_features = data.shape
print(
    f"Embedding {n_samples} samples with dim {n_features} using various "
    "random projections"
)

n_components_range = np.array([300, 1_000, 10_000])
dists = euclidean_distances(data, squared=True).ravel()

# select only non-identical samples pairs
nonzero = dists != 0
dists = dists[nonzero]

for n_components in n_components_range:
    t0 = time()
    rp = SparseRandomProjection(n_components=n_components)
    projected_data = rp.fit_transform(data)
    print(
        f"Projected {n_samples} samples from {n_features} to {n_components} in "
        f"{time() - t0:0.3f}s"
    )
    if hasattr(rp, "components_"):
        n_bytes = rp.components_.data.nbytes
        n_bytes += rp.components_.indices.nbytes
        print(f"Random matrix with size: {n_bytes / 1e6:0.3f} MB")

    projected_dists = euclidean_distances(projected_data, squared=True).ravel()[nonzero]

    plt.figure()
    min_dist = min(projected_dists.min(), dists.min())
    max_dist = max(projected_dists.max(), dists.max())
    plt.hexbin(
        dists,
        projected_dists,
        gridsize=100,
        cmap=plt.cm.PuBu,
        extent=[min_dist, max_dist, min_dist, max_dist],
    )
    plt.xlabel("Pairwise squared distances in original space")
    plt.ylabel("Pairwise squared distances in projected space")
    plt.title("Pairwise distances distribution for n_components=%d" % n_components)
    cb = plt.colorbar()
    cb.set_label("Sample pairs counts")

    rates = projected_dists / dists
    print(f"Mean distances rate: {np.mean(rates):.2f} ({np.std(rates):.2f})")

    plt.figure()
    plt.hist(rates, bins=50, range=(0.0, 2.0), edgecolor="k", density=True)
    plt.xlabel("Squared distances rate: projected / original")
    plt.ylabel("Distribution of samples pairs")
    plt.title("Histogram of pairwise distance rates for n_components=%d" % n_components)

    # TODO: compute the expected value of eps and add them to the previous plot
    # as vertical lines / region

plt.show()

Embedding 300 samples with dim 130107 using various random projections
Projected 300 samples from 130107 to 300 in 0.262s
Random matrix with size: 1.293 MB
Mean distances rate: 1.01 (0.16)
Projected 300 samples from 130107 to 1000 in 0.850s
Random matrix with size: 4.333 MB
Mean distances rate: 1.00 (0.09)
Projected 300 samples from 130107 to 10000 in 7.971s
Random matrix with size: 43.239 MB
Mean distances rate: 1.01 (0.03)

Видно, что при малых значениях n_components распределение широкое, содержит много искажённых пар и асимметрично (из-за жёсткой нижней границы отношения, равной нулю: расстояния всегда положительны). При больших значениях n_components искажение контролируется, и случайная проекция хорошо сохраняет расстояния.

Примечания

Согласно лемме JL, для проецирования 300 образцов без значительного искажения потребуется не менее нескольких тысяч измерений, независимо от числа признаков исходного набора данных.

Поэтому использование случайных проекций для набора данных цифр, содержащего всего 64 признака во входном пространстве, не имеет смысла: в данном случае оно не позволяет уменьшить размерность.

С другой стороны, для набора данных twenty newsgroups размерность можно уменьшить с 56 436 до 10 000, сохранив попарные расстояния с приемлемой точностью.

Общее время выполнения скрипта: (0 минут 11.304 секунды)

Граница Джонсона — Линденштраусса для вложения с помощью случайных проекций

Теоретические границы

Эмпирическая проверка

Примечания

Похожие примеры