Alink's Talk (11): L-BFGS Оптимизация линейной регрессии

0x00 сводка

Alink — это платформа алгоритмов машинного обучения нового поколения, разработанная Alibaba на основе вычислительного движка реального времени Flink.Это первая в отрасли платформа машинного обучения, которая поддерживает как пакетные, так и потоковые алгоритмы. В этой статье рассказывается, как L-BFGS-оптимизация линейной регрессии реализована в Alink, и я надеюсь, что ее можно будет использовать в качестве дорожной карты для всех, чтобы увидеть код линейной регрессии.

Поскольку общедоступная информация Alink слишком мала, нижеследующее - все предположения, и обязательно будут упущения и ошибки. Я надеюсь, что все укажут, и я обновлю их в любое время.

0x01 Обзор

До сих пор ввод был обработан, см.Alink's Talk (10): Предварительная обработка данных для реализации линейной регрессиис последующей оптимизацией. Основная цель оптимизации - найти направление, после движения параметров в этом направлении значение функции потерь может быть уменьшено Это направление часто получается различными комбинациями частных производных первого порядка или частных производных второго порядка. Итак, давайте еще раз рассмотрим основную идею.

1.1 Основная идея оптимизации

Для модели линейной регрессии f(x) = w'x+e мы строим функцию стоимости (функцию потерь)J(θ), и найдяJ(θ)Минимальное значение функции, коэффициент w линейной модели можно определить.

Окончательная функция оптимизации: min(L(Y, f(x)) + J(x)), которая предназначена для оптимизации эмпирического риска и структурного риска, и эта функция называетсяцелевая функция.

Что нам нужно сделать, так это выбрать оптимальное θ на основе нашего тренировочного набора и сделать f(x) как можно более близким к истинному значению в нашем тренировочном наборе. Мы определяем функцию для описания «разрыва между f (x) и истинным значением y», эта функция называется целевой функцией, и выражение выглядит следующим образом:

