Эффективная онлайн линейная регрессия

53

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

Я предполагаю простую модель линейной многомерной регрессии, т.е.

y=Ax+b+e

Каков наилучший алгоритм для создания постоянно обновляемой оценки параметров линейной регрессии и ?бAb

Идеально:

  • Мне бы хотелось, чтобы алгоритм с наибольшим пространством и временем на обновление, где - размерность независимой переменной ( ), а - размерность зависимой переменной ( ).N x M yO(NM)NxMy
  • Я хотел бы иметь возможность указать какой-либо параметр, чтобы определить, насколько параметры обновляются каждым новым образцом, например, 0,000001 будет означать, что следующий образец даст одну миллионную часть оценки параметра. Это дало бы некоторый экспоненциальный спад влияния образцов в далеком прошлом.
mikera
источник
2
Посмотрите (1) Гибкая линейная регрессия, (2) Фильтры Калмана.
Джейс

Ответы:

31

Майндональд описывает последовательный метод, основанный на вращениях Гивенса . (Вращение Гивенса - это ортогональное преобразование двух векторов, которое обнуляет данную запись в одном из векторов.) На предыдущем шаге вы разложили матрицу проекта в треугольную матрицу через ортогональное преобразование так что . (Получить результаты регрессии быстро и просто из треугольной матрицы.) Присоединив новую строку ниже , вы эффективно расширяете ненулевой строкой тоже скажиT Q Q X = ( T , 0 ) v X ( T , 0 ) t T T t T QXTQQX=(T,0)vX(T,0)t, Задача состоит в том, чтобы обнулить эту строку, сохраняя записи в позиции диагонали . Последовательность вращений Гивенса делает это: вращение с первой строкой обнуляет первый элемент ; затем вращение со второй строкой обнуляет второй элемент и так далее. Эффект заключается в предварительном умножении серией вращений, что не меняет ее ортогональности.TTtTQ

Когда матрица проектирования имеет столбцы (что имеет место при регрессии на переменных плюс константа), необходимое количество вращений не превышает и каждый поворот меняет два вектора. Память, необходимая для составляет . Таким образом, этот алгоритм имеет вычислительную стоимость как во времени, так и в пространстве.p p + 1 p + 1 T O ( ( p + 1 ) 2 ) O ( ( p + 1 ) 2 )p+1pp+1p+1TO((p+1)2)O((p+1)2)

Подобный подход позволяет определить влияние регрессии на удаление строки. Майндональд дает формулы; как и Белсли, Кух и Уэльс . Таким образом, если вы ищете движущееся окно для регрессии, вы можете сохранить данные для окна в круговом буфере, примыкая к новому элементу данных и отбрасывая старый при каждом обновлении. Это удваивает время обновления и требует дополнительной памяти для окна шириной . Похоже, что будет аналогом параметра влияния.k 1 / kO(k(p+1))k1/k

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

Рекомендации

JH Maindonald, Статистические вычисления. J. Wiley & Sons, 1984. Глава 4.

Д. А. Белсли, Е. Кух, Р. Уэлш. Регрессионная диагностика: идентификация влиятельных данных и источников коллинеарности. J. Wiley & Sons, 1980.

Whuber
источник
1
Является ли метод, описанный Майндональдом, связанным или таким же, как алгоритм Джентльмена? jstor.org/stable/2347147
онстоп
6
В этом случае см. Также расширения Алана Миллера jstor.org/stable/2347583 . Архив его программного сайта на Фортране
onestop
5
Явный алгоритм появляется в нижней части р. 4 из saba.kntu.ac.ir/eecd/people/aliyari/NN%20%20files/rls.pdf . Это может быть найдено Гуглингом "рекурсивные наименьшие квадраты". Это не похоже на улучшение подхода Джентльмена / Майндональда, но, по крайней мере, оно четко и явно описано.
whuber
2
Последняя ссылка выглядит как метод, который я собирался предложить. Используемая ими матричная идентичность известна в других местах как идентичность Шермана - Моррисона - Вудбери. Это также достаточно эффективно для реализации, но может быть не таким стабильным, как вращение Гивенса.
кардинал
2
@suncoolsu Хм ... Книга Майндональда была недавно опубликована, когда я начал ее использовать :-).
whuber
8

Я думаю, что преобразование вашей модели линейной регрессии в модель пространства состояний даст вам то, что вам нужно. Если вы используете R, вы можете использовать пакет dlm и взглянуть на книгу-компаньон Petris et al.

Ф. Туселл
источник
может быть, я запутался, но, похоже, это относится к модели временного ряда? моя модель на самом деле проще в том смысле, что выборки не являются временными рядами (фактически они являются независимыми (x-> y) выборками, они просто накапливаются в больших объемах с течением времени)
mikera
1
Да, в общем случае это используется для временных рядов с независимыми наблюдениями; но вы всегда можете предположить несоответствие между последовательными наблюдениями, что представляет особый интерес для вас.
Ф. Туселл
7

Вы всегда можете просто выполнить градиентный спуск по сумме квадратов стоят WRT параметров вашей модели . Просто возьмите его градиент, но не используйте решение для закрытой формы, а только для направления поиска.EW

Пусть является стоимость i - й обучающей выборки , учитывая параметры . Ваше обновление для параметра j'th тогдаE(i;W)W

WjWj+αE(i;W)Wj

где - это скорость шага, которую вы должны выбрать с помощью перекрестной проверки или надлежащей меры.α

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

Конечно, могут быть применены более сложные алгоритмы оптимизации (импульс, сопряженный градиент, ...).

