Monday, November 28, 2011

Регуляризация линейной регрессии с примерами на R

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

Дабы сильно не усложнять задачу, сегодня я воспользуюсь синтетическими данными. Как обычно, пример будет на языке программирования R. Итак, данные:

"gen" <- function(x){
    return(3*x^2 + 2*x + 2 + rnorm(length(x))*0.5)
}

#Генерируем тренировочные данные
X <- runif(6)
Y <- gen(X)
#Данные для кросс-валидации
Xcv <- runif(50)
Ycv <- gen(Xcv)
Тут все просто, данные распределены согласно некой выдуманной формуле 3*x2+2*x+2 плюс некоторая погрешность, случайным образом распределенная. Тренировочных данных было умышленно взято очень мало, для наглядности (в реальных задачах их все равно, как правило, не хватает, так что это ничуть не мешает обобщению результатов).

Допустим, мы не знаем как на самом деле распределены наши данные, но мы хотим построить какую-нибудь модель. Простейшая линейная регрессия по переменной x нам не подойдет, т.к. линейная регрессия предполагает только линейные зависимости. Поэтому, мы попробуем построить полиномиальную модель. И для этого добавим в наши данные все возможные степени x от 1 до 5.

#Линейная регрессия по полиномиальному набору данных
train <- data.frame(Y, X, X^2, X^3, X^4, X^5)
colnames(train) <- c('Y', 'X', 'X2', 'X3', 'X4', 'X5')
simple <- lm(Y ~ X+X2+X3+X4+X5, train)
error <- sum((predict(simple, train)-Y)^2)/length(Y)
cat("Train error: ",error,"\n")

cv <- data.frame(Xcv, Xcv^2, Xcv^3, Xcv^4, Xcv^5)
colnames(cv) <- c('X', 'X2', 'X3', 'X4', 'X5')
error <- sum((predict(simple, cv)-Ycv)^2)/length(Ycv)
cat("Cross-validation error: ",error,"\n")

#Построим кривую, описывающую наше решение
x <- (1:100)/100
test = data.frame(x, x^2, x^3, x^4, x^5)
names(test) <- c('X', 'X2', 'X3', 'X4', 'X5')
y0 <- predict(simple, test)

plot(X, Y, ylim=range(y0), xlim=c(0,1))
lines(x,y0, col='red')
В результате, ошибка на тренировочных данных близка к нулю, а на данных для кросс-валидации очень велика.
Train error:  4.951098e-13 
Cross-validation error:  6.751967 
Это, так называемый, overfitting. И он абсолютно ожидаемый, ведь через 6 точек на плоскости всегда можно провести полином 5-ой степени.

Тут, читатель может заметить, а к чему, собственно, нам полином 5-ой степени, при таком количестве точек разумней было бы использовать меньше степеней. Согласен, случай выдуманный, но аналогичная проблема возникает в случае сильно-многомерных данных. Например, при обработке текстов. Где переменная X запросто может быть вектором с 10000 элементов.

Что делать?

Для начала, посмотрим на формулу, которой мы пытаемся описать наши данные:
yt = a0 + a1*x + a2*x + a3*x + a4*x + a5*x
Мы пытаемся подобрать эти коэффициенты ai таким образом, чтобы минимизировать ошибку т.е. минимизировать выражение:
L = Σ(yt-y)2
А что если в это выражение добавить еще какой-нибудь член, чтобы как-то уменьшить величины коэффициентов ai? Например, вместо L можно было бы минимизировать выражение:

L1 = Σ(yt-y)2 + λ*Σ|ai|

или другой вариант

L2 = Σ(yt-y)2 + λ*Σ(ai)2

Что мы здесь сделали? Мы добавили в наше выражение штраф за большие значения ai. И величина этого штрафа пропорциональна величине параметра λ, с помощью которого мы теперь можем настраивать наш алгоритм.

Первый вариант носит название L1-регуляризация (в английской литературе LASSO regression), второй вариант L2-регуляризация или регуляризация Тихонова (в английской литературе ridge regression)

Оба варианта реализованы в R. К сожалению, реализация немного не последовательна в используемых входных параметрах и возвращаемых значениях. Тем не менее, мы можем посмотреть её использование чтобы сравнить результат работы.

L1-регуляризация или Lasso regression

L1 регуляризация реализована в пакете lars.

#Lasso regression
#Требуется установить пакет lars если его нет
#install.packages('lars')
library(lars)

train <- cbind(X, X^2, X^3, X^4, X^5)
colnames(train) <- c('X', 'X2', 'X3', 'X4', 'X5')
lasso <- lars(train, Y, type='lasso')

cv <- cbind(Xcv, Xcv^2, Xcv^3, Xcv^4, Xcv^5)
colnames(cv) <- c('X', 'X2', 'X3', 'X4', 'X5')

x <- (1:100)*(max(X)*1.1)/100
test <- cbind(x, x^2, x^3, x^4, x^5)
colnames(test) <- c('X', 'X2', 'X3', 'X4', 'X5')