J(θ)≈\frac{1}{2} \sum_{i≈1}^m(f_θ(x^{(i)} - y^{(i)})^2\\

Целевой функцией здесь является знаменитаяфункция наименьших квадратов.

Мы хотим выбрать оптимальное θ такое, чтобы f(x) было ближе всего к истинному значению. Эта задача трансформируется в решение оптимального θ, при котором целевая функция J(θ) принимает минимальное значение.

1.2 Различные методы оптимизации

найти правильныйWПусть целевая функцияf(W)Минимум — это задача оптимизации без ограничений.Общий способ решения этой задачи — задать случайным образом начальноеW0, итерируя, вычисляя направление убывания целевой функции на каждой итерации и обновляяW, пока целевая функция не стабилизируется в точке минимума.

Отличие разных алгоритмов оптимизации заключается в направлении убывания целевой функции.Dtрасчет. Нисходящее направление определяется целевой функцией в текущемWЗатем найдите обратную величину первого порядка (Градиент, Градиент) и найдите производную второго порядка (Матрица Гессе, Матрица Гессе). Общие алгоритмы включают градиентный спуск, метод Ньютона и метод квазиньютона.

Метод градиентного спуска напрямую использует целевую функцию в текущейWВ качестве нисходящего направления принимается противоположное направление градиента.
В настоящее время метод НьютонаWЗатем используйте квадратичное расширение Тейлора, чтобы аппроксимировать целевую функцию, а затем используйте аппроксимированную функцию, чтобы решить нисходящее направление целевой функции. вBtявляется целевой функциейf(W)существуетWtМатрица Гессе при . Это направление поиска также называют направлением Ньютона.
Квазиньютоновский метод требует только вычисления градиента целевой функции на каждой итерации и нахождения приблизительной матрицы Гессе для вычисления ньютоновского направления путем подгонки.
L-BFGS (BFGS с ограниченной памятью) устраняет необходимость сохранения после каждой итерации в BFGS.N*NДля задачи обратной матрицы Гессе вам нужно сохранить только два набора векторов и один набор скаляров для каждой итерации.

В Alink UnaryLossObjFunc — это целевая функция, SquareLossFunc — это функция потерь, а для оптимизации модели используется алгоритм L-BFGS..

То есть направление оптимизации задаетсяКвазиньютоновский метод L-BFGSСделайте это (как это сделать, чтобы увидеть разложение Тейлора второго порядка для f (x)), минимальное значение функции потерь равноквадрат функции потерьвычислять.

0x02 Основные понятия

Поскольку алгоритм L-BFGS является разновидностью квазиньютоновского метода, мы начнем с расширения Тейлора сути метода Ньютона.

2.1 Расширение Тейлора

Расширение Тейлора заключается в использовании набора простых степенных функций для аппроксимации локального значения сложной функции f(x) в этом интервале на основе расширения точки x_0 в определенном интервале. Сценарий применения разложения Тейлора, например: трудно напрямую получить значение sin(1), но мы можем получить приблизительное значение sin(1) через ряд Тейлора sin, и чем больше членов разложения, тем выше точность. Компьютеры обычно вычисляют sin(x) по разложению Тейлора.

Почему Тейлор изобрел эту формулу?

Потому что в то время исследования и применение простых функций в математическом сообществе созрели, а сложные функции, такие как: f(x) = sin(x)ln(1+x), на первый взгляд были головной болью, и были и кривые выражения просто не нашел. Трудно сделать что-либо, кроме как подставить x, чтобы получить его y. Итак, ученики Тейлора приняли вызов! Решил позволить этим формулам показать свою первоначальную форму и упростить их.

Чтобы сделать сложную функцию простой, можно ли преобразовать ее в другое выражение? Формула Тейлора описывается одним предложением: она использует полиномиальную функцию для аппроксимации гладкой функции. То есть по математической идее «подставить прямо музыку и разбить целое на ноль» была выведена формула Тейлора.

Формула ТейлораПутем преобразования (перезаписи) [выражения произвольной функции] в [полиномиальную] форму, это чрезвычайно мощная функцияИнструмент аппроксимации. Почему говорят, что он мощный?

Многочлены очень [дружественные], их легко вычислить, легко вывести, легко интегрировать.
Чувство геометрии и расчета очень интуитивное, типа парабола и несколько квадратов, умножь основание на себя и умножь на себя.

Как рассуждать простым языком?

Что делает формула Тейлора:Оцените (приблизительно) значения вокруг, используя полиномиальное выражение.

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

Физики пришли к выводу: Применяя жизненный опыт «имитации» к кинематическим задачам, если вы хотите имитировать кривую, вы должны, во-первых, убедиться, что начальная точка кривой одинакова, а во-вторых, обеспечить скорость изменения перемещение во времени в начальной точке То же самое (одинаковая скорость), еще раз следует убедиться, что первые два равны, а скорость изменения второго порядка по времени одинакова (одинаковое ускорение). , Если каждая скорость изменения (каждая производная) во времени одинакова, то эти две кривые определенно эквивалентны.

Итак, если Тейлор нашел способ не прикасаться к таким функциям, как sin(x), то естьЗамените такие функции.По изображению таких функций можно имитировать изображение, подобное исходному, такое поведение в математике называется аппроксимацией. Не говорите об этом термине. Расскажите о том, как подделать изображение.

Первый шаг в подделке состоит в том, чтобы кривая подделки также прошла эту точку.

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

После прохождения второго шага теперь отправная точка та же, и общая тенденция изменений аналогична, что может показаться немного интересным. Для большей точности следует учитывать вогнутость и выпуклость. Средняя школа: Параметр, характеризующий вогнутость и выпуклость изображения, является «производной производной». Итак, следующий шаг — сделать производные двух производных равными.

Начальная точка та же, увеличение и уменьшение одинаковы, вогнутость и выпуклость одинаковы, а имитируемая функция более похожа. Если вы продолжите его уточнять, он должен быть бесконечно близок. Так думает Тейлор».Чтобы имитировать кривую, сначала убедитесь, что начальная точка та же самая, затем убедитесь, что производная здесь такая же, и продолжайте следить за тем, чтобы производная производной здесь была такой же..."

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

То есть естьИсходная функция f(x), я воссоздаю изображение, похожее на исходное функциональное изображениеПолиномиальная функция g(x), чтобы обеспечить сходство, мне нужно только убедиться, что эти две функции находятся в определенной точкеНачальное значение равно, первая производная равна, вторая производная равна и n-я производная равна.

2.2 Метод Ньютона

Основная идея метода Ньютона заключается в следующем:Выполните разложение Тейлора второго порядка для f (x) рядом с существующей минимальной оценкой, чтобы найти следующую оценку минимума.. Основная идея состоит в том, чтобы использовать производную первого порядка (градиент) и производную второго порядка (матрицу Гессена) в точке итерации x_k для приближения квадратичной функции к целевой функции, а затем использовать точку минимума квадратичной модели как новой итерационной точки, и непрерывно Этот процесс повторяется до тех пор, пока не будет получено приближенное минимальное значение, удовлетворяющее точности.

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

То, что мы хотим оптимизировать, — это многомерные функции, а x часто является не действительным числом, а вектором. Следовательно, первая производная f(x) также является вектором, а вторая производная производной является матрицей, которая является матрицей Гессе.

2.2.1 Разложение Тейлора первого порядка

Метод Ньютона для корней можно расширить по первому порядку Тейлора. Например, для уравнения f(x) = 0 мы выполняем разложение Тейлора первого порядка в любой точке x0:

f(x) = f(x_0) + (x - x_0)f^{ '}(x_0)

Пусть f(x) = 0, и подведя его к приведенной выше формуле, мы можем получить:

x = x_0 - \frac{f(x_0)}{f^{'}(x_0)}

Обратите внимание, что здесь используется приближение, и полученный x не является корнем уравнения, а является приблизительным решением. Однако x ближе к корням уравнения, чем x0.

Затем решите его итеративно, где x равно x0, чтобы найти следующую позицию ближе к корню уравнения. Итеративную формулу можно записать в виде:

x_{n+1} = x_n - \frac{f(x_n)}{f^{'}(x_n)}

2.2.2 Разложение Тейлора второго порядка

В машинном обучении и глубоком обучении проблема оптимизации функции потерь обычно основана на градиентном спуске по первой производной. Теперь, с другой точки зрения, попытка минимизировать функцию потерь на самом деле является максимальной задачей, соответствующей первой производной функции f'(x) = 0. То есть, если мы находим точку x, где f'(x) = 0, функция потерь минимизируется, и цель оптимизации модели достигается.

Теперь цель состоит в том, чтобы вычислить x, соответствующий f'(x) = 0, корень f'(x) = 0. Преобразовав в корневую проблему, вы можете использовать метод Ньютона из предыдущего раздела.

В отличие от предыдущего раздела, во-первых, выполните разложение Тейлора второго порядка для f (x) в точке x0:

f(x) = f(x_0) + (x - x_0)f^{'}(x_0) + \frac{1}{2}(x-x_0)^2f^{''}(x_0)

Приведенная выше формула верна, если f(x) приблизительно равно f(x0). Пусть f(x) = f(x0) и возьмем производную от (x - x0), получим:

x = x_0 - \frac{f^{'}(x_0)}{f^{''}(x_0)}

Кроме того, хотя x не является оптимальной точкой решения, x ближе к корню f'(x) = 0, чем x0. Таким образом, можно получить оптимизированную итерационную формулу:

x_{n+1} = x_n - \frac{f^{'}(x_n)}{f^{''}(x_n)}

С помощью итерационной формулы можно непрерывно находить приблизительный корень f'(x) = 0, тем самым достигая цели оптимизации по минимизации функции потерь.

2.2.3 Многомерное пространство

В задаче оптимизации машинного обучения все, что мы хотим оптимизировать, — это многомерные функции, поэтому x часто является не действительным числом, а вектором. Следовательно, когда в машинном обучении используется метод поиска корня Ньютона, x является вектором, производная первого порядка f(x) также является вектором, а производная второго порядка производной является матрицей, которая является Матрица Гессе.

В многомерном пространстве мы заменяем первую производную на градиент, а вторую производную на матрицу Гессе, и итерационная формула метода Ньютона остается неизменной:

x_{k+1} = x_k - [Hf(x_k)^{-1}].J_f(x_k)

где J определяется какМатрица Якоби, соответствующая частной производной первого порядка.

Вывод следующий:

Мы предполагаем, что f(x) — дифференцируемая действительная функция второго порядка, а Тейлор разлагает f(x) в точке xk и принимает приближение второго порядка как

\begin{выровнено}&f(x)≈f(x^k)+∇f(x^k)T(x−x^k)+\frac{1}{2}(x−x^k)^T ∇^2f(x^k)(x−x^k)\\&x^k — текущая минимальная оценка\\&∇f(x^k) — оценка f(x) в точке x^k Производная \ &∇^2f(x) — матрица Гессена функции f(x) в точке x^k. \\\конец{выровнено}

Наша цель — найти минимальное значение f(x), а точка с производной 0, скорее всего, будет крайней точкой, поэтому здесь мы берем производную f(x) и делаем ее производную равной 0, то есть , ∇f(x )=0, мы можем получить

∇f(x)=∇f(x^k)+∇^2f(x^k)(x−x^k)=0

Предполагая обратимость ∇2f(x), итеративную формулу метода Ньютона можно получить из (2)

\begin{выровнено}&x^{k+1}=x^k-∇^2f(x^k)^{-1}∇f(x^k) \\&d= -∇^2f(x^k) ^{−1}∇f(x^k) называется ньютоновским направлением\\\end{aligned}

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

2.2.4 Основной процесс метода Ньютона

Подытожим (моделируем) этапы нахождения корней по методу Ньютона:

a.已知函数的情况下,随机产生点.

b.由已知的点按照的公式进行k次迭代.

c.如果与上一次的迭代结果相同或者相差结果小于一个阈值时,本次结果就是函数的根.

Псевдокод выглядит следующим образом

def newton(feature, label, iterMax, sigma, delta):
    '''牛顿法
    input:  feature(mat):特征
            label(mat):标签
            iterMax(int):最大迭代次数
            sigma(float), delta(float):牛顿法中的参数
    output: w(mat):回归系数
    '''
    n = np.shape(feature)[1]
    w = np.mat(np.zeros((n, 1)))
    it = 0
    while it <= iterMax:
        g = first_derivativ(feature, label, w)  # 一阶导数
        G = second_derivative(feature)  # 二阶导数
        d = -G.I * g
        m = get_min_m(feature, label, sigma, delta, d, w, g)  # 得到步长中最小的值m
        w = w + pow(sigma, m) * d
        it += 1       
    return w

2.2.5 Проблемы и решения

Метод Ньютона требует вычисления не только гессиана, но и обратного гессиана. Когда объем данных относительно невелик, скорость работы не сильно пострадает. Однако, когда количество данных велико, особенно в глубоких нейронных сетях, вычисление матрицы Гессе и ее обратной требует очень много времени. В целом скорость оптимизации метода Ньютона не такая высокая, как у алгоритма градиентного спуска. Поэтому большинство современных стратегий оптимизации функций потерь нейронной сети основаны на градиентном спуске.

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

Когда размерность x очень велика, нам очень трудно найти вторую производную от f(x), а метод Ньютона для нахождения критической точки является итеративным алгоритмом, поэтому нам приходится сталкиваться с этой трудностью бесконечно много раз, что приводит к методу Ньютона.Застой Фацю нельзя использовать в машинном обучении. Есть ли решение?

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

Идея квазиньютоновского метода на самом деле очень проста: как разница между двумя точками значения функции может аппроксимировать производную, так и разница между двумя точками первой производной может также аппроксимировать вторую производную. Геометрический смысл в том, чтобы найти «секанс» первой производной, при переходе к пределу секанс будет приближаться к тангенсу, а матрица В будет приближаться к матрице Гессе.

Метод квазиньютона, как следует из названия, имитирует метод Ньютона, заменяя ∇^2f(x)^{−1} матрицей H_{k+1}, которая аппроксимирует ∇^2f(x)^{−1}

2.3 Квазиньютоновский метод

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

В различных квазиньютоновских методах используются итерационные методы для аппроксимации обратной матрицы Гессе и самой себя соответственно;

В различных квазиньютоновских методах общая стратегия построения Hk+1 заключается в том, что H1 обычно выбирает любую симметричную положительно определенную матрицу n-го порядка (обычно I), а затем дает Hk+1, непрерывно корректируя Hk, т. е.

H_{k+1}=H_k+ΔH_k \\ ΔH_k называется корректирующей матрицей

Например, метод BFGS должен обновлять матрицу H (обратную матрицу матрицы Гессе) каждый раз, когда k-й шаг итерационной разности точек s и градиентной разности y, а k+1-й шаг H эквивалентен необходимости от начала до k-го шага и значение y.

Нам нужно найти корень функции с помощью метода стационарной точки Ньютона и алгоритма BFGS Оба алгоритма требуют итерации, поэтому мы просто позволяем им выполнять итерацию вместе. Оба алгоритма приближаются к корню функции медленно, поэтому после k итераций полученное решение является корнем производной функции целевой функции в машинном обучении. Этот метод расчета двух алгоритмов, повторяющихся вместе, мы называем «на лету».

На первом шаге итерации алгоритма BFGS единичная матрица умножается на градиент g, который равен градиенту g, формально совпадающему с выражением для градиентного спуска. Следовательно, алгоритм BFGS можно понимать как алгоритм, который постепенно переходит от градиентного спуска к методу Ньютона для поиска решений функций.

Хотя мы использовали алгоритм BFGS для постепенной аппроксимации матрицы H единичной матрицей, матрица D должна сохраняться каждый раз при ее вычислении, а также размер матрицы D. Шерстяная ткань. Предполагая, что наш набор данных имеет 100 000 измерений (не слишком много), результат хранения D-матрицы за итерацию составляет 74,5 ГБ. Мы не можем сохранить такой огромный матричный контент, как это решить? Используйте алгоритм L-BFGS.

2.4 Алгоритм L-BFGS

Основная идея алгоритма L-BFGS заключается в том, что алгоритм только сохраняет и использует информацию о кривизне последних m итераций для построения приближенной матрицы матрицы Гессе.

Хотя мы использовали алгоритм BFGS для постепенной аппроксимации матрицы H единичной матрицей, матрица D должна сохраняться каждый раз при ее вычислении, а также размер матрицы D. Шерстяная ткань. Предполагая, что наш набор данных имеет 100 000 измерений (не слишком много), результат хранения D-матрицы за итерацию составляет 74,5 ГБ. Мы не можем сохранить такой огромный матричный контент, как это решить? Используйте алгоритм L-BFGS.

Каждая итерация нашей D-матрицы получается путем итеративного вычисления sk и yk. Когда наша память не может их сохранить, мы можем только выбросить некоторые данные, которые не могут их сохранить. Предполагая, что мы установили количество сохраненных векторов равным 100, когда s и y повторяются более 100, первые s и y отбрасываются. Самые передние s и y отбрасываются соответственно для каждой дополнительной итерации. Таким образом, хотя точность и теряется, можно гарантировать, что решение функции будет получено алгоритмом BFGS с использованием ограниченной памяти. Следовательно, алгоритм L-BFGS можно понимать как другое приближение или приближение к алгоритму BFGS.

Математические рассуждения здесь не вводятся, потому что в Интернете есть много отличных статей, а здесь представлена только инженерная реализация. Общие шаги обобщения алгоритма L-BFGS следующие:

Step1:  选初始点x_0，存储最近迭代次数m；
Step2:  k=0, H_0=I, r=f(x_0)；
Step3:  根据更新的参数计算梯度和损失值，如果达到阈值，则返回最优解x_{k+1}，否则转Step4；
Step4:  计算本次迭代的可行梯度下降方向 p_k=-r_k；
Step5:  计算步长alpha_k，进行一维搜索;
Step6:  更新权重x;
Step7:  只保留最近m次的向量对；
Step8:  计算并保存 s_k, y_k
Step9:  用two-loop recursion算法求r_k；
k=k+1，转Step3。

0x03 Модель оптимизации -- алгоритм L-BFGS

вернуться к коду,BaseLinearModelTrainBatchOp.optimizeВызов функции

return new Lbfgs(objFunc, trainData, coefficientDim, params).optimize();

Возвращает коэффициенты линейной модели после оптимизации.

/**
 * optimizer api.
 *
 * @return the coefficient of linear problem.
 */
@Override
public DataSet <Tuple2 <DenseVector, double[]>> optimize() {

   /**
    * solving problem using iteration.
    * trainData is the distributed samples.
    * initCoef is the initial model coefficient, which will be broadcast to every worker.
    * objFuncSet is the object function in dataSet format
    * .add(new PreallocateCoefficient(OptimName.currentCoef)) allocate memory for current coefficient
    * .add(new PreallocateCoefficient(OptimName.minCoef))     allocate memory for min loss coefficient
    * .add(new PreallocateLossCurve(OptimVariable.lossCurve)) allocate memory for loss values
    * .add(new PreallocateVector(OptimName.dir ...))          allocate memory for dir
    * .add(new PreallocateVector(OptimName.grad))             allocate memory for grad
    * .add(new PreallocateSkyk())                             allocate memory for sK yK
    * .add(new CalcGradient(objFunc))                         calculate local sub gradient
    * .add(new AllReduce(OptimName.gradAllReduce))            sum all sub gradient with allReduce
    * .add(new CalDirection())                                get summed gradient and use it to calc descend dir
    * .add(new CalcLosses(objFunc, OptimMethod.GD))           calculate local losses for line search
    * .add(new AllReduce(OptimName.lossAllReduce))            sum all losses with allReduce
    * .add(new UpdateModel(maxIter, epsilon ...))             update coefficient
    * .setCompareCriterionOfNode0(new IterTermination())             judge stop of iteration
    */
   DataSet <Row> model = new IterativeComQueue()
      .initWithPartitionedData(OptimVariable.trainData, trainData)
      .initWithBroadcastData(OptimVariable.model, coefVec)
      .initWithBroadcastData(OptimVariable.objFunc, objFuncSet)
      .add(new PreallocateCoefficient(OptimVariable.currentCoef))
      .add(new PreallocateCoefficient(OptimVariable.minCoef))
      .add(new PreallocateLossCurve(OptimVariable.lossCurve, maxIter))
      .add(new PreallocateVector(OptimVariable.dir, new double[] {0.0, OptimVariable.learningRate}))
      .add(new PreallocateVector(OptimVariable.grad))
      .add(new PreallocateSkyk(OptimVariable.numCorrections))
      .add(new CalcGradient())
      .add(new AllReduce(OptimVariable.gradAllReduce))
      .add(new CalDirection(OptimVariable.numCorrections))
      .add(new CalcLosses(OptimMethod.LBFGS, numSearchStep))
      .add(new AllReduce(OptimVariable.lossAllReduce))
      .add(new UpdateModel(params, OptimVariable.grad, OptimMethod.LBFGS, numSearchStep))
      .setCompareCriterionOfNode0(new IterTermination())
      .closeWith(new OutputModel())
      .setMaxIter(maxIter)
      .exec();

   return model.mapPartition(new ParseRowModel());
}

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

3.1 Как распределена реализация

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

Давайте посмотрим на шаги, которые можно распределить и распараллелить в L-BFGS:

Рассчитать градиентЕго можно распараллелить, например, вычислить градиент 1 на машине 1, вычислить градиент 2 на машине 2.... и затем объединить их в полный вектор градиента с помощью Reduce.
Рассчитать направлениеЕго можно распараллелить, а также разделить на разделы данных, тогда Map вычисляет значения на соответствующих разделах, а Reduces объединяются для получения направления.
Рассчитать потериЕго можно распараллелить, а потери можно получить и путем партиционирования данных, тогда Map вычисляет значения на соответствующих партициях, а Reduce комбинируется.

В Alink функция AllReduce используется вместо собственной функции Map/Reduce Flink для выполнения параллельных вычислений и передачи данных по трем вышеуказанным точкам.

3.2 CalcGradient

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

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

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

CalcGradient является производным классом оператора итерации Alink.Функции резюмируются следующим образом:

получить обработанные входные данные с меткойVectors,
Получить «параметр итерации coef», «функция оптимизации objFunc» и «градиент» из статической памяти;
Вычислить локальный градиент objFunc.calcGradient(labledVectors, coef, grad.f0) Здесь вызывается связанный с градиентом API целевой функции;objFunc.calcGradientГрадиенты рассчитываются на основе точек выборки. Alink здесь состоит в том, чтобы заменить x, y, и найти градиент функции потерь, чтобы найти частную производную по отношению к Coef. В частности, мы упоминали об этом ранее.
- Для каждой выборки в расчетных выборках расчетные градиенты различных признаков рассчитываются отдельно.
После умножения вновь рассчитанного градиента на вес он сохраняется в статической памяти gradAllReduce.
Впоследствии градиенты признаков всех рассчитанных выборок будут накапливаться с помощью функции агрегирования для получения кумулятивного градиента и потери каждого признака.

public class CalcGradient extends ComputeFunction {

    /**
     * object function class, it supply the functions to calc local gradient (or loss).
     */
    private OptimObjFunc objFunc;

    @Override
    public void calc(ComContext context) {
        Iterable<Tuple3<Double, Double, Vector>> labledVectors = context.getObj(OptimVariable.trainData);
      
// 经过处理的输入数据      
labledVectors = {ArrayList@9877}  size = 4
 0 = {Tuple3@9895} "(1.0,16.8,1.0 1.0 1.4657097546055162 1.4770978917519928)"
 1 = {Tuple3@9896} "(1.0,6.7,1.0 1.0 -0.338240712601273 -0.7385489458759964)"
 2 = {Tuple3@9897} "(1.0,6.9,1.0 1.0 -0.7892283294029703 -0.3692744729379982)"
 3 = {Tuple3@9898} "(1.0,8.0,1.0 1.0 -0.338240712601273 -0.3692744729379982)"
      
        // get iterative coefficient from static memory.
        Tuple2<DenseVector, Double> state = context.getObj(OptimVariable.currentCoef);
        int size = state.f0.size();
      
// 是Coef，1.7976931348623157E308是默认最大值
state = {Tuple2@9878} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)"
 f0 = {DenseVector@9879} "0.001 0.0 0.0 0.0"
 f1 = {Double@9889} 1.7976931348623157E308
  
        DenseVector coef = state.f0;
        if (objFunc == null) {
            objFunc = ((List<OptimObjFunc>)context.getObj(OptimVariable.objFunc)).get(0);
        }
      
// 变量如下       
objFunc = {UnaryLossObjFunc@9882} 
 unaryLossFunc = {SquareLossFunc@9891} 
 l1 = 0.0
 l2 = 0.0
 params = {Params@9892} "Params {featureCols=["f0","f1","f2"], labelCol="label", predictionCol="linpred"}"      
      
        Tuple2<DenseVector, double[]> grad = context.getObj(OptimVariable.dir);

// 变量如下      
grad = {Tuple2@9952} "(0.0 0.0 0.0 0.0,[0.0, 0.1])"
 f0 = {DenseVector@9953} "0.0 0.0 0.0 0.0"
 f1 = {double[2]@9969}       
coef = {DenseVector@9951} "0.001 0.0 0.0 0.0"
 data = {double[4]@9982} 
      
        // calculate local gradient，使用目标函数
        Double weightSum = objFunc.calcGradient(labledVectors, coef, grad.f0);

        // prepare buffer vec for allReduce. the last element of vec is the weight Sum.
        double[] buffer = context.getObj(OptimVariable.gradAllReduce);
        if (buffer == null) {
            buffer = new double[size + 1];
            context.putObj(OptimVariable.gradAllReduce, buffer);
        }

        for (int i = 0; i < size; ++i) {
            buffer[i] = grad.f0.get(i) * weightSum;
        }

      /* the last element is the weight value */
        buffer[size] = weightSum;
    }
}

// 最后结果是
buffer = {double[5]@9910} 
 0 = -38.396
 1 = -38.396
 2 = -14.206109929253465
 3 = -14.364776997288134
 4 = 0.0

3.3 AllReduce

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

.add(new AllReduce(OptimVariable.gradAllReduce))

Подробнее о том, как работает AllReduce, читайте в статье

Alink's Talk (3) Обмен коммуникационной моделью AllReduce

Alink Talks (3) AllReduce Коммуникационная модель AllReduce

3.4 CalDirection

Градиент, полученный в это время, уже агрегирован, поэтому можно рассчитать направление.

3.4.1 Предварительное распределение

OptimVariable.grad предварительно выделен.

public class PreallocateSkyk extends ComputeFunction {
    private int numCorrections;

    /**
     * prepare hessian matrix of lbfgs method. we allocate memory fo sK, yK at first iteration step.
     *
     * @param context context of iteration.
     */
    @Override
    public void calc(ComContext context) {
        if (context.getStepNo() == 1) {
            Tuple2<DenseVector, double[]> grad = context.getObj(OptimVariable.grad);
            int size = grad.f0.size();
            DenseVector[] sK = new DenseVector[numCorrections];
            DenseVector[] yK = new DenseVector[numCorrections];
            for (int i = 0; i < numCorrections; ++i) {
                sK[i] = new DenseVector(size);
                yK[i] = new DenseVector(size);
            }
            context.putObj(OptimVariable.sKyK, Tuple2.of(sK, yK));
        }
    }
}

3.4.2 Направление расчета

В процессе расчета необходимо непрерывно вычислять и хранить историческую матрицу Гессе.В алгоритме L-BFGS предполагается, что только последние m итераций информации могут быть сохранены для соответствия матрице Гессе. В алгоритме L-BFGS полный Hk больше не хранится, но сохраняются векторные последовательности {sk} и {yk} Когда требуется матрица, Hk можно вычислить, используя векторную последовательность {sk} и {yk} и последовательность векторов. Не все {sk} и {yk} нужно сохранять, просто сохраните последний вектор с m-шагами.

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

Важно отметить, что различные варианты использования данных dir:

dir = {Tuple2@9931} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[4.0, 0.1])"
 f0 = {DenseVector@9954} "-9.599 -9.599 -3.5515274823133662 -3.5911942493220335" //梯度
 f1 = {double[2]@9938} 
  0 = 4.0 //权重
  1 = 0.1 //学习率 learning rate，0.1是初始化数值，后续UpdateModel时候会更新

Резюме кода выглядит следующим образом:

@Override
public void calc(ComContext context) {
   Tuple2 <DenseVector, double[]> grad = context.getObj(OptimVariable.grad);
   Tuple2 <DenseVector, double[]> dir = context.getObj(OptimVariable.dir);
   Tuple2 <DenseVector[], DenseVector[]> hessian = context.getObj(OptimVariable.sKyK);
   int size = grad.f0.size();
   // gradarr是上一阶段CalcGradient的结果
   double[] gradarr = context.getObj(OptimVariable.gradAllReduce);
// 变量为
gradarr = {double[5]@9962} 
 0 = -38.396
 1 = -38.396
 2 = -14.206109929253465
 3 = -14.364776997288134
 4 = 4.0
   
   if (this.oldGradient == null) {
      oldGradient = new DenseVector(size);
   }
   // hessian用来当作队列，存储sK，yK，只保留最近m个
   DenseVector[] sK = hessian.f0;
   DenseVector[] yK = hessian.f1;
   for (int i = 0; i < size; ++i) {
      //gradarr[size]是权重
      grad.f0.set(i, gradarr[i] / gradarr[size]); //size = 4
   }
// 赋值梯度，这里都除以权重
grad = {Tuple2@9930} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[0.0])"
 f0 = {DenseVector@9937} "-9.599 -9.599 -3.5515274823133662 -3.5911942493220335"
  data = {double[4]@9963} 
   0 = -9.599
   1 = -9.599
   2 = -3.5515274823133662
   3 = -3.5911942493220335
 f1 = {double[1]@9961} 
  0 = 0.0  
  
   dir.f1[0] = gradarr[size]; //权重
   int k = context.getStepNo() - 1;

   if (k == 0) { //首次迭代
      dir.f0.setEqual(grad.f0); // 梯度赋予在这里
      oldGradient.setEqual(grad.f0);
   } else {
      yK[(k - 1) % m].setEqual(grad.f0);
      yK[(k - 1) % m].minusEqual(oldGradient);
      oldGradient.setEqual(grad.f0);
   }
   // copy g_k and store in qL

   dir.f0.setEqual(grad.f0); //拷贝梯度到这里
//  
dir = {Tuple2@9931} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[4.0, 0.1])"
 f0 = {DenseVector@9954} "-9.599 -9.599 -3.5515274823133662 -3.5911942493220335" //梯度
 f1 = {double[2]@9938} 
  0 = 4.0 //权重
  1 = 0.1 //学习率 learning rate，0.1是初始化数值
  
   // compute H^-1 * g_k
   int delta = k > m ? k - m : 0;
   int l = k <= m ? k : m; // m = 10
   if (alpha == null) {
      alpha = new double[m];
   }
   // two-loop的过程，通过拟牛顿法计算Hessian矩阵       
   for (int i = l - 1; i >= 0; i--) {
      int j = (i + delta) % m;
      double dot = sK[j].dot(yK[j]);
      if (Math.abs(dot) > 0.0) {
         double rhoJ = 1.0 / dot;
         alpha[i] = rhoJ * (sK[j].dot(dir.f0)); // 计算alpha
         dir.f0.plusScaleEqual(yK[j], -alpha[i]); // 重新修正d
      }
   }
   for (int i = 0; i < l; i++) {
      int j = (i + delta) % m;
      double dot = sK[j].dot(yK[j]);
      if (Math.abs(dot) > 0.0) {
         double rhoJ = 1.0 / dot;
         double betaI = rhoJ * (yK[j].dot(dir.f0)); // 乘以rho
         dir.f0.plusScaleEqual(sK[j], (alpha[i] - betaI));// 重新修正d
      }
   }
}

//最后是存储在 OptimVariable.dir

3.5 CalcLosses

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

Основная логика CalcLosses следующая:

1) Получите этот шаг Double beta = dir.f1[1] / numSearchStep;
- Скорость обучения следующего вычисления будет обновлена в последующей модели обновления, например, dir.f1[1] *= 1,0/(numSearchStep * numSearchStep) или dir.f1[1] = Math.min(dir.f1 [1 ], число шагов поиска);
2) Вызвать calcSearchValues целевой функции для расчета убытка, соответствующего текущему коэффициенту;
3) calcSearchValues обходит входные вектора меток и вычисляет потери в соответствии с шагами линейного поиска для каждого вектора меток. vec[i] += weight * this.unaryLossFunc.loss(etaCoef - i * etaDelta, labelVector.f1); Внутренняя часть цикла выглядит следующим образом:
- 3.1) Y рассчитывается по x-vec и coef, etaCoef = getEta(labelVector, coefVector);
- 3.2) deltaY вычисляется по x-vec и dirVec, etaDelta = getEta(labelVector, dirVec) * beta;
- 3.3) Рассчитать потери по шагам линейного поиска. vec[i] += weight * this.unaryLossFunc.loss(etaCoef - i * etaDelta, labelVector.f1); Из функции потерь мы знаем, что etaCoef - i * etaDelta, labelVector.f1 — это отклонение между предсказанным значением обучающие данные и актуальная категория;
4) подготовить данные для последующей агрегации lossAllReduce;