bayerj
источник
Кажется, очень похож на этот документ eprints.pascal-network.org/archive/00002147/01/… . Это было реализовано в проекте с открытым исходным кодом под названием jubatus.
сахарин
3

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

L(β)=(yXβ)t(yXβ)
Градиент становится и гессианом:
L(β)=2Xt(yXβ)
2L(β)=XtX

Теперь вы получили некоторые прошлые данные, выполнили линейную регрессию и работаете с вашими параметрами ( ). Градиент в этой точке равен нулю по определению. Гессиан, как указано выше. Новая точка данных ( ) прибывает. Вы просто рассчитываете градиент для новой точки через:βxnew,ynew

Lnew(β)=2xnew(ynewxnewTβ)
и это станет вашим общим градиентом (так как градиент от существующих данных был нулевым) , Гессен для новой точки данных:

2Lnew=xnewxnewT
.

Добавьте это к старому гессиану, указанному выше. Затем просто сделайте шаг Ньютона Рафсона.

βnew=βold+(2L)1Lnew

И вы сделали.

ryu576
источник
1
Мне нравится идея из-за ее простоты, но (а) ради того, чтобы не сбивать с толку читателей, я предпочел бы увидеть явное определение " " и (б) считаю, что вам нужно правильно вычислить градиент (или показать почему быть в 2 раза не имеет значения). Это было бы более убедительно, если бы вы могли привести небольшой пример, демонстрирующий, что это правильно. Для больших было бы целесообразно оценить вычислительные усилия. Разве инверсия гессиана не занимает времени? Lnewp,O(p3)
whuber
Спасибо, добавлю больше деталей чуть позже сегодня. Да, инвертирование гессиана требует для больших . Вы также можете попытаться сохранить обратное гессиану и обновить его напрямую, используя степенные ряды ( ). Если у вас есть миллион параметров, градиентный спуск в любом случае является более или менее единственным вариантом. O(p3)p(IA)1=I+A+A2+
ryu576
2

Стандартный метод наименьших квадратов дает коэффициенты регрессии

β=(XTX)1XTY

где X - это матрица из M значений для каждой из N точек данных, и имеет размер NXM. Y - матрица выходов NX1. конечно матрица коэффициентов MX1. (Если вы хотите перехват, просто сделайте один набор х равным всегда 1.)β

Предположительно, чтобы сделать это онлайн, вам просто нужно отслеживать и , то есть одну матрицу MXM и одну матрицу MX1. Каждый раз, когда вы получаете новую точку данных, вы обновляете эти элементы , а затем снова вычисляете , что стоит вам инверсии матрицы MXM и умножения матрицы MXM и матрицы MX1.XTXXTYM2+Mβ

Например, если M = 1, то один коэффициент

β=i=1Nxiyii=1Nxi2

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

Если вы хотите геометрически ослабить более ранние оценки, я полагаю, что вы могли бы взвешивать и на каждый раз перед добавлением нового члена, где - это небольшое число.X T Y ( 1 - λ ) λXTXXTY(1λ)λ

Марк Хиггинс
источник
2
Приятно видеть объяснение этого простого случая. Вы заметили, однако, что вопрос конкретно задает о многомерной регрессии? В таком случае обновить знаменатель не так просто ! β
whuber
Я думаю, что мой ответ работает до сих пор: т.е. вам нужно отслеживать матрицу MxM и матрицу Mx1 . Каждый элемент этих матриц является суммой, как в примере M = 1. Или я что-то упустил? X T YXTXXTY
Марк Хиггинс
6
Да: в дополнение к вычислению матричного произведения и применению матрицы к вектору, теперь вам нужно инвертировать на каждом шаге. Это дорого. Весь смысл алгоритмов онлайн состоит в том, чтобы заменить оптовые дорогие шаги более дешевыми процедурами обновления. XX
whuber
1
@whuber Кстати, быстрый онлайновый способ оценки для изменяющейся матрицы и вектора дан Schraudolph, NN (2002). Матрично-векторные произведения быстрой кривизны для градиентного спуска второго порядка. По сути, вы позволяете и как . C x z t + 1 = z t + x - C z t z C - 1 x t C1xCxzt+1=zt+xCztzC1xt
Нил Г
1

Проблему легче решить, если немного переписать вещи:

Y = y

X = [x, 1]

тогда

Y = A * X

Одноразовое решение найдено путем расчета

V = X '* X

а также

C = X '* Y

обратите внимание, что V должен иметь размер N-by-N, а C - размер N-by-M. Параметры, которые вы ищете, затем определяются как:

A = inv (V) * C

Поскольку и V, и C рассчитываются путем суммирования по вашим данным, вы можете рассчитать A для каждой новой выборки. Однако это имеет временную сложность O (N ^ 3).

Поскольку V является квадратным и полуопределенным положительным, LU-разложение существует, что делает инвертирование V численно более устойчивым. Существуют алгоритмы для обновления ранга 1 для обратной матрицы. Найдите их, и вы получите эффективную реализацию, которую вы ищете.

Алгоритмы обновления ранга 1 можно найти в «Матричных вычислениях» Голуба и Ван Лоан. Это сложный материал, но он имеет исчерпывающий обзор таких алгоритмов.

Примечание: метод выше дает оценку наименьших квадратов на каждом шаге. Вы можете легко добавлять веса к обновлениям X и Y. Когда значения X и Y становятся слишком большими, они могут масштабироваться одним скаляром, не влияя на результат.

Мистер Уайт
источник