Использование квантового отжига для выбора признаков в scikit-learn

Выбор функций для моделей scikit-learn для наборов данных с множеством функций с использованием квантовой обработки

Выбор функций — обширная тема машинного обучения. Если все сделано правильно, это может помочь уменьшить переоснащение, повысить интерпретируемость, снизить вычислительную нагрузку и т. Д. Для выполнения выбора признаков используются многочисленные методы. Их всех объединяет то, что они смотрят на набор функций и пытаются отделить функции, которые приводят к хорошим результатам (точные прогнозы, интерпретируемые модели и т. д.), от функций, которые не помогают в достижении этой цели.

Особенно сложными являются случаи, когда количество признаков очень велико. Исчерпывающее исследование всех комбинаций функций часто непомерно вычислительно. Такие алгоритмы, как пошаговый поиск по функции regsubsets() в R, могут попытаться вывести многообещающие комбинации признаков, добавляя/удаляя признаки и сравнивая результаты. Но в конечном счете, когда количество признаков велико, остается компромисс между степенью успеха поиска и усилиями, затраченными на поиск наилучшего сочетания признаков.

Как показывает опыт, когда лучшее решение проблемы включает в себя поиск по большому количеству комбинаций, возможно, стоит изучить квантовый отжиг. Я покажу пример выбора функций для набора данных с сотнями функций с использованием плагина scikit-learn, недавно опубликованного D-Wave.

D-Wave и scikit-learn

Имейте в виду, что это не универсальные квантовые вычисления модели шлюза. Это алгоритм, который, по сути, похож на имитацию отжига тем, что есть целевая функция, и что-то вроде имитации отжига используется для нахождения комбинации значений, минимизирующей цель. За исключением того, что отжиг не моделируется — вместо этого реальная система программируется таким образом, чтобы физическая энергия системы соответствовала целевой функции. Энергия системы снижается до минимума (отжиг), а затем решением является просто состояние системы, которое считывается и возвращается пользователю.

D-Wave создает квантовые отжиги, которые могут эффективно решать многие задачи с высокой комбинаторной сложностью. Пока вы можете свести проблему к бинарной квадратичной модели (BQM), или BQM с ограничениями (CQM), или некоторому дискретному обобщению вышеизложенного (DQM), проблема может быть передана квантовым решателям. Подробнее в документации.

Теоретически вы могли бы сформулировать алгоритм выбора признаков в терминах BQM, где наличие признака — это бинарная переменная со значением 1, а отсутствие признака — переменная, равная 0, но это требует некоторых усилий. D-Wave предоставляет плагин scikit-learn, который можно подключить непосредственно к пайплайнам scikit-learn и упростить процесс.

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

Но если все, что вам нужно, это выбрать лучшие функции в наборе данных, достаточно простого вызова метода SelectFromQuadraticModel(). Это сворачивает весь алгоритм в одну строку кода. Мы покажем это в конце.

Проблема и метод

Рассмотрим Набор данных распознавания сцен, загруженный из OpenML, первоначально созданный Бутелл, М., Луо, Дж., Шен, X., Браун, К. (2004) и распространяемый его авторами по лицензии CC BY. Это набор данных бинарной классификации с одной зависимой переменной (Городской) и имеет около 300 признаков. Значение переменной ответа (бинарное: городской или не городской) необходимо прогнозировать на основе комбинации признаков. Простая модель классификации, такая как RandomForestClassifier(), должна работать хорошо, если признаки выбраны правильно.

>>> dataset.get_data()[0]

         attr1     attr2     attr3     attr4     attr5     attr6     attr7     attr8  ...   attr293   attr294  Beach  Sunset  FallFoliage  Field  Mountain  Urban
0     0.646467  0.666435  0.685047  0.699053  0.652746  0.407864  0.150309  0.535193  ...  0.014025  0.029709      1       0            0      0         1      0
1     0.770156  0.767255  0.761053  0.745630  0.742231  0.688086  0.708416  0.757351  ...  0.082672  0.036320      1       0            0      0         0      1
2     0.793984  0.772096  0.761820  0.762213  0.740569  0.734361  0.722677  0.849128  ...  0.112506  0.083924      1       0            0      0         0      0
3     0.938563  0.949260  0.955621  0.966743  0.968649  0.869619  0.696925  0.953460  ...  0.049780  0.090959      1       0            0      0         0      0
4     0.512130  0.524684  0.520020  0.504467  0.471209  0.417654  0.364292  0.562266  ...  0.164270  0.184290      1       0            0      0         0      0
...        ...       ...       ...       ...       ...       ...       ...       ...  ...       ...       ...    ...     ...          ...    ...       ...    ...
2402  0.875782  0.901653  0.926227  0.721366  0.795826  0.867642  0.794125  0.899067  ...  0.254413  0.134350      0       0            0      0         0      1
2403  0.657706  0.669877  0.692338  0.713920  0.727374  0.750354  0.684372  0.718770  ...  0.048747  0.041638      0       0            0      0         0      1
2404  0.952281  0.944987  0.905556  0.836604  0.875916  0.957034  0.953938  0.967956  ...  0.017547  0.019734      0       0            0      0         0      1
2405  0.883990  0.899004  0.901019  0.904298  0.846402  0.858145  0.851362  0.852472  ...  0.226332  0.223070      0       0            0      0         0      1
2406  0.974915  0.866425  0.818144  0.936140  0.938583  0.935087  0.930597  1.000000  ...  0.025059  0.004033      0       0            0      0         0      1