код показывает, как показано ниже:

public class CalcLosses extends ComputeFunction {

    @Override
    public void calc(ComContext context) {
        Iterable<Tuple3<Double, Double, Vector>> labledVectors = context.getObj(OptimVariable.trainData);
        Tuple2<DenseVector, double[]> dir = context.getObj(OptimVariable.dir);
        Tuple2<DenseVector, Double> coef = context.getObj(OptimVariable.currentCoef);
        if (objFunc == null) {
            objFunc = ((List<OptimObjFunc>)context.getObj(OptimVariable.objFunc)).get(0);
        }
        /**
         *  calculate losses of current coefficient.
         *  if optimizer is owlqn, constraint search will used, else without constraint.
         */
        Double beta = dir.f1[1] / numSearchStep;
        double[] vec = method.equals(OptimMethod.OWLQN) ?
            objFunc.constraintCalcSearchValues(labledVectors, coef.f0, dir.f0, beta, numSearchStep)
            : objFunc.calcSearchValues(labledVectors, coef.f0, dir.f0, beta, numSearchStep);

// 变量为
dir = {Tuple2@9988} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[4.0, 0.1])"
coef = {Tuple2@9989} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)"
beta = {Double@10014} 0.025
vec = {double[5]@10015} 
 0 = 0.0
 1 = 0.0
 2 = 0.0
 3 = 0.0
 4 = 0.0  
      
        // prepare buffer vec for allReduce.
        double[] buffer = context.getObj(OptimVariable.lossAllReduce);
        if (buffer == null) {
            buffer = vec.clone();
            context.putObj(OptimVariable.lossAllReduce, buffer);
        } else {
            System.arraycopy(vec, 0, buffer, 0, vec.length);
        }
    }
}

