ML: Линейная модель

---

Введение

В этом документе рассматривается многомерная линейная модель. Она имеет самостоятельное значение в некоторых задачах машинного обучения и лежит в основе сложных архитектур глубоких нейронных сетей. Обсуждаются математические формулировки модели и аналитический поиск её параметров. Приведены примеры кода на Python которые можно найти в файле: ML_Line_Model.ipynb.

Рассмотрим $n_x$-мерное пространство признаков. Каждый объект в этом пространстве описывается вектором $\mathbf{x}$. Задача состоит в преобразовании исходного пространства признаков в другое, $n_y$-мерное пространство $\mathbf{y}$. Величины $\mathbf{y}$, сами по себе, могут иметь практический смысл и тогда это регрессия или быть новыми признаками объектов, определяющими вероятности в задаче классификации.

Линейная модель является простейшей связью входа $\mathbf{x}$ и выхода $\mathbf{y}$. Её матричная и индексная записи имеют вид:

$$\mathbf{y} = \mathbf{b} + \mathbf{x}\cdot\mathbf{w}, ~~~~~~~~~~~~~~~~~~~~ y_\alpha = b_\alpha + \sum_\beta (x_\beta\, w_{\beta\alpha})$$

В процессе обучения необходимо подобрать оптимальные параметры модели: матрицу вектор $b_\beta$ и $w_{\alpha\beta}$. Линейная модель является первым нетривиальным приближением разложения в ряд Тейлора произвольной функции $\mathbf{y}=f(\mathbf{x}\, \boldsymbol{\omega})$.

---

Линейная модель как гиперплоскость

Рассмотрим геометрическую интерпретацию линейной модели. В $n$-мерном пространстве каждая точка описывается $n$ вещественными координатами $\mathbf{x}=\{x_1,...,x_n\}$. Как и в обычном 3-мерном пространстве, плоскость задаётся вектором нормали $\boldsymbol{\omega}=\{w_1,...,w_n\}$ (перпендикуляр к плоскости) и произвольной точкой $\mathbf{x}_0=\{x_{01},...,x_{0n}\}$, лежащей в этой плоскости (фиксирующей смещение плоскости вдоль нормали).
Когда $n > 3$ плоскость принято называть гиперплоскостью.

Расстояние $d$ от гиперплоскости до некоторой точки $\mathbf{x}=\{x_1,...,x_n\}$ вычисляется по формуле ($\boldsymbol{\omega}^2=1$):

$$d = b+\mathbf{x}\cdot\boldsymbol{\omega}, ~~~~~~\text{где}~~~~~~~~~~ b ~=~ -\mathbf{x}_0\cdot\boldsymbol{\omega} ~=~ -( x_{01}w_1 + ... + x_{0n}w_{n} ).$$

При этом $d > 0$, если точка $\mathbf{x}$ лежит с той стороны гиперплоскости, куда указывает вектор $\boldsymbol{\omega}$ и $d < 0$, если с противоположной. Когда $d=0$ - точка $\mathbf{x}$ лежит в гиперплоскости ($b+\mathbf{x}\cdot\boldsymbol{\omega} = 0$ - её уравнение).

Изменение параметра $b$ сдвигает плоскость параллельным образом в пространстве. Если $b$ уменьшается, то плоскость смещается в направлении вектора $\boldsymbol{\omega}$ (расстояние $d$ меньше), а если $b$ увеличивается - плоскость смещается против вектора $\boldsymbol{\omega}$. Это непосредственно следует из приведенной выше формулы.

Если пространство имеет $n$ измерений, то гиперплоскость это $(n-1)$-мерный объект.
Она делит всё пространство на две части. Для наглядности рассмотрим 2-мерное пространство. Гиперплоскостью в нём будет прямая линия (одномерный объект). Справа на рисунке кружок изображает одну точку пространства, а крестик - другую. Они расположены по разные стороны от линии (гиперплоскости). Если длина вектора $\boldsymbol{\omega}$ много больше единицы, то и расстояния $d$ от точек к плоскости по модулю будут существенно большими единицы.

---

Регрессии (предсказание вещественного числа)

Простейшим примером регрессии является одномерная линейная функция $y = b + w \,x$, в которой параметр $w$ характеризует наклон прямой, а параметр $b$ - сдвиг (смещениие) по оси $y$ от начала координат. Наличие подобной связи между $x$ и $y$, обычно, легко обнаружить визуально на графике даннных.