[2407 rows x 300 columns]

Метод выбора признаков, описанный в статье Милна, Раунда и Годдарда (2018), хорошо подходит для квантовых отжигов и может быть резюмирован следующим образом:

Разделите набор данных на матрицу признаков, а вектор ответа — знакомый по scikit-learn кортеж (X, y). Используя их, постройте матрицу корреляции C, где Cii представляют собой корреляции признаков с откликом, а Cij — взаимные корреляции признаков. Из функций, которые сильно коррелируют с ответом, мы хотим захватить как можно больше. Из функций, которые коррелируют друг с другом, мы хотим захватить как можно меньше, но не затрагивая функции, которые сильно коррелируют с откликом. Очевидно, мы хотим найти баланс между всеми этими критериями.

Пусть Xi будет бинарными переменными, которые указывают, выбран ли объект i для окончательного набора данных. Если выбран признак i, то Xi = 1, иначе Xi = 0. Тогда проблема становится эквивалентной поиску экстремума:

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

Целевая функция буквально представляет собой бинарную квадратичную модель — BQM. Он бинарный, потому что Xi может быть либо 0, либо 1. Он квадратичный, потому что членами высшего порядка являются квадратичные взаимодействия. Эту проблему можно решить с помощью квантового отжига. Единственное ограничение, которое нам нужно будет применить, заключается в том, что количество переменных Xi, равных 1, не может быть больше, чем общее количество признаков, которые мы готовы принять.

Трудный путь выбора функций

Давайте решим проблему. Блок кода ниже импортирует все необходимые библиотеки, загружает набор данных и создает экземпляр модели классификатора.

import numpy as np
import openml
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
import plotly.express as px
from plotly.subplots import make_subplots
import dimod
from dwave.system import LeapHybridCQMSampler

dataset = openml.datasets.get_dataset(312)

X, y, categorical_indicator, attribute_names = dataset.get_data(
    target=dataset.default_target_attribute, dataset_format="dataframe"
)
X = X.astype(float)
y = y.values.astype(int)

clf = RandomForestClassifier()

Нас интересуют две вещи:

«релевантность» каждой функции, которая представляет собой силу ее корреляции с ответом
«избыточность» признаков, которая представляет собой просто набор недиагональных элементов в корреляционной матрице (веса квадратичных членов взаимодействия)

Давайте создадим пару функций для визуализации релевантности и избыточности, применим их ко всему набору данных без выбора признаков, а затем проверим (в 5 раз) производительность модели на всем наборе данных.

def evaluate_model(m, X, y, indices=None):
    if not indices is None:
        X_filtered = X.iloc[:, indices]
    else:
        X_filtered = X
    acc = np.mean(cross_val_score(clf, X_filtered, y, cv=5))
    return acc


def show_relevance_redundancy(X, y, indices=None, title=""):
    if not indices is None:
        X = X.iloc[:, indices].copy()
    y = y
    fig = make_subplots(
        rows=1,
        cols=2,
        column_widths=[0.68, 0.32],
        column_titles=["relevance", "redundancy"],
    )
    trace_rel = px.bar(np.array([abs(np.corrcoef(x, y)[0, 1]) for x in X.values.T]))
    trace_red = px.imshow(abs(np.corrcoef(X.values, rowvar=False)))
    fig.add_trace(trace_rel.data[0], row=1, col=1)
    fig.add_trace(trace_red.data[0], row=1, col=2)
    fig.update_layout(width=1200, height=480, title=title)
    fig.show()


show_relevance_redundancy(
    X,
    y,
    None,
    f"full dataset: acc={evaluate_model(clf, X, y, None):.4f}",
)

Вот результат:

Точность модели 0,9082. Характеристики значительно различаются по релевантности. Есть очаги избыточности, где, по-видимому, существует сильная корреляция между функциями.

Чтобы минимизировать целевую функцию, давайте создадим квадратичную модель с ограничениями и отправим ее в квантовый отжиг. Код будет объяснен ниже.