stats <- NULL
lambda <- 0
plot(X, Y, ylim=c(-5, 10), pch=20, col='red')
points(Xcv, Ycv, pch=20, col='blue')
ls <- NULL
cs <- NULL
for(i in 1:15){
    Yp <- predict(lasso, train, s=lambda, type='fit', mode='lambda')
    trainError <- sum((Yp$fit-Y)^2)/length(Y)

    Yp <- predict(lasso, cv, s=lambda, type='fit', mode='lambda')
    cvError <- sum((Yp$fit-Ycv)^2)/length(Ycv)

    stats <- rbind(stats, c(lambda, trainError, cvError))
    if (i%%5==1){
        #Построим кривую, описывающую наше решение
        y0 <-  predict(lasso, test, lambda, type='fit', mode='lambda')
        lines(x,y0$fit, col=i)
        ls <- c(ls, paste("lambda=",lambda,sep=''))
        cs <- c(cs, i)
    }

    lambda <- ifelse(lambda==0,0.00000001, lambda*10)
}
legend("topleft", c("Train", "Cross-validation"), pch=20, col=c('red', 'blue'))
legend("bottomright",inset=0.05, legend=ls, pch=20, col=cs, text.col=cs, bty="n")



plot(stats[,2], ylim=range(stats[,2:3]), type='l', col='red', ylab="error", xlab="log(lambda)")
lines(stats[,3], col='blue')

#Оптимальное значение lamda=1
coef(lasso, s=1, mode='lambda')

На графике показано несколько решений, для разных значений λ. При λ=0 получается обычная линейная регрессия, которую я описал выше. При очень больших значениях λ получается горизонтальная линия, т.е. решение при котором коэффициенты ai=0, для всех i>0. Оба крайних случая очевидно не являются оптимальными решениями. Как же найти оптимальное? Для этого построим график ошибки.

Здесь, на графике, красным цветом отмечена ошибка для тренировочного множества, а синим цветом ошибка для cross-validation множества. Видно, что ошибка тренировочного множества сначала очень мала (overfitting), но постепенно возрастает (underfitting). В то же время ошибка cross-validation имеет выраженный минимум, который для нас и будет оптимальным значением λ В нашем случае, это оптимальное значение оказалось равным 10.

У L1 регуляризации есть одна полезная особенность. Многие коэффициенты оказываются в точности равны 0. Например, в нашем случае:
X        X2        X3        X4        X5 
0.2989565 0.0000000 0.0000000 0.0000000 5.2679886 
Т.е. L1 регуляризация может быть использована еще и как алгоритм feature selection.

Регуляризация Тихонова или Ridge regression

Этот алгоритм реализован в пакете MASS. К сожалению, реализация не совсем стандартная, например, пакет MASS не предоставляет нам функции predict которой так привычно пользоваться тем, кто знаком с R. Тем не менее, не представляет никакой сложности самостоятельная реализация этой функции через умножение матриц.

#Ridge regression
#Требуется установить пакет MASS если его нет
#install.packages('MASS')
library(MASS)

train <- data.frame(Y, X, X^2, X^3, X^4, X^5)
colnames(train) <- c('Y', 'X', 'X2', 'X3', 'X4', 'X5')

stats <- NULL
lambda <- 0
x <- (1:100)/100
test <- as.matrix(cbind(1, x, x^2, x^3, x^4, x^5))

plot(X, Y, ylim=c(-5,15), pch=20, col='red')
points(Xcv, Ycv, pch=20, col='blue')
ls <- NULL
cs <- NULL
for(i in (1:15)) {

    ridge <- lm.ridge(Y ~ X+X2+X3+X4+X5, train, lambda=lambda)
    Yp <- as.matrix(cbind(1, X, X^2, X^3, X^4, X^5)) %*% as.matrix(coef(ridge))
    trainError <- sum((Yp - Y)^2)/length(Y)

    cv <- as.matrix(cbind(1, Xcv, Xcv^2, Xcv^3, Xcv^4, Xcv^5))
    Yp <- cv %*% as.matrix(coef(ridge))
    cvError <- sum((Yp-Ycv)^2)/lentgh(Ycv)

    stats <- rbind(stats, c(lambda, trainError, cvError))
    if (i%%5==1){
        #Построим кривую, описывающую наше решение
        y0 <- test %*% as.matrix(coef(ridge))
        lines(x,y0, col=i)
        ls <- c(ls, paste("lambda=",lambda,sep=''))
        cs <- c(cs, i)
    }
    lambda <- ifelse(lambda==0,0.00000001, lambda*10)

}
legend("topleft", c("Train", "Cross-validation"), pch=20, col=c('red', 'blue'))
legend("bottomright",inset=0.05, legend=ls, pch=20, col=cs, text.col=cs, bty="n")


#Рисуем изменение ошибки CV
plot(stats[,2], ylim=range(stats[,2:3]), type='l', col='red', ylab="error", xlab="log(lambda)")
lines(stats[,3], col='blue')