Поиск обеспечивается унарной целевой функцией, которая, в свою очередь, вызывает функцию потерь.

public class UnaryLossObjFunc extends OptimObjFunc {
   /**
     * Calculate loss values for line search in optimization.
     *
     * @param labelVectors train data.
     * @param coefVector   coefficient of current time.
     * @param dirVec       descend direction of optimization problem.
     * @param beta         step length of line search.
     * @param numStep      num of line search step.
     * @return double[] losses.
     */
    @Override
    public double[] calcSearchValues(Iterable<Tuple3<Double, Double, Vector>> labelVectors,
                                     DenseVector coefVector,
                                     DenseVector dirVec,
                                     double beta,
                                     int numStep) {
        double[] vec = new double[numStep + 1];
      
// labelVector是三元组Tuple3<weight, label, feature vector>      
labelVectors = {ArrayList@10007}  size = 4
 0 = {Tuple3@10027} "(1.0,16.8,1.0 1.0 1.4657097546055162 1.4770978917519928)"
 1 = {Tuple3@10034} "(1.0,6.7,1.0 1.0 -0.338240712601273 -0.7385489458759964)"
 2 = {Tuple3@10035} "(1.0,6.9,1.0 1.0 -0.7892283294029703 -0.3692744729379982)"
 3 = {Tuple3@10036} "(1.0,8.0,1.0 1.0 -0.338240712601273 -0.3692744729379982)"
coefVector = {DenseVector@10008} "0.001 0.0 0.0 0.0"
dirVec = {DenseVector@10009} "-9.599 -9.599 -3.5515274823133662 -3.5911942493220335"
beta = 0.025
numStep = 4
vec = {double[5]@10026} 
 0 = 0.0
 1 = 0.0
 2 = 0.0
 3 = 0.0
 4 = 0.0     
   
        for (Tuple3<Double, Double, Vector> labelVector : labelVectors) {
            double weight = labelVector.f0;
            //用x-vec和coef计算出来的 Y
            double etaCoef = getEta(labelVector, coefVector); 
            //以x-vec和dirVec计算出来的 deltaY
            double etaDelta = getEta(labelVector, dirVec) * beta;
weight = 1.0
etaCoef = 0.001
etaDelta = -0.7427013482280431          
            for (int i = 0; i < numStep + 1; ++i) {
                //labelVector.f1就是label y
                //联系下面损失函数可知，etaCoef - i * etaDelta, labelVector.f1 是 训练数据预测值 与 实际类别 的偏差
                vec[i] += weight * this.unaryLossFunc.loss(etaCoef - i * etaDelta, labelVector.f1); 
            }
        }
        return vec;
    }  
  
