Ограничение Джонсона-Линденбаума для встраивания с помощью случайных проекций

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

Искажение, введенное случайной проекцией p , подтверждается тем, что p определяет вложение eps с хорошей вероятностью, как определено:

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

Минимальное количество компонентов, гарантирующих вложение eps, задается формулой:

На первом графике показано, что с увеличением количества выборок 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()

На втором графике показано, что увеличение допустимого искажения 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()

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

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

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

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

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

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

• 2D распределение пар образцов с парами расстояний в исходном и проецированном пространствах соответственно как оси x и y.
• 1D гистограмму отношения этих расстояний (проектированное / исходное).

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.224s
Random matrix with size: 1.293 MB
Mean distances rate: 1.01 (0.17)
Projected 300 samples from 130107 to 1000 in 0.782s
Random matrix with size: 4.323 MB
Mean distances rate: 1.00 (0.09)
Projected 300 samples from 130107 to 10000 in 7.572s
Random matrix with size: 43.268 MB
Mean distances rate: 1.01 (0.03)

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

Замечания

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

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

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

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

Download Jupyter notebook: plot_johnson_lindenstrauss_bound.ipynb

Download Python source code: plot_johnson_lindenstrauss_bound.py

Download zipped: plot_johnson_lindenstrauss_bound.zip

Связанные примеры

Обучение многообразию на рукописных цифрах: Локально-линейное вложение, Isomap…

Демонстрация алгоритма кластеризации HDBSCAN

Оценка Ледоита-Вольфа против OAS

Преобразование данных в нормальное распределение