#Оптимальное значение lamda=10
ridge <- lm.ridge(Y ~ X+X2+X3+X4+X5, train, lambda=10)
coef(ridge)

Аналогично случаю L1 регрессии, здесь мы тоже можем определить оптимальные значения параметра λ

Оптимальное значение λ в данном случае оказалось равным 1. И соответствующие коэффициенты, как видим, совсем не равны 0.
X        X2        X3        X4        X5 
3.0257674 0.7951164 0.7637522 0.8631815 1.0352908 1.2643954 

Естественно методы регуляризации не ограничиваются L1 и L2-регуляризациями. Существует еще несколько популярных подходов, например алгоритм elastic net, включающий в себя как L1 так и L2.

10 comments:

  1. Спасибо за пост!

    Вопрос: подойдет ли lars для feature selection в SVM with RBF kernel с 25 признаками, или его линейность станет помехой? Про взять и попробовать не так просто, т.к. вся реализация на питоне, а не на R.

    P.S. Полагаю, в AI-class все на 100%? :)

    ReplyDelete
  2. А реализация SVM на питоне самописная или библиотека какая-то? Если это, например, liblinear то она, насколько мне известно, поддерживает L1-регуляризацию. Так что можно было бы действительно взять и попробовать.

    Вообще, L1 SVM активно используются для feature selection, но подойдет ли оно в данном конкретном случае сказать не берусь. Можно поискать дополнительной информации, например тут http://olivier.chapelle.cc/pub/jsm_fsel.pdf

    > P.S. Полагаю, в AI-class все на 100%? :)
    Что-то близкое к этому, но пока не всё так радужно. Когда пришлют сертификат (или что там у них будет) выставлю на обозрение.

    ReplyDelete
  3. LIBSVM (для liblinear маловато признаков, вроде бы). Регуляризация через 'С' есть и там, но она, насколько я понимаю, не влияет на пропорции признаков.
    Под L1-SVM имеется в виду SVM с линейным kernel? Если да, то остается в силе предыдущий вопрос - правильно ли использовать feature selection с помощью линейных классификаторов на 25 признаках, среди которых предполагаются многочисленные нелинейные закономерности?

    > Что-то близкое к этому, но пока не всё так радужно. Когда пришлют сертификат (или что там у них будет) выставлю на обозрение.
    У меня тоже проскальзывают дурацкие ошибки, в результате половина работ на 90-95%. Зато midterm на 100. :)

    ReplyDelete
  4. Под L1-SVM я имел в виду L1-регуляризацию в SVM. LIBSVM вроде бы её не поддерживает, а жаль.

    Если предполагаются нелинейные зависимости, то использовать линейный классификатор для feature selection, конечно, не стоит. С другой стороны, если признаков всего 25 то почему бы не воспользоваться обертками (wrapper) вокруг имеющейся реализации SVM?

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

    ReplyDelete
  5. Сергей,

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

    Заранее спасибо!

    ReplyDelete
  6. Илья,

    Спасибо за замечание. Там была небольшая ошибка (data вместо x) но дело еще в том, что функции runif и rnorm это генераторы случайных чисел. Поэтому в точности мой результат повторить запуская код повторно не получится. Но вы можете заменить строчки:
    X <- runif(6)
    Y <- gen(X)

    на
    X <- c(0.1, 0.35, 0.42, 0.5, 0.82, 0.88)
    Y <- c(2.5, 3, 4, 2.5, 5, 7.5)

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

    ReplyDelete
  7. Сергей,

    Спасибо большое за помощь!
    Странно, что на несуществующую переменную R не выругался, мне странным сначала показалось, что за data там, подумал какая-нить константа :)

    ReplyDelete
  8. А как на счёт более стандратных эконометрических критериев?
    1) t-статистика о значимости данного коэффициента
    2) информационные критерии, котоыре также дают штраф за большое количество переменных

    ReplyDelete
  9. Сколько текстов на эту тему просмотрено, а влияние постоянного сдвига на лямбда и коэффициенты не попалось в обсуждении. А ведь (утрирую), проводя зависимость года рождения от номера ряда в аудитории, результат (лямбда, параметры) будет сильно зависеть от того, будете ли записывать этот год, как напр. 1991, 91, 1 или 0. Частичная ортогонализация базисных функций (хотя бы вычесть из x среднее, и из y тоже) уже чуть улучшила бы ситуацию (хотя бы в смысле однозначности решения). Но и в этому случае, напр. коэффициенты полинома будут зависеть от того, измеряем мы время в годах, месяцах, сутках или секундах. Ну бессмысленно складывать первые или вторые степени от коэффициентов (метр, метр/сек, метр/сек^2 и т.д.). Так что в качестве регуляризирующего параметра уже давно применяются однородные данные, характеризующие гладкость, напр. (y' )^2, (y")^2 и даже (y''')^2. Тогда нет зависимости ни от нуль-пункта x,y, ни от единиц измерения x

    ReplyDelete
  10. Спасибо за статью!
    Вы спасли несколько нерадивых студентов =)

    ReplyDelete