Processing math: 5%

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


Введение

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

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

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

y=b+xw,                    yα=bα+β(xβwβα)

В процессе обучения необходимо подобрать оптимальные параметры модели: матрицу вектор bβ и wαβ. Линейная модель является первым нетривиальным приближением разложения в ряд Тейлора произвольной функции \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}. Это непосредственно следует из приведенной выше формулы.

◄ Запишем вектор \mathbf{x}-\mathbf{x}_0, начинающийся в точке \mathbf{x}_0 (лежащей в плоскости) и направленный в точку \mathbf{x} (см. рисунок справа; векторы складываются по правилу треугольника). Положение точки \mathbf{x}_0 выбрано в основании вектора \boldsymbol{\omega}, поэтому \boldsymbol{\omega} и \mathbf{x}-\mathbf{x}_0 коллинеарны (лежат на одной прямой). Если вектор \boldsymbol{\omega} единичный (\boldsymbol{\omega}^2=1), то скалярное произведение векторов \mathbf{x}-\mathbf{x}_0 и \boldsymbol{\omega} равно расстоянию точки \mathbf{x} до плоскости:

d = (\mathbf{x}-\mathbf{x}_0)\cdot \boldsymbol{\omega} ~=~ -\mathbf{x}_0\cdot\boldsymbol{\omega} + \mathbf{x}\cdot\boldsymbol{\omega} ~=~ b + \mathbf{x}\cdot\boldsymbol{\omega}.

Если длина \omega=|\boldsymbol{\omega}| вектора \boldsymbol{\omega} нормали к гиперплоскости отлична от единицы, то d в \omega раз больше (\omega > 1) или меньше (\omega < 1) евклидового расстояния в n-мерном пространстве. Когда векторы \mathbf{x}-\mathbf{x}_0 и \boldsymbol{\omega} направлены в противоположные стороны: d < 0. ►

Если пространство имеет 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генерим сначала модельные данные в которых выходы линейно связаны с входами и добавлен небольшой случайный шум:

import numpy as np                                    # работа с тензорами
 
N      = 100                                          # число точек  (примеров)
nX, nY = 2, 3                                         # число входов (признаков) и выходов
 
w = np.array([[  135],                          # матрица w формы (nX,nY)
              642] ])
b = np.array( [ -3, -2, -1]  )                        # вектор  b формы (nY, )
 
= np.random.random ((N, nX))                       # равномерные случ.точки в [0...1]
= np.dot(X, w) + b                                 # выходы модели
Y += np.random.normal (0, 0.01, (N, nY))              # добавляем случайный гаусовый шум

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

from numpy.linalg import inv, pinv                    # обратные для квадрат. и прямоуг.
 
XX = np.concatenate( (X, np.ones((N, 1))), axis=1)    # добавляем столбик единиц
 
W1 = np.dot( inv(np.dot(XX.T, XX)), np.dot(XX.T, Y) ) # первая формула
W2 = np.dot( pinv(XX), Y)                             # вторая формула


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

Для решения этой-же задачи, на практике удобнее пользоваться готовой библиотекой машинного обучения Scikit-learn (кратко sklearn):
from sklearn.linear_model import LinearRegression
 
lr = LinearRegression().fit(X,Y)  
 
print(lr.coef_.T)                   #  [[1.007  3.001  5.001]  матрица W
                                    #   [6.003  4.     1.998]]
print(lr.intercept_)                #  [-3.002 -2.001 -0.999]  вектор b

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

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

mse =  np.mean((X @ lr.coef_.T + lr.intercept_ - Y) ** 2)
 
print(f"mse:{mse:.5f}, sqrt(mse):{np.sqrt(mse):.5f} ")
#       mse:0.00010,   sqrt(mse):0.01005 


Классификация на 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} до трёх плоскостей. Расстояния положительны, если вектор нормали "смотрит" на точку и отрицательны - если в противоположную сторону. Метод наименьших квадратов будет стремится развернуть эти векторы, чтобы, по крайней мере, совпадали знаки выхода модели и целевых значений. Это даст необходимый классификатор.

Подготовим модельные данные в виде трёх кластеров с нормальным разбросом точек вокруг их центров:
nx, ny, N = 2, 3, 300                    # число признаков, классов, объектов (точек)
 
C = np.zeros(N)                          # номера правильных классов
Y = np.full((N,ny), -1)                  # матрица N x ny для выходов, заполненная -1
X = np.random.normal (0, 0.4, (N,nx))    # облако N точек с центром в начале координат
 
X0 = [[-1,1], [1,1], [0,-1]]             # центры кластеров каждого класса
nc = int(N/ny)                           # число примеров каждого класса
 
for i, x0 in enumerate(X0):
    C[i*nc:(i+1)*nc]    = i              # номер i-го класса
    Y[i*nc:(i+1)*nc, i] = 1              # целевые выходы
    X[i*nc:(i+1)*nc]   += np.array(x0)   # сдвигаем кластер i-го класса

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

class Model:
    def __init__(self):
        "Конструктор"
        self.lr = LinearRegression()
         
    def __call__(self, x):
        "Функциональный вызов объекта"
        y = x @ self.lr.coef_.T + self.lr.intercept_
        return np.argmax(y, axis=1)
     
    def fit(self, x,y):
        "Обучение по примерам x,y"
        self.lr.fit(x,y)

Процесс обучения и тестирования выглядит следующим образом:
model = Model()       # экземпляр модели
 
model.fit(X,Y)        # обучаем
 
y = model(X)          # поучаем номера классов
A = (y==C).mean()     # сравниваем с точными (аккуратность)

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

.   def plot(self, x1_min=0, x1_max=1, x2_min=0, x2_max=1, n=201):
        # сетка точек 201 x 201 в интервале x1,x2=[xMin,xMax]
        x1,x2 = np.meshgrid(np.linspace(x1_min, x1_max, n),
                            np.linspace(x2_min, x2_max, n))
        grid  = np.c_[x1.ravel(), x2.ravel()]     # shape = (n*n, 2)       
             
        yp = self(grid).reshape(x1.shape)         # (n*n, ) -> (n, n)
         
        fig,ax = plt.subplots(figsize=(5, 5))     # размеры картинки (квадрат)       
        plt.axis([x1_min,x1_max, x2_min,x2_max])  # диапазон изменения осей   
         
        plt.imshow(yp, interpolation='nearest',   # цвет. карту для 2D массива yp
                  extent=(x1.min(), x1.max(), x2.min(), x2.max()),
                  aspect='auto', origin='lower', alpha=0.2,
                  vmin=0, vmax=yp.max())       

С помощью этого метода рисуем разделяющие поверхности и обучающие объекты:
model.plot(-2, 2-2, 2)
plt.scatter(X[:,0], X[:,1], c=C, s=30, edgecolors='black')
plt.show()

Введение в машинное обучение <
> Градиентный метод