[Из песочницы] MLBootCamp «Оценка производительности». Очень простой и быстрый вариант решения
1. Загрузка и подготовка данных
Коротко о соревновании, вот как поставлена задача:
«Мы предлагаем вам научить компьютер предсказывать, сколько секунд будут умножаться две матрицы размера mｘk и kｘn на данной вычислительной системе, если известно, сколько решалась эта задача на других вычислительных системах с другими размерами матриц».
Стоит сразу отметить, что строго говоря это не так. Нет ни одной вычислительной системы в тестовой выборке, которой бы не было в тренировочной. В качестве метрики используется средняя процентная ошибка (MAPE).
А теперь немного подробнее. В обеих выборках для каждого наблюдения имеется 951 признак. Первые три из них это размерности умножаемых матриц, остальные — параметры систем, на которых проводились измерения. Оказывается, что если отбросить параметры умножаемых матриц, то в тренировочной выборке 92 уникальных конфигурации использовавшихся вычислительных систем, т.е. оставшиеся 948 параметров имеют всего лишь 92 комбинации. Заменив эти 948 параметров одним, идентификатором конфигурации, мы получим, что в наших данных всего лишь 4 признака для каждого наблюдения. Подготовим данные для дальнейшей работы добавив еще произведение всех размеров матриц. Теперь в данных есть признаки, m, k, n, m*k*n и conf (идентификатор компьютера). Также, тренировочная выборка содержит целевую переменную time — время, за которое было выполнено умножение матриц.
library(quantreg)
library(dplyr)
# загрузка исходных данных
y <- read.csv('y_train.csv')[,1]
X.train <- read.csv('x_train.csv', stringsAsFactors = FALSE)
X.test <- read.csv('x_test.csv', stringsAsFactors = FALSE)
# создадим таблицу с уникальными конфигурациями
X.all <- rbind(X.train, X.test)
X.uni <- unique(X.all[,4:951])
# создадим датафрейм с уникальными конфигурациями и вектор с их именами
X.uni <- cbind(conf=paste0('conf', formatC(1:nrow(X.uni), 2, flag = '0')), X.uni)
confs <- unique(X.uni$conf)
# разметим тренировочную выборку именами и выкинем все ненужные столбцы
# преобразуем целый тип к вещественным и добавим целевую переменную
cols <- c('m', 'k', 'n', 'conf')
X.train <- dplyr::inner_join(X.train, X.uni, by=colnames(X.uni)[2:ncol(X.uni)])[cols]
X.train[,1:3] <- lapply(X.train[,1:3], as.numeric) # int->numeric
X.train$mkn <- X.train$m * X.train$k * X.train$n
X.train$time <- y
# то же для тестовой выборки
X.test <- dplyr::inner_join(X.test, X.uni, by=colnames(X.uni)[2:ncol(X.uni)])[cols]
X.test[,1:3] <- lapply(X.test[,1:3], as.numeric) # int->numeric
X.test$mkn <- X.test$m * X.test$k * X.test$n
2. Предварительный анализ данных и выбор рабочей модели
Теперь можно посмотреть на данные и увидеть ожидаемую зависимость целевой переменной time от m*k*n. Например, конфигурация 39 (слева). Она хорошо аппроксимируется линейной регрессией. Однако, не для всех конфигураций линейная модель может дать хороший результат, например, конфигурация 11 (справа):
Шум и выбросы вносят свой весомый отрицательный вклад модель, поэтому один из вариантов борьбы с ними — удалять такие наблюдения. Недостаток этого способа это решить, насколько наблюдение должно быть «неправильным», чтобы его можно было удалить. Причем для некоторых систем очень мало данных и выкидывать даже зашумленные данные расточительно.
Альтернативный вариант — пользоваться моделями, которые устойчивы к выбросам и шуму. Примеры таких моделей это квантильная и робастная регрессии.
В R их можно найти в пакетах quantreg и MASS соответственно. Каждая из них имеет свои настройки, и/или более продвинутые вариации. Здесь эти методы использовались с параметрами по умолчанию. Вот как будет выглядеть график со всеми регрессиями для шумной конфигурации 11:
3. Обучение и тест моделей
Разброс значений целевой переменной варьируется в широких пределах для каждой из систем, поэтому напрашивается решение: обучить стек моделей, по одной на каждую конфигурацию. Пусть зависимость в каждой модели будет немного сложнее, а именно время затраченное от умножение матриц будет зависеть не только от произведения их размеров, но и от каждого из них в отдельности, т.е. time ~ m + k + n + mkn:
models <- list()
for (i in 1:length(confs)){
# выбор конфигурации
X.conf <- subset(X.train, conf==confs[i])
X.conf$conf <- NULL
# обучаем модель и добавляем ее в стек
fit <- rq(time ~ ., data=X.conf)
models[[confs[i]]] <- fit
}
Теперь легко рассчитать предсказания перебирая модели и выбирая из датасета данные, соответствующие каждой отдельно взятой конфигурации и ее модели.
y.test <- rep(1, nrow(X.test))
y.train <- rep(1, nrow(X.train))
# предсказания на тренировочной выборке
for (i in 1:length(confs)){
X.conf <- subset(X.train, conf==confs[i])
y.train[as.numeric(rownames(X.conf))] <- predict(models[[confs[i]]])
}
# предсказания на тестовой выборке
for (i in 1:length(confs)){
X.conf <- subset(X.test, conf==confs[i])
y.test[as.numeric(rownames(X.conf))] <- predict(models[[confs[i]]], newdata = X.conf)
}
В силу неидеальности исходных данных и, следовательно, полученных моделей в предсказанных значениях будут числа меньше 1, которые по условию задачи быть не должны. Здесь оказывается полезным воспользоваться предсказаниями на тренировочной выборке и заменить все неподходящие значения. Простым и эффективным методом является замена всех чисел, меньших определенной константы, отсечки, ей самой. На графике представлен график метрики (MAPE) в зависимости от значения отсечки cutoff. Горизонтальная линия — метрика без замены, вертикальная — отсечка, при котором метрика принимает наименьшее значение.
Используя такую замену можно немного улучшить значение метрики. Как показала практика, уровень отсечки на тренировочной и тестовой выборках приблизительно одинаковый.
Очевидный недостаток приведенного выше кода это то, что модели строились без кросс-валидации, поэтому полученному значению метрики на тренировочной выборке доверять нельзя, как и оценке на публичном лидерборде. Недостаток ликвидируется дополнительным соответствующим кодом.
Однако имеется и явное преимущество — в скорости (у автора скрипт выполнился за полсекунды), что позволяет в короткое время оценить множество различных вариантов комбинаций моделей. Для сравнения, модели с помощью градиентного бустинга, случайных лесов или нейронных сетей строились от как минимум нескольких десятков минут.
Тем не менее, алгоритм с выбором между квантильной и робастной регрессией для каждой конфигурации, показал результат MAPE=5.23% на публичном лидерборде. По итогам соревнования метрики нет, но на основе данных по кросс-валидации ожидается MAPE около 5.4%.
4. Заключение
Таким образом, обучение стека моделей и прогноз целевой переменной выполняются очень быстро, меньше секунды и при этом достигается хорошая точность. Описанный способ можно улучшить за счет более детального внимания к отдельно взятой конфигурации, например, подбирать квантиль для каждой регрессии или параметры при использовании робастных методов, удалять явные выбросы, подбирать отдельно отсечку и т.д., что, как правило, не сильно сказывается на времени выполнения скрипта.
По скромному мнению автора, использование квантильной (или другой, устойчивой к выбросам) регрессии для предсказания времени выполнения компьютерных экспериментов является весьма практичным, точным и куда менее вычислительно затратным, чем способ приведенный в статье к соревнованию, где использовались линейная регрессия с удалением значений (собственный велосипед?) и случайные леса.