k = 30
# Pearson correlation
correlation_matrix = abs(np.corrcoef(np.hstack((X, y[:, np.newaxis])), rowvar=False))
# fix the alpha parameter from the Milne paper
# to account for the numerous quadratic terms that are possible
beta = 0.5
alpha = 2 * beta * (k - 1) / (1 - beta + 2 * beta * (k - 1))
# generate weights for linear and quadratic terms, per Milne algorithm
Rxy = correlation_matrix[:, -1]
Q = correlation_matrix[:-1, :-1] * (1 - alpha)
np.fill_diagonal(Q, -Rxy * alpha)
# create binary quadratic model from the linear and quadratic weights
bqm = dimod.BinaryQuadraticModel(Q, "BINARY")
# create constrained quadratic model
cqm = dimod.ConstrainedQuadraticModel()
# the objective function of the CQM is the same as BQM
cqm.set_objective(bqm)
# constraint: limit the number of features to k
cqm.add_constraint_from_iterable(
    ((i, 1) for i in range(len(cqm.variables))), "==", rhs=k
)
# the sampler that will be used is the hybrid sampler in the DWave cloud
sampler = LeapHybridCQMSampler()
# solve the problem
sampleset = sampler.sample_cqm(cqm)

Там много чего происходит. Давайте объясним каждый шаг.

Мы ограничиваем количество признаков k=30. Это основное ограничение модели.

Мы немного отклоняемся от целевой функции, описанной в статье. Вместо прямого определения альфа мы используем эквивалентный параметр бета, который служит той же цели. И затем мы определяем альфу таким образом, чтобы держать под контролем вклад условий взаимодействия — если количество признаков чрезвычайно велико, это гарантирует, что условия взаимодействия не преобладают над линейными членами в целевой функции.

Мы формируем корреляционную матрицу таким образом, чтобы ее можно было использовать для непосредственного построения BQM. Мы ограничиваем BQM условием, что можно использовать не более k = 30 признаков, и поэтому мы получаем ограниченную квадратичную модель (CQM). Мы отправляем CQM в квантовый отжиг и собираем результаты.

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

Мы еще не закончили. Аппаратное обеспечение D-Wave обычно производит выборку пространства решений и возвращает большое количество потенциально возможных решений. Это связано с тем, как работает аппаратное обеспечение: одно событие отжига выполняется быстро и легко, поэтому имеет смысл запускать отжиг несколько раз, что здесь происходит автоматически. Если сгенерировано достаточно образцов, некоторые из них будут лучшим решением.

Поэтому нам нужно отсортировать решения и выбрать лучшее. «Лучший» означает — у него самое низкое значение целевой функции, которую D-Wave называет «энергией», потому что это действительно физическая энергия квантового процессора.

# postprocess results
feasible = sampleset.filter(lambda s: s.is_feasible)
if feasible:
    best = feasible.first
else:
    assert len(cqm.constraints) == 1
    best = sorted(
        sampleset.data(),
        key=lambda x: (list(cqm.violations(x.sample).values())[0], x.energy),
    )[0]
assert list(best.sample.keys()) == sorted(best.sample.keys())
is_selected = np.array([bool(val) for val in best.sample.values()])
features = np.array([i for i, val in enumerate(is_selected) if val])
best_energy = best.energy

features — это массив с индексами признаков, выбранных квантовым отжигом. Это решение процесса выбора функций. Очевидно, его длина будет k=30.

Давайте измерим точность модели после выбора признаков:

show_relevance_redundancy(
    X,
    y,
    features,
    f"explicit optimization: acc={evaluate_model(clf, X, y, features):.4f}",
)

Точность модели после выбора признаков составляет 0,9381. Мы получили ровно столько функций, сколько хотели. Большинство из них имеют высокую актуальность. И значения корреляции между признаками, как правило, низкие. Модель работает лучше, ее легче интерпретировать, а выбор правильных функций занял не так много времени.

Но процесс долгий и утомительный, если все, что вы хотите, это выбрать некоторые функции. К счастью, сейчас есть гораздо более простой способ.

Простой способ выбора функций

Если у вас установлен плагин D-Wave scikit-learn, все, что вам нужно сделать, это:

X_new = SelectFromQuadraticModel(num_features=30, alpha=0.5).fit_transform(X, y)

Это буквально так. За кулисами код создает квадратичную модель, накладывает на нее ограничения, отправляет ее в квантовый решатель, возвращает результаты, анализирует их и выбирает лучший результат для возврата в формате массива NumPy, как это нравится scikit-learn. . Время работы примерно одинаковое. Но давайте посмотрим, как выглядят результаты:

X_new_df = pd.DataFrame(data=X_new, columns=list(range(X_new.shape[1])))

show_relevance_redundancy(
    X_new_df,
    y,
    None,
    f"plugin optimization: acc={evaluate_model(clf, X_new_df, y, None):.4f}",
)

Точность модели составляет 0,9369, что практически совпадает с производительностью, которую мы получили при явной оптимизации (в пределах типичных колебаний различных случайных компонентов).

Набор выбранных функций немного отличается. Это можно объяснить небольшими различиями между моим ручным процессом и автоматизированной реализацией в библиотеке.

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

Дополнительные темы для изучения

Альфа-параметр будет смещать алгоритм выбора признаков либо в сторону меньшей избыточности, но потенциально худшего качества (альфа=0), либо в сторону более высокого качества, но за счет повышенной избыточности (альфа=1). Какое значение является лучшим, зависит от проблемы, которую вы пытаетесь решить.

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

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