В многомерном случае $y=b + w_1x_1+...w_n x_n$ визуализировать данные сложнее. Чтобы убедится в адекватности линейной модели, необходимо найти оптимальные параметры и вычислить величину ошибки (сравнив её с ошибками других моделей, в т.ч. с тривиальной $y=b=\text{const}$).

---

Логистическая регрессия (задача классификации)

В задаче классификации объект $\mathbf{x}$ необходимо отнести к одному или нескольким классам, общее число $C$ которых фиксированно. Обычно, результатом классификации является $C$ условных вероятностей $p(c|\mathbf{x})$ того, что объект с признаками $\mathbf{x}$ принадлежит $c$-тому классу: $c\in [0...C-1]$.

При использовании в задаче классификации линейной модели, её выходы $y_\alpha$ "пропускают" через нелинейную функцию $z_\alpha=f(y_\alpha)$, значения которой находятся в диапазоне $[0...1]$. Чем ближе $z_\alpha$ к единице, тем увереннее модель относит объект к классу с номером $\alpha$.

$$z_{\alpha} ~=~f(\mathbf{y}) ~=~f\Bigr(b_\alpha + \sum_{\mu} x_{\mu}\, w_{\mu \alpha} \Bigr)$$

Если объект может принадлежать к нескольким классам (классы пересекаются), то в качестве $f(z)$ удобна сигмоидная функция. Если на нескольких выходах модели получаются значения близкие к единице, то объект принадлежит всем этим классам. При больших $y_\alpha$ сигмоид стремиться к единице, а при $y_\alpha\to -\infty$ - к нулю: $\sigma(y_\alpha) = 1/(1+e^{-y_\alpha})$

Для непересекающихся классов $C > 2$ (объект $\mathbf{x}$ принадлежит только одному классу), обычно, применяют функцию softmax. На её выходе получаются положительные значения, сумма которых равна единице. Поэтому их можно интерпретировать как распределение вероятностей $p(c|\mathbf{x})$: $$z_\alpha = \mathrm{sm}(y_\alpha) = \displaystyle \frac{e^{y_\alpha}}{\sum_\beta e^{y_\beta}},~~~~~~~~~~~~~~~~~~~~~\sum_\alpha z_\alpha = 1,~~~~~~~~~z_\alpha > 0.$$

Если разделяющие плоскости для $\{y_0,y_1,y_2\}$ проведены верно и классы линейно разделимы, то для точки, относящейся к классу $0$, значения выходов будут стремиться к $\{1,0,0\}$, для класса $1$ - $\{0,1,0\}$ и т.д.

---

Формула Байеса и софтмакс

Существует замечательный математический факт, который подводит "теоретическую основу" под использование функции софтмакса. Пусть примеры каждого класса $c=[0...,C-1]$ в пространстве признаков $\mathbf{x}$ образуют "кластеры" с нормальными распределениями, которые имеют одинаковую дисперсию $\mathbf{D}$ и различные средние $\bar{\mathbf{x}}_c$. Тогда условная вероятность $p(c|\mathbf{x})$ принадлежности объекта с признаками $\mathbf{x}$ к классу $c$, описывается линейной моделью с софтмакс-функцией на выходе: $$P(\mathbf{x} | c) \sim e^{-\frac{1}{2}\,(\mathbf{x}-\bar{\mathbf{x}}_c)\,\mathbf{D}^{-1}(\mathbf{x}-\bar{\mathbf{x}}_c)}~~~~~~~~~~~~~ \Rightarrow~~~~~~~~~~~~~ p(c|\mathbf{x}) = \frac{e^{y_c}}{\sum_\alpha e^{y_\alpha}},~~~~~~~\mathbf{y} = \mathbf{b} + \mathbf{x}\mathbf{w}.$$

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

---

🔥 Из определения условной вероятности следует формула Байеса: $$p(c|\mathbf{x}) = \frac{p(c,\mathbf{x})}{p(\mathbf{x})} = \frac{p(c)\,P(\mathbf{x}|c)}{p(\mathbf{x})} = \frac{p(c)\,P(\mathbf{x}|c)}{\sum_\alpha p(\alpha)\,P(\mathbf{x}|\alpha)},$$ где $p(c)$ - вероятности принадлежности случайно выбранного объекта к классу $c$. Распишем произведение $(\mathbf{x}-\bar{\mathbf{x}}_c)\,\mathbf{D}^{-1}(\mathbf{x}-\bar{\mathbf{x}}_c)$ под экспонентой в $P(\mathbf{x} | c)$. Слагаемое $-(1/2)\,\mathbf{x}\,\mathbf{D}^{-1}\mathbf{x}$, присутствующее и в числителе, и в знаменателе формулы Байеса сокращается. Остальные слагаемые и множитель $p(c)$ можно факторизовать следующим образом (обратная матрица дисперсии $\mathbf{D}^{-1}$ симметрична): $$p(c)\,P(\mathbf{x}|c) \sim e^{b_c + \mathbf{x}\,\mathbf{w}_c},~~~~~~~~~ b_c = -\frac{1}{2}\,\bar{\mathbf{x}}_c\,\mathbf{D}^{-1}\bar{\mathbf{x}}_c + \log\, p(c),~~~~~~~~~~ \mathbf{w}_c = \mathbf{D}^{-1}\bar{\mathbf{x}}_c.$$