    private double getEta(Tuple3<Double, Double, Vector> labelVector, DenseVector coefVector) {
      //labelVector.f2 = {DenseVector@9972} "1.0 1.0 1.4657097546055162 1.4770978917519928"
        return MatVecOp.dot(labelVector.f2, coefVector);
    }
}

vec = {double[5]@10026} 
 0 = 219.33160199999998
 1 = 198.85962729259512
 2 = 179.40202828917856
 3 = 160.95880498975038
 4 = 143.52995739431051

Просмотрите функцию потерь следующим образом.

/**
 * Squared loss function.
 * https://en.wikipedia.org/wiki/Loss_functions_for_classification#Square_loss
 */
public class SquareLossFunc implements UnaryLossFunc {
   public SquareLossFunc() { }

   @Override
   public double loss(double eta, double y) {
      return 0.5 * (eta - y) * (eta - y);
   }

   @Override
   public double derivative(double eta, double y) {
      return eta - y;
   }

   @Override
   public double secondDerivative(double eta, double y) {
      return 1;
   }
}

3.6 UpdateModel

Этот модуль делает две вещи

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

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

OptimVariable.dir — скорректированный результат градиента, рассчитанного с помощью CalcGradient.
OptimVariable.lossAllReduce Это изменится, на этот раз это потеря, рассчитанная на предыдущем этапе.