---

Батчи и перемножение матриц

Пусть есть 2 признака: $\{x_0,x_1\}$ (вход модели), которые определяют три целевых величины: $\{y_0,y_1,y_2\}$ (выход). Линейную связь между выходом и входом: $$\left\{\begin{array}{lcl} y_0 &=& b_0+ x_0 \,w_{00} + x_1\, w_{10} \\ y_1 &=& b_1+ x_0 \,w_{01} + x_1\, w_{11} \\ y_2 &=& b_2+ x_0 \,w_{02} + x_1\, w_{12} \end{array} \right. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \mathbf{w} = \overbrace{ \begin{array}{|c|c|c|} \hline w_{00} & w_{01} & w_{02} \\ \hline w_{10} & w_{11} & w_{12} \\ \hline \end{array} }^{\displaystyle n_y} \left. \phantom{ \begin{array}{c} \\ \\ \end{array} } \right\} n_x$$ можно записать в табличном виде (точка - это dot-умножение строки на столбец):

$$\begin{array}{|c|c|} \hline b_{0} & b_{1} & b_{2} \\ \hline \end{array} ~~~+~~~ \begin{array}{|c|c|} \hline x_0 & x_1 \\ \hline \end{array} ~\cdot~ \begin{array}{|c|c|c|} \hline w_{00} & w_{01} & w_{02} \\ \hline w_{10} & w_{11} & w_{12} \\ \hline \end{array} ~~~=~~~ \begin{array}{|c|c|} \hline y_{0} & y_{1} & y_{2} \\ \hline \end{array}$$

В алгоритмах машинного обучения вычисление происходит не с одним примером, а сразу с некоторым их множеством - пачкой (batch). Поэтому представим $\mathbf{x}$ в виде 2D матрицы, расположив примеры построчно (первый индекс, начиная с нуля - номер примера, а второй - номер признака в данном примере).
Ниже приведено $N=4$ примера, $n_x=2$ признака и $n_y=3$ выходов модели:

$$\begin{array}{|c|c|c|} \hline b_{0} & b_{1} & b_{2}\\ \hline \end{array} ~~~+~~~ \begin{array}{|c|c|} \hline x_{00} & x_{01} \\ \hline x_{10} & x_{11} \\ \hline x_{20} & x_{21} \\ \hline x_{30} & x_{31} \\ \hline \end{array} ~\cdot~ \begin{array}{|c|c|c|} \hline w_{00} & w_{01} & w_{02} \\ \hline w_{10} & w_{11} & w_{12} \\ \hline \end{array} ~~~=~~~ \left. \begin{array}{|c|c|c|} \hline y_{00} & y_{01} & y_{02} \\ \hline y_{10} & y_{11} & y_{12} \\ \hline y_{20} & y_{21} & y_{22} \\ \hline y_{30} & y_{31} & y_{32} \\ \hline \end{array} \right\} \text{примеры}$$

Размерности индексов (атрибут shape в numpy) матрицы произведения $\mathbf{x}\cdot\mathbf{w}$ равны: $$(N,n_x) . (n_x,n_y) = (N,n_y).$$ При свёртке матриц: $(\mathbf{x}\cdot\mathbf{w})_{i\alpha} = \sum_\mu x_{i\mu}\,w_{\mu\alpha}$ действует принцип "строка на столбец" и размерность среднего индекса сокращается. К каждой строке матрицы $\mathbf{x}\cdot\mathbf{w}$ формы (N,ny) добавляется вектор $\mathbf{b}$ формы (ny,).
В библиотеке numpy это соответствует правилу расширения (broadcasting).

---

Тензорный анализ

Далее будут вычисляться производные по тензорам и использоваться матричные оозначения. Поэтому напомним некоторые важные моменты.

---

☝ Матрица, обратная к матрице $\mathbf{A}$, обозначается как $\mathbf{A}^{-1}$, а транспонированная как $\mathbf{A}^{\top}$. У транспонированной матрицы колонки и строки (индексы) переставлены местами: $A^\top_{ij}=A_{ji}$, а произведение обратной на исходную даёт единичную матрицу (символ Кронекера). Они обладают следующими свойствами: $$\mathbf{A}^{-1}\cdot\mathbf{A}=\mathbf{A}\cdot\mathbf{A}^{-1}=\mathbf{1},~~~~~~~~~~ (\mathbf{A}\cdot\mathbf{B})^{-1} = \mathbf{B}^{-1}\cdot\mathbf{A}^{-1},~~~~~~~~~~ (\mathbf{A}\cdot\mathbf{B})^{\top} = \mathbf{B}^{\top}\cdot\mathbf{A}^{\top}.$$

---

☝ Справедлива формула для частной производной (многоточие - это другие индексы):

$$\frac{\partial W_{\mu\nu...}}{\partial W_{\alpha\beta...}} = \delta_{\mu\alpha}\,\delta_{\nu\beta}\,...,~~~~~~~~~~~~~~ \sum_{\mu} T_{...\mu...}\,\delta_{\mu\alpha} = T_{...\alpha...},$$

где $\delta_{\alpha\beta}$ - единичная квадратная матрица (символ Кронекера), равная $0$, если $\alpha\neq\beta$ и $1$, если $\alpha=\beta$. При суммировании с символом Кронекера (по одному из его индексов), сумма и символ Кронекера убираются, а суммационный индекс заменяется на второй индекс символа Кронекера (выше вторая формула).

Формула для $\partial W_{\mu\nu}/\partial W_{\alpha\beta}$, означает, что производная равна нулю, если $\mu\neq\alpha$ и $\nu\neq\beta$. Например, $\partial W_{12}/\partial W_{13} = 0$ так это разные переменные (типа $\partial y/\partial x = 0$), а $\partial W_{12}/\partial W_{12} = 1$ (типа $\partial x/\partial x = 1$).

---

Аналитическое решение

В линейной модели с mse-ошибкой аналитически несложно найти оптимальные параметры. Вычисления оказываются заметно проще, если вектор $\mathbf{b}$ и матрицу $\mathbf{w}$ объединить в одну матрицу $\mathbf{W}$, добавив в $\mathbf{w}$, в качестве ещё одной строки, вектор $\mathbf{b}$. Дополнительно в матрицу $\mathbf{x}$ добавляется колонка, состоящая из единиц, что даёт матрицу $\mathbf{X}$:

$$\mathbf{X}\cdot\mathbf{W} ~=~ \begin{array}{|c|c|c|} \hline x_{00} & x_{01} & ~1~ \\ \hline x_{10} & x_{11} & 1 \\ \hline x_{20} & x_{21} & 1 \\ \hline x_{30} & x_{31} & 1 \\ \hline \end{array} ~\cdot~ \begin{array}{|c|c|} \hline w_{00} & w_{01} & w_{02} \\ \hline w_{10} & w_{11} & w_{12} \\ \hline b_{0} & b_{1} & b_{2} \\ \hline \end{array} ~=~ \begin{array}{|c|c|c|} \hline y_{00} & y_{01} & y_{02} \\ \hline y_{10} & y_{11} & y_{12} \\ \hline y_{20} & y_{21} & y_{22} \\ \hline y_{30} & y_{31} & y_{32} \\ \hline \end{array} ~=~ \mathbf{y}.$$ Такая запись означает введение ещё одного фиктивного признака всегда равного единице: $\mathbf{x}=\{x_0,...,x_{n-1},1\}$.

Теперь модель имеет вид: $\mathbf{y} = \mathbf{X}\cdot\mathbf{W}$, что, как несложно проверить, эквивалентно исходной записи $\mathbf{y}=\mathbf{b} + \mathbf{x}\mathbf{w}$. Найдём минимум mse-ошибки $L$ (опускаем постоянный множитель):

$$L = \sum_{i,\nu} (y_{i\nu} -\hat{y}_{i\nu} )^2,~~~~~~~~~~~~~~~ y_{i\nu} = \sum_{\mu} X_{i\mu} W_{\mu\nu}.$$