Логика кода примерно такая:

1) Получить минимальную позицию pos "последнего убытка"
2) Получить скорость обучения бета = dir.f1[1] / numSearchStep;
3) Отдельная обработка в соответствии с оценкой «последнего значения убытка» и последнего минимального значения последнего значения убытка:
- 3.1) Если все потери больше последнего значения потерь, то
  - 3.1.1) Уменьшите скорость обучения, умножив 1,0 / (numSearchStep * numSearchStep)
  - 3.1.2) Установите eta на 0, это размер шага
  - 3.1.3) Установите текущую потерю на последнее значение потери
- 3.2) Если loss[numSearchStep] меньше последнего значения проигрыша, то
  - 3.2.1) Увеличьте скорость обучения, умножив numSearchStep
  - 3.2.2) Установите eta как наименьшее значение pos, eta = beta * pos; это eta — размер шага
  - 3.2.3) Установите текущие потери на минимальное значение текущих потерь loss[numSearchStep]
- 3.3) еще
  - 3.3.1) Скорость обучения не меняется
  - 3.3.2) Установите eta как наименьшее значение pos, eta = beta * pos; это eta представляет собой размер шага
  - 3.3.3) Установите текущие потери равными минимальным значениям текущих потерь loss[pos]
4) Исправьте матрицу Гессе
5) Используйте вектор направления и размер шага для обновления вектора коэффициентов curCoef.f0.plusScaleEqual(dir.f0, -eta);
6) Если текущий убыток меньше минимального убытка, обновите минимальный убыток текущим убытком.
- 6.1) minCoef.f1 = curCoef.f1 обновить минимальный убыток
- 6.2) minCoef.f0.setEqual(curCoef.f0); Обновите Coef, соответствующий минимальным потерям, то есть последний коэффициент, требуемый линейной моделью

Искал здесь размер шага, не нашел использования критерия Вольфа-Пауэлла, Alink сделал какую-то оптимизацию. Если есть друзья, которые могут видеть, как используется Вольф-Пауэлл, пожалуйста, оставьте сообщение для руководства, спасибо.

код показывает, как показано ниже:

public class UpdateModel extends ComputeFunction {
    @Override
    public void calc(ComContext context) {
        double[] losses = context.getObj(OptimVariable.lossAllReduce);
        Tuple2<DenseVector, double[]> dir = context.getObj(OptimVariable.dir);
        Tuple2<DenseVector, double[]> pseGrad = context.getObj(OptimVariable.pseGrad);
        Tuple2<DenseVector, Double> curCoef = context.getObj(OptimVariable.currentCoef);
        Tuple2<DenseVector, Double> minCoef = context.getObj(OptimVariable.minCoef);

        double lossChangeRatio = 1.0;
        double[] lossCurve = context.getObj(OptimVariable.lossCurve);
        for (int j = 0; j < losses.length; ++j) {
            losses[j] /= dir.f1[0]; //dir.f1[0]是权重
        }
        int pos = -1;
        //get the min value of losses, and remember the position.
        for (int j = 0; j < losses.length; ++j) {
            if (losses[j] < losses[0]) {
                losses[0] = losses[j];
                pos = j;
            }
        }

        // adaptive learningRate strategy
        double beta = dir.f1[1] / numSearchStep;
        double eta;
      
// 变量如下      
losses = {double[5]@10001} 
 0 = 35.88248934857763
 1 = 49.71490682314878
 2 = 44.85050707229464
 3 = 40.239701247437594
 4 = 35.88248934857763
dir = {Tuple2@10002} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[4.0, 0.1])"
pseGrad = null
curCoef = {Tuple2@10003} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)"
minCoef = {Tuple2@10004} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)"
lossChangeRatio = 1.0
lossCurve = {double[100]@10005} 
pos = 4
beta = 0.025
curCoef.f1 = {Double@10006} 1.7976931348623157E308 
      
        if (pos == -1) {
            /**
             * if all losses larger than last loss value. we'll do the below things:
             * 1. reduce learning rate by multiply 1.0 / (numSearchStep*numSearchStep).
             * 2. set eta with zero.
             * 3. set current loss equals last loss value.
             */
            eta = 0;
            dir.f1[1] *= 1.0 / (numSearchStep * numSearchStep);  // 学习率 
            curCoef.f1 = losses[0]; // 最小loss
        } else if (pos == numSearchStep) {
            /**
             * if losses[numSearchStep] smaller than last loss value. we'll do the below things:
             * 1. enlarge learning rate by multiply numSearchStep.
             * 2. set eta with the smallest value pos.
             * 3. set current loss equals smallest loss value.
             */
            eta = beta * pos;
            dir.f1[1] *= numSearchStep;
            dir.f1[1] = Math.min(dir.f1[1], numSearchStep); // 学习率
            lossChangeRatio = Math.abs((curCoef.f1 - losses[pos]) / curCoef.f1);
            curCoef.f1 = losses[numSearchStep]; // 最小loss
          
// 当前数值
numSearchStep = 4 是NUM_SEARCH_STEP缺省值，是线性搜索的参数，Line search parameter, which define the search value num of one step.        
lossChangeRatio = 1.0
pos = 4
eta = 0.1
curCoef.f1 = {Double@10049} 35.88248934857763   
dir.f1[1] = 0.4  
  
        } else {
            /**
             * else :
             * 1. learning rate not changed.
             * 2. set eta with the smallest value pos.
             * 3. set current loss equals smallest loss value.
             */
            eta = beta * pos;
            lossChangeRatio = Math.abs((curCoef.f1 - losses[pos]) / curCoef.f1);
            curCoef.f1 = losses[pos]; // 最小loss
        }

		/* update loss value in loss curve at this step */
        lossCurve[context.getStepNo() - 1] = curCoef.f1;
      
lossCurve = {double[100]@9998} 
 0 = 35.88248934857763
 1 = Infinity
   
        if (method.equals(OptimMethod.OWLQN)) {
           .....          
        } else if (method.equals(OptimMethod.LBFGS)) {          
            Tuple2<DenseVector[], DenseVector[]> sKyK = context.getObj(OptimVariable.sKyK);
            int size = dir.f0.size();
            int k = context.getStepNo() - 1;
            DenseVector[] sK = sKyK.f0;
            for (int s = 0; s < size; ++s) {
                // 修正矩阵
                sK[k % OptimVariable.numCorrections].set(s, dir.f0.get(s) * (-eta));
            }
            curCoef.f0.plusScaleEqual(dir.f0, -eta); // 这里是用方向向量和步长来更新系数向量
          
sKyK = {Tuple2@10043} "([0.9599000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0],[0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0])"
 f0 = {DenseVector[10]@10044} 
  0 = {DenseVector@10074} "0.9599000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034"
    
        } 
      
        /**
         * if current loss is smaller than min loss, then update the min loss and min coefficient by current.
         */
        if (curCoef.f1 < minCoef.f1) {
            minCoef.f1 = curCoef.f1; // 最小loss
            minCoef.f0.setEqual(curCoef.f0); // 最小loss所对应的Coef，即线性模型最后需要的系数
        }

curCoef = {Tuple2@9996} "(0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034,35.88248934857763)"
 f0 = {DenseVector@10059} "0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034"
 f1 = {Double@10048} 35.88248934857763
minCoef = {Tuple2@9997} "(0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034,35.88248934857763)"
 f0 = {DenseVector@10059} "0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034"
 f1 = {Double@10048} 35.88248934857763
      
        // judge the convergence of iteration.
        filter(dir, curCoef, minCoef, context, lossChangeRatio);
    }
}

Фильтр оценивает, сходится он или нет, а журнал печати внутри ясно иллюстрирует логику функции.

/**
 * judge the convergence of iteration.
 */