Запишем производные ошибки и выхода модели: $$\frac{\partial L}{\partial y_{i\nu}}~=~ 2(y_{i\nu} -\hat{y}_{i\nu} ),~~~~~~~~ \frac{\partial y_{i\nu}}{\partial W_{\alpha\beta}} = \sum_{\mu} X_{i\mu}\,\delta_{\mu\alpha}\,\delta_{\nu\beta}. ~=~X_{i\alpha} \,\delta_{\nu\beta}.$$ Теперь несложно найти производную ошибки $L$ по элементам матрицы $W_{\alpha\beta}$: $$\frac{\partial L}{\partial W_{\alpha\beta}} = \sum_{i,\nu} \frac{\partial L}{\partial y_{i\nu}}\,\frac{\partial y_{i\nu}}{\partial W_{\alpha\beta}} ~=~ 2 \sum_{i} (y_{i\beta} -\hat{y}_{i\beta} )\, X_{i\alpha} ~=~2 \bigr[\mathbf{X}^\top\cdot(\mathbf{y}-\hat{\mathbf{y}} )\bigr]_{\alpha\beta}.$$

Приравняв производную нулю (условие экстремума и в данном случае минимума $L$), получаем матричное уравнение $\mathbf{X}^\top\cdot(\mathbf{X}\cdot\mathbf{W}-\hat{\mathbf{y}} ) = 0$, которое легко решается:

$$\mathbf{W} = (\mathbf{X}^\top\cdot \mathbf{X})^{-1}\cdot (\mathbf{X}^\top\cdot \hat{\mathbf{y}}) ~~~~~~~~\text{или}~~~~~~~~~~~~~~ \mathbf{W}= \mathbf{X}^{-1} \cdot \hat{\mathbf{y}}.$$

В принципе, можно пользоваться любой из этих формул. Первая формула оказывается быстрее, если примеров N существенно больше чем признаков nX = число входов $\mathbf{x}$. Искать обратную матрицу для квадратной матрицы $\mathbf{X}^\top\cdot \mathbf{X}$ формы (nX+1, nX+1) проще, чем для матрицы $\mathbf{X}^{-1}$ формы (nX+1, N) во второй формуле, хотя при этом добавляется перемножение матриц $\mathbf{X}^\top\cdot \mathbf{X}$.

---

Пример на numpy

Проведём вычисление по полученным формулам при помощи библиотеки numpy. Краткое введение в numpy можно найти в документе NN_Base_Numpy.

Cгенерим сначала модельные данные в которых выходы линейно связаны с входами и добавлен небольшой случайный шум:

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

---

Регрессия на sklearn

Для решения этой-же задачи, на практике удобнее пользоваться готовой библиотекой машинного обучения Scikit-learn (кратко sklearn):

Вызов LinearRegression() создаёт экземпляр класса данного метода обучения, которому через функцию fit передаются обучающие данные X,Y. Эта функция снова возвращает экземпляр класса, в атрибутах которого находятся параметры модели. Обратим внимание, что матрица $\mathbf{w}$ (атрибут lr.coef_) в sklearn хранится в транспонированном виде (по сравнению с принятым выше).

Получив параметры модели, можно теперь вычислить её mse-ошибку:

---

Классификация на sklearn

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

Рассмотрим для наглядности двумерное пространство признаков $\mathbf{x}=\{x_0, x_1\}$ и три класса $c=\{0,1,2\}$, к одному из которых необходимо отнести объект. У модели будет три выхода: $n_y=3$, на которых, в качестве целевых значений будем требовать для первого, второго и третьего классов соответственно: $$\mathbf{y}=\{+1,-1,-1\},~~~~~\mathbf{y}=\{-1,+1,-1\},~~~~~\mathbf{y}=\{-1,-1,+1\}.$$ По-прежнему будем использовать среднеквадратичную ошибку, для которой есть точное решение. Естественно, линейная модель не сможет для всех примеров добиться целевых значений. Однако, три выхода модели $\mathbf{y}=\{d_0,d_1,d_2\}$ - это три знакопеременных расстояния от точки $\mathbf{x}$ до трёх плоскостей. Расстояния положительны, если вектор нормали "смотрит" на точку и отрицательны - если в противоположную сторону. Метод наименьших квадратов будет стремится развернуть эти векторы, чтобы, по крайней мере, совпадали знаки выхода модели и целевых значений. Это даст необходимый классификатор.

Подготовим модельные данные в виде трёх кластеров с нормальным разбросом точек вокруг их центров:

Модель, которая проводит классификацию оформим в виде класса:

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

Чтобы получить график как в начале раздела, необходимо добавить к классу Model следующий метод:

С помощью этого метода рисуем разделяющие поверхности и обучающие объекты:

---