public void filter(Tuple2<DenseVector, double[]> grad,
                   Tuple2<DenseVector, Double> c,
                   Tuple2<DenseVector, Double> m,
                   ComContext context,
                   double lossChangeRatio) {
    double epsilon = params.get(HasEpsilonDv0000001.EPSILON);
    int maxIter = params.get(HasMaxIterDefaultAs100.MAX_ITER);
    double gradNorm = ((Tuple2<DenseVector, double[]>)context.getObj(gradName)).f0.normL2();
    if (c.f1 < epsilon || gradNorm < epsilon) {
        printLog(" method converged at step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio);
        grad.f1[0] = -1.0;
    } else if (context.getStepNo() > maxIter - 1) {
        printLog(" method stop at max step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio);
        grad.f1[0] = -1.0;
    } else if (grad.f1[1] < EPS) {
        printLog(" learning rate is too small, method stops at step : ", c.f1, m.f1, grad.f1[1], gradNorm,
            context, lossChangeRatio);
        grad.f1[0] = -1.0;
    } else if (lossChangeRatio < epsilon && gradNorm < Math.sqrt(epsilon)) {
        printLog(" loss change ratio is too small, method stops at step : ", c.f1, m.f1, grad.f1[1], gradNorm,
            context, lossChangeRatio);
        grad.f1[0] = -1.0;
    } else {
        printLog(" method continue at step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio);
    }
}

3.7 OutputModel

.closeWith(new OutputModel())Эта часть предназначена для временного вывода модели в конце каждой итерации, преобразования данных в общий тип строки Flink и передачи состояния в строки модели.

public class OutputModel extends CompleteResultFunction {

   @Override
   public List <Row> calc(ComContext context) {
      // get the coefficient of min loss.
      Tuple2 <DenseVector, double[]> minCoef = context.getObj(OptimVariable.minCoef);
      double[] lossCurve = context.getObj(OptimVariable.lossCurve);

      int effectiveSize = lossCurve.length;
      for (int i = 0; i < lossCurve.length; ++i) {
         if (Double.isInfinite(lossCurve[i])) {
            effectiveSize = i;
            break;
         }
      }

      double[] effectiveCurve = new double[effectiveSize];
      System.arraycopy(lossCurve, 0, effectiveCurve, 0, effectiveSize);

      Params params = new Params();
      params.set(ModelParamName.COEF, minCoef.f0);// 重点在这里，minCoef是我们真正需要的
      params.set(ModelParamName.LOSS_CURVE, effectiveCurve);
      List <Row> model = new ArrayList <>(1);
      model.add(Row.of(params.toJson()));
      return model;
   }
}

0x04 Подготовить метаданные модели

Здесь параллелизм установлен на 1.

// Prepare the meta info of linear model.
DataSet<Params> meta = labelInfo.f0
    .mapPartition(new CreateMeta(modelName, linearModelType, isRegProc, params))
    .setParallelism(1);

специальный код

public static class CreateMeta implements MapPartitionFunction<Object, Params> {
        @Override
        public void mapPartition(Iterable<Object> rows, Collector<Params> metas) throws Exception {
            Object[] labels = null;
            if (!this.isRegProc) {
                labels = orderLabels(rows);
            }

            Params meta = new Params();
            meta.set(ModelParamName.MODEL_NAME, this.modelName);
            meta.set(ModelParamName.LINEAR_MODEL_TYPE, this.modelType);
            meta.set(ModelParamName.LABEL_VALUES, labels);
            meta.set(ModelParamName.HAS_INTERCEPT_ITEM, this.hasInterceptItem);
            meta.set(ModelParamName.VECTOR_COL_NAME, vectorColName);
            meta.set(LinearTrainParams.LABEL_COL, labelName);
            metas.collect(meta);
        }  
}

// 变量为
meta = {Params@9667} "Params {hasInterceptItem=true, vectorColName=null, modelName="Linear Regression", labelValues=null, labelCol="label", linearModelType="LinearReg"}"

0x05 Построить модель

Когда итерационный цикл завершен, Alink строит модель на основе данных Coef.

/**
 * build the linear model rows, the format to be output.
 */
public static class BuildModelFromCoefs extends AbstractRichFunction implements
        @Override
        public void mapPartition(Iterable<Tuple2<DenseVector, double[]>> iterable,
                                 Collector<Row> collector) throws Exception {
            for (Tuple2<DenseVector, double[]> coefVector : iterable) {
                LinearModelData modelData = buildLinearModelData(meta,
                    featureNames,
                    labelType,
                    meanVar,
                    hasIntercept,
                    standardization,
                    coefVector);

                new LinearModelDataConverter(this.labelType).save(modelData, collector);
            }
        }  
}

Получены модельные данные, где coef равен w, b в f(x)=w^Tx+b. наконец, используется для расчета.

modelData = {LinearModelData@10584} 
 featureNames = {String[3]@9787} 
 featureTypes = null
 vectorColName = null
 coefVector = {DenseVector@10485} "-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847"
 coefVectors = null
 vectorSize = 3
 modelName = "Linear Regression"
 labelName = null
 labelValues = null
 linearModelType = {LinearModelType@4674} "LinearReg"
 hasInterceptItem = true
 lossCurve = {double[12]@10593} 
  0 = 35.88248934857763
  1 = 12.807366842002144
  2 = 0.5228366663917704
  3 = 0.031112070740366038
  4 = 0.01098914933042993
  5 = 0.009765757443537283
  6 = 0.008750523231785415
  7 = 0.004210085397869248
  8 = 0.0039042232755530704
  9 = 0.0038821509860327537
  10 = 0.003882042680010676
  11 = 0.0038820422536391033
 labelType = {FractionalTypeInfo@9790} "Double"

0x06 Использовать предсказание модели

При прогнозировании используется LinearModelMapper, и некоторые его внутренние переменные выводятся следующим образом, и можно увидеть данные модели.

this = {LinearModelMapper@10704} 
 vectorColIndex = -1
model = {LinearModelData@10597} 
 featureNames = {String[3]@10632} 
 featureTypes = null
 vectorColName = null
 coefVector = {DenseVector@10612} "-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847"
 coefVectors = null
 vectorSize = 0
 modelName = "Linear Regression"
 labelName = null
 labelValues = {Object[0]@10613} 
 linearModelType = {LinearModelType@10611} "LinearReg"
 hasInterceptItem = true
 lossCurve = null
 labelType = null

Конкретный прогнозLinearModelMapper.predictЗавершение завершено, следующим образом:

соответствует исходному вводуRow.of("$3$0:1.0 1:7.0 2:9.0", "1.0 7.0 9.0", 1.0, 7.0, 9.0, 16.8),,
пройти черезFeatureLabelUtil.getFeatureVectorПосле обработки,
Полученная четверка"1.0 1.0 7.0 9.0", где передается первый 1.0aVector.set(0, 1.0);Специально установленное фиксированное значение. Например, модель f(x) = ax + b, фиксированное значение 1,0 — это значение инициализации b, а b будет получено в процессе оптимизации. Так что для прогнозов все еще нужна версия 1.0.
Коэффициенты модели:"-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847".
Результат скалярного произведения кватерниона и коэффициентов модели равенdotValue = 16.814789059973744.

Это показывает, как используются коэффициенты модели.

public Object predict(Vector vector) throws Exception {
   double dotValue = MatVecOp.dot(vector, model.coefVector);

   switch (model.linearModelType) {
      case LR:
      case SVM:
      case Perceptron:
         return dotValue >= 0 ? model.labelValues[0] : model.labelValues[1];
      case LinearReg:
      case SVR:
         return dotValue;
   }
}

vector = {DenseVector@10610} "1.0 1.0 7.0 9.0"
 data = {double[4]@10619} 
  0 = 1.0
  1 = 1.0
  2 = 7.0
  3 = 9.0
model.coefVector = {DenseVector@10612} "-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847"
 data = {double[4]@10620} 
  0 = -3.938937407856857
  1 = 4.799499941426075
  2 = 0.8929571907809862
  3 = 1.078169576770847

ссылка 0xFF

Вуху.. Mommycode.com/info-detail…)

Анализ реализации алгоритма машинного обучения - алгоритм L-BFGS liblbfgs

BFGS и L-BFGS из серии Deep Machine Learning

Вывод формулы квазиньютоновского метода и реализация кода на Python (1)

Легко понять - разложение Тейлора

Расширение Тейлора — Расширение Тейлора

От интерполяции Ньютона к формуле Тейлора

Как лучше понять и запомнить разложения Тейлора?

Алгоритм оптимизации - метод Ньютона

Алгоритм оптимизации - алгоритм L-BFGS квазиньютоновского метода

Одна статья для понимания алгоритма L-BFGS

«Алгоритмы распределенного машинного обучения, теория и практика — Лю Теян»

LogisticRegressionWithLBFGS — логистическая регрессия

Метод оптимизации логистической регрессии — L-BFGS

Машинное обучение и оптимизация исследования операций (3) От метода Ньютона к L-BFGS

Серия углубленного машинного обучения 17-BFGS и L-BFGS

Серия "Основы методов оптимизации" - Методы точного одномерного поиска

Серия "Основы методов оптимизации" - неточные методы одномерного поиска

[Оригинал] Объясните критерии Армихо-Гольдштейна и критерия Вулфа-Пауэлла в неточном линейном поиске с помощью «человеческих слов»

Ууху. Call.com/question/36…

Введение и реализация алгоритмов одномерного поиска

Численная оптимизация | Примечания отделочные (1) - введение, поиск строки: условия выбора шага

zhuanlan.zhihu.com/p/32821110

Поиск линии (1): выбор размера шага

Задача оптимизации - одномерный поиск (2)

Принцип поиска строки CRF L-BFGS и анализ кода

Размер шага и скорость обучения

Алгоритм оптимизации машинного обучения L-BFGS и его распределенная реализация

Подробное объяснение алгоритма L-BFGS (алгоритм оптимизации по умолчанию для логистической регрессии)

Принцип и практика линейной регрессии (метод Ньютона)

Линейная регрессия, градиентный спуск

Алгоритм L-BFGS

Параллельная логистическая регрессия

★★★★★★Думая о жизни и технологиях★★★★★★

Публичный аккаунт WeChat: мысли Росси

Если вы хотите получать своевременные новости о статьях, написанных отдельными лицами, или хотите видеть технические материалы, рекомендованные отдельными лицами, обратите внимание.