16 Apr 2022
Алгоритмы и структуры данных
http://mech.math.msu.su/~vvb/
public_html/index.html
Мы будем использовать Python3.
Про Питон: это скриптовый язык
Сам по себе Питон -- это объектно-ориентированный язык, значит,
очень медленный, в сотни раз медленнее, чем С, С++, Fortran
Понятно, что Питон задумывался как игрушечный язык (Гвидо Россумом).
У него нашлось много последователей, они называются питонистами.
Язык стал настолько популярен, что его стали использовать
в реальных приложениях. Питон позволяет подключать модули, написанные
на другом языке -- Fortran, C, C++. В этих модулях все операции
выполняются быстро.
На самом Питоне пишется только верхний уровень программы.
Например, мы используем модуль numpy (самый популярный модуль в Питоне),
в нем мы задаем матрицы. При этом сама команда, например, умножения матриц,
вычисления обратной матрицы, решения системы линейных уравнений и т.п.
выполняется уже внутри модуля numpy. Аналогично вся работа с нейросетями
выполняется внутри модуля tensorflow, и т.п.
Стиль программирования на Питоне: как правило, большинство алгоритмов
программисты не реализуют самостоятельно, а используются модули, в которых эти
алгоритмы уже реализованы.
Пример: решение системы линейных уравнений или нахождение
обратной матрицы. Вся работа с матрицами организована внутри модуля
numpy, алгоритмы линейной алгебры реализованы внутри подмодуля
numpy.linalg (от слова линейная алгебра).
Как установить Питон на свой компьютер?
1. Установить сам по себе Питон.
2. Альтернатива -- можно использовать Anaconda для установки Python
и разных модулей для Питона.
3. Есть удобная оболочка, которая раньше называлась
Jupiter Notebook (этот проект завершен, к настоящему
моменту устарел), теперь называется Jupiter-lab. Это
оболочка, которая объединяет текстовый редактор и интерпретатор
Питона. Все это работает под браузером.
ПЕРЕРЫВ 10 минут до 15:00
Moscow
Мы реализовали алгоритм Гаусса приведения матрицы к ступенчатому
виду.
Время работы алгоритма Гаусса O(n^3), где n -- размер матрицы.
Если вычислять обратную матрицу с помощью алгебраических дополнений
и вычисления определителя, то время работы будет O(n^5).
Мораль: этот алгоритм неприемлем!
Обратная матрица вычисляется с помощью сначала алгоритма Гаусса
приведения матрицы к ступенчатому виду, а затем обратного прохода.
Рассмотрим вычисление обратной матрицы.
Приписываем справа к исходной матрице единичную матрицу
1 2 3
4 5 6 --> приписываем справа единичную матрицу
7 8 0
1 2 3 | 1 0 0
4 5 6 | 0 1 0 --> приводим расширенную матрицу к ступенчатому виду
7 8 0 | 0 0 1
* * * | * * *
0 * * | * * * --> домножаем строки на элементы, обратные к
0 0 * | * * * диагональным
1 * * | * * *
0 1 * | * * * --> выполняем обратный проход
0 0 1 | * * *
1 0 0 | * * * в правой части мы получаем обратную матрицу!
0 1 0 | * * *
0 0 1 | * * *
A --> E
ek ... e3 e2 e1 A = E =>
ek ... e3 e2 e1 = inv(A)
Домашнее задание
Студенты с нечетными номерами делают задачу 1,
с четными 2
1. Реализовать вычисление обратной матрицы (для numpy-матриц):
def invMatrix(a, eps = 1e-8):
# Матрица a не изменяется
b = a.copy()
. . .
return c
2. Решить систему линейных уравнений
def solveLinearSystem(a, b, eps = 1e-8):
# a -- квадратная невырожденная матрица,
# b -- линейный массив
# исходные данные не должны меняться
# Return: массив x решения системы
# a x = b
Решения присылайте на почту
vladimir_borisen@mail.ru
В теме письма (поле Subject:)
записывайте Душанбе
================================================
19 Apr 2022
Метод наименьших квадратов в задаче линейной
регрессии. Решение этой задачи с помощью псевдообратной
матрицы Мура-Пенроуза
Задача линейной регрессии.
Пусть у нас есть некоторая независимая величина x <- R^n
x = (x1, x2, ..., x_n) и зависящая от нее величина y <- R
Пусть зависимость описывается линейной функцией
y = <w, x> + b, w <- R^n
(угловые скобки означают скалярное произведение). Через координаты
y = w1*x1 + w2*x2 + ... + w_n*x_n + b
Мы рассматриваем задачу машинного обучения: параметры зависимости w, b
нам неизвестны, мы хотим восстановить их по заданной выборке:
v1, v2, ..., v_m
y1, y2, ..., y_m y_i = <w, v_i> + b
Пусть m > n, мы получаем линейную систему, в которой число уравнений
больше, чем число неизвестных:
<w, v1> = y1
<w, v2> = y2
. . .
<w, vm> = ym
Если записать через координаты.
Пусть
v_i = (x_i1, x_i2, ..., x_in)
Тогда система записывается следующим образом:
w1*x_11 + w2*x_12 + ... + w_n*x_1n + b = y1
w1*x_21 + w2*x_22 + ... + w_n*x_2n + b = y2
. . .
w1*x_m1 + w2*x_m2 + ... + w_n*x_mn + b = y_m
Здесь w_i -- неизвестные, а элементы матрицы -- это x_ij
[x_11 x_12 ... x_1n]
A = [x_21 x_22 ... x_2n]
[. . . ]
[x_m1 x_m2 ... x_mn]
A -- это матрица размера mxn, где m >= n (число строк не меньше,
чем число столбцов).
Система уравнений записывается как
AV + B = Y B = (b, b, ..., b)^T, Y = (y1, y2, ..., y_m)
Вектор B -- это свободный член. Обычно от него избавляются с помощью следующего
приема:
к каждому объекту x = (x1, x2, ..., x_n) мы добавляем еще одну координату,
тождественно равную единице
x' = (x1, x2, ..., x_n, 1)
v'_i = (x_i1, x_i2, ..., x_in, 1)
Тогда для этих "расширенных" объектов зависимость запиисывается
совсем просто:
y = <w', x'> w' = (w', b)
Мы избавились от свободного члена. Опустим штрих, чтобы не загромождать
обозначения:
y = <w, x>
Система уравнений имеет вид:
AX = Y
где A -- матрица размера mxn, m >= n
Вообще говоря, эта система не имеет решения.
В методе наименьших кадратов мы минимизируем суммарную квадратичную
ошибку. Мы ищем вектор W такой, что
||AW - Y||^2 --> min
Сумма_i=1^m |<v_i, w> - y_i|^2 --> min
A: R^n --> R^m, m >= n
Y <- R^m
Мы ищем такой вектор W, что AW максимально близко к Y ==>
вектор AW есть проекция Y на подпространство Im(A) = A R^n ==>
это значит, что разность векторов Y - AW перпендикулярна
пространству Im(A) ==> для всякого вектора X скалярное
произведение равно нулю:
<AW - Y, AX> = 0 для всякого вектора X
В матричном виде
(AW - Y)^T AX = 0
((AW - Y)^T A) X = 0
Поскольку это верто для любого X, то матрица слева равна нулю:
(AW - Y)^T A = 0 (U V)^T = V^T U^T
(AW - Y)* A = 0
(W* A* - Y*) A = 0
W* A* A = Y* A
Транспонируем матрицы в обеих частях:
A* A W = A* Y
Мы получили так называемое нормальное уравнение в методе наим. квадратов.
В линейной алгебре справедливо следующее утверждение:
пусть мы имеем прямоугольную матрицу, у которой столбцы линейно независимы.
Тогда матрица A* A обратима (это квадратная матрица; ее элементы --
это попарные скалярные произведения векторов-столбцов матрицы A).
Домножив равенство на обратную матрицу inv(A* A), получим
W = (inv(A* A) A*) Y
Итак, мы нашли решение нашей задачи -- вектор W
Матрица
(A* A)^(-1) A* называется псевдообратной матрицей (Мура-Пенроуза)
к матрице A.
Для произвольной прямоугольной матрицы псевдообратная матрица
вычисляется через СИНГУЛЯРНОЕ РАЗЛОЖЕНИЕ произвольной матрицы
(SVD -- Singular Value Decomposition):
A = U D V
U, V -- ортогональные квадратные матрицы
U* U = 1, V* V = 1
D -- диагональная (прямоугольная матрица)
Пример: 1 0 0 0 0 0
0 2 0 0 0 0
0 0 5 0 0 0
Пусть A = U D V -- сингулярное разложение, тогда псевдообратная матрица
pinv(A) = V* D' U*
D' получается из D
1) транспонированием;
2) заменой ненулевых диагональных элементов на обратные.
В Питоне всё это реализовано в модуле numpy.linalg
Почему задача нахождения линейной зависимости называется регрессией?
В XIX веке ученый исследовал зависимость между ростом отца и
ростом его взрослого сына.
По набору точек (рост отца, рост сына) получилась прямая,
угол наклона которой меньше 45 градусов (и которая не
проходит через точку (0, 0), а пересекает ось y в точке (0, y) где y > 0).
В среднем получилось, что у высокого отца сын ниже него, у низкого -- выше,
т.е. как бы рост по мере смены поколений стремится к среднему значению =>
отсюда возник термин "регрессия".
=====================================================
20 Apr 2022
Метод сопряженных градиентов
Имеется квадратичная функция от n переменных
f(x) = 1/2 <x, Ax> + <b, x> --> min
x <- R^n, b <- R^n
Ищем минимум функции по x
A -- симметричная положительно определенная матрица
A* = A, для всякого x != 0 <x, Ax> > 0
У нас есть обычное скалярное произведение <x, y>
Свойства:
1) линейность по аргументам <x1 + x2, y> = <x1, y> + <x2, y>
<a x, y> = a <x, y>
2) симметрия <x, y> = <y, x>
3) положительная определенность:
для любого x != 0 <x, x> > 0 (квадрат длины вектора положителен).
Любое скалярное произведение, удовлетворяющее этим свойствам,
задается матрицей грамма:
a_ij = <e_i, e_j>
A = [ a_ij ]
Тогда скалярное произведение, заданное матрицей грамма A,
вычисляется с помощью формулы
<x, y>_A = <x, Ay> = <Ax, y>
Два вектора называются сопряженными относительно матрица A (симметричной и
положительно определенной), если <x, Ay> = 0 (аналог перпендикулярности
векторов).
Вернемся к нашей задаче:
f(x) = 1/2 <x, Ax> + <b, x> --> min
Минимум единственный. В точке минимума градиент функции равен нулю.
Выпишем формулу для градиента функции f:
f'(x) = A x + b
Обозначим x* -- точка минимума.
f'(x*) = 0
A x* + b = 0
x* -- решение линейного уравнения A x* = -b
Нахождение точки минимума эквивалентно решению системы линейных уравнений
с симметричной положительно определенной матрицей A.
Получаем ответ:
x* = - A^(-1) b
Почему мы на этом не останавливаемся???
Потому что обращение матрицы -- довольно сложная вычислительная,
причем она неустойчива, если мера обусловленности матрицы большая.
Мера обусловленности матрицы -- это
|A| |A^(-1)|
Норма матрицы |A| -- это во сколько раз она может удлинить вектор
max_v |Av|/|v|
Для хороших матриц это есть
максимальное собственное значение/минимальное собственное значение
Итак, рассмотрим метод сопряженных градиентов для поиска минимума
квадратичной функции
f(x) = 1/2 <x, Ax> + <b, x>
Основная идея метода:
возьмем произвольное начальное приближение x0.
Пусть имеется набор из n сопряженных векторов s1, s2, ..., s_n.
<s_i, A s_j> = 0 при i != j.
Легко доказать, что векторы s_i линейно независимы.
Мы ищем минимум нашей функции в виде
x* = x0 + a1*s1 + a2*s2 + ... + a_n*s_n,
где коэффициенты a_i -- это вещественные числа
Как выражаются коэффициенты a_i?
Мы последовательно вычисляем точки
x1 = x0 + a1*s1
x2 = x1 + a2*s2
. . .
x_i = x_i-1 + a_i*s_i
. . .
x_n = x_n-1 + a_n*s_n = x*
Коэффициент a_i находится из задачи одномерной минимизации:
x_i = x_i-1 + a_i*s_i
f(x_i) --> min по переменной a_i (ищем минимум функции f на прямой
x_i-1 + a_i*s_i)
d/da f(x_i-1 + a*s_i) = 0
Производная функции от многих переменных по некоторому направлению
равна скалярному произведению градиента функции на направление
Grad f(x) = A x + b
<Grad f(x_i-1 + a_i*s_i), s_i> = 0
Отсюда вычисляем a_i:
<A(x_i-1 + a_i*s_i) + b, s_i> =
= <Ax_i-1, s_i> + a_i <A s_i, s_i> + <b, s_i> = 0
a_i = (-<A x_i-1, s_i> - <b, s_i>) / <A s_i, s_i>
a_i = <-f'(x_i-1), s_i>/<A s_i, s_i>
Мораль:
если мы имеем систему сопряженных векторов
s1, s2, ..., s_n,
то мы в точности за n шагов находим точку минимума функции f.
Откуда взять эти направления?
Мы их тоже вычисляем последовательно.
1) Начинаем с вектора s1 = -f'(x0)
Вычисляем коэффициент a1 из условия минимума функции f(x) на прямой
x0 + a1*s1
Получаем точку x1 = x0 + a1*s1, где a1 = <-f'(x0), s1>/<A s1, s1>
2) i-й шаг. Мы уже вычислили x_i-1 и направление s_i-1
Вычисляем направление s_i из условия того, что <s_i, A s_i-1> = 0
(т.е. s_i сопряжено с s_i-1). На самом деле, получится и то,
что s_i будет сопряжено со всеми предыдущими векторами s_j, j < i.
Ищем s_i в виде
s_i = -Grad f(x_i-1) + beta*s_i-1
0 = <s_i, A s_i-1> = <-f'(x_i-1) + beta*s_i-1, A s_i-1> =
= <-f'(x_i-1), A s_i-1> + beta<s_i-1, A s_i-1>
Отсюда выражаем beta
beta = <-f'(x_i-1), A s_i-1> / <s_i-1, A s_i-1>
И дальше вычисляем коэффициент a_i по той формуле, которую мы уже
вывели раньше:
a_i = <-f'(x_i-1), s_i>/<A s_i, s_i>
x_i = x_i-1 + a_i s_i =
= x_i-1 - <f'(x_i-1), s_i>/<A s_i, s_i> * s_i
Выпишем алгоритм (метод Хестенса-Штифеля):
x0 -- начальное приближение
цикл для k = 1, 2, ..., n выполнять
если f'(x_k-1) = 0 то ответ = x_k-1; завершить алгоритм
если k = 1 (первая итерация), то
s1 = -f'(x0)
иначе если k > 1, то
beta_k-1 = <f'(x_k-1>, A s_k-1> / <s_k-1, A s_k-1>
s_k = -f'(x_k-1) + beta_k-1*s_k-1
конец если
x_k = x_k-1 - ( <f'(x_k-1), s_k>/<s_k, A s_k> )*s_k
конец цикла
ПЕРЕРЫВ 10 минут до 13:55 (Москва)
Метод Хестенса-Штифеля применим только к квадратичной функции.
А как быть с произвольной функцией? Ее можно разложить в ряд Тейлора и
в качестве матрицы A рассматривать Гессиан (матрицу вторых производных).
Однако вычисление Гессиана в каждой новой точке -- очень затратная операция.
Идея: давайте в квадратичном случае преобразуем наши формулы так,
чтобы в них не участвовала матрица A.
Я не буду выписывать преобразования, а только укажем результат:
beta_k = <f'(x_k>, A s_k> / <s_k, A s_k> = . . . =
= ||f'(x_k)||^2 / ||f'(x_k-1)||^2
Это формула Флетчера-Ривса
s_k+1 = -f'(x_k) + beta_k*s_k
Получили формулы для сопряженных направлений, в которых матрица гессиана A
не участвует! Коэффициент a_k+1 вычисляется как результат ОДНОМЕРНОЙ
минимизации функции f(x) на прямой x = x_k + a_k+1*s_k+1
a_k+1 = argmin_a f(x_k + a*s_k+1)
Эти формулы в точности соответствуют методу Хестенса-Штифеля для
квадратичной функции, но они уже применимы к произвольной выпуклой
функции. Мы получаем метод минимизации Флетчера-Ривса.
Это итеративный метод, но, в отличие от метода Хестенса-Штифеля,
он не находит сразу за n шагов точный минимум.
x0 ---> n шагов метода Флетчера-Ривса ---> x_n1
x_n1 ---> n шагов метода Флетчера-Ривса ---> x_n2
. . .
Этот процесс сходится очень быстро, для хороших функций этот метод
за n итераций получает точность порядка O(1/n^2).
========================================================================
Домашнее задание 2
Реализовать графическую оконную функцию, иллюстрирующую
метод наименьших квадратов.
Пользователь задает точки кликами мыши:
p0, p1, p2, ..., p_m-1
С помощью слайдера задается задается степень многочлена d.
Надо вычислить методом наименьших квадратов многочлен f степени d
такой, что сумма квадратов отклонений значений многочлена
в точке p.x от значения p.y была бы минимальной:
Sum_i (f(p[i].x) - p[i].y)^2 --> min
Минимум находится по коэффициентам многочлена (w0, w1, ..., w_d):
f(x) = w0 x^d + w1 x^d-1 + ... + w_d-1 x + w_d
f(x) = <w, [x^d, x^d-1, ..., x, 1]>
import numpy as np
def computeLeastSquaresPol(x, y, degree):
m = len(x)
assert(m == len(y))
a = np.vander(x, degree + 1)
pinv_a = np.linalg.pinv(a)
w = pinv_a @ y
return w
Этот код надо вставить в оконное приложение на Python3
с использованием модуля tkinter.
В качестве образца дается программа plotFunc.py
http://mech.math.msu.su/~vvb/
============================================
23 Apr 2022
Метод опорных векторов
Метод опорных векторов (support vectors machine, SVM)
относится к задаче классификации объектов.
Пример: объекты -- это потенциальные клиенты банка
Каждый объект описывается набором признаков (features)
Признаки -- в простейшем случае просто числа.
Например, для клиента признаки:
1) средний доход в месяц;
2) возраст;
3) сумма невыплаченных на настоящее время кредитов;
4) месячная суммы выплат по текущим кредитам;
5) суммарная стоимость собственности: автомобили,
квартиры и т.п.;
6) семейное положение (число близких родственников);
7) число детей.
Объект описывается точкой 7-мерного пространства R^7
В методе опорных векторов строится линейный классификатор.
x -- объект, x <- R^n
Классификатор -- это линейная функция
f(x) = <w, x> - b, w <- R^n, b <- R
Если f(x) > 0, то клиент надежный, ему можно выдать кредит.
Если f(x) < 0, то клиент ненадежный, кредит выдавать ему нельзя.
Если f(x) = 0, то непонятно.
Уравнение f(x) = 0 определяет гиперплоскость в n-мерном пространстве,
которая разделяет два класса объектов (отделяет надежных клиентов от
ненадежных).
Когда мы составляем тренировочную выборку, то нам удобно сделать
разметку симметричной: отмечаем надежных клиентов значением y = 1,
ненадежных -- значением -1
В общем случае в тренировочной выборке есть точки двух классов,
точки первого класса отмечаются значением y = 1,
точки второго класса -- значением y = -1
x1, x2, x3, ..., x_m <- R^n
y1, y2, y3, ..., y_m <- { -1, 1 }
Мы строим линейный классификатор
f(x) = <w, x> - b
который разделяет точки x_i двух классов
y_i = 1 ==> f(x_i) > 0
y_i = -1 ==> f(x_i) < 0
Гиперплоскость f(x) = 0 должна разделять эти два множества точек.
Она задается уравнением
<w, x> - b = 0
Отметим, что классификатор дает правильный ответ на объекте x_i <==>
y_i (<w, x_i> - b) > 0
Когда ответ неправильный, то y_i и f(x_i) имеют разные знаки.
Рассмотрим сначала разделимый случай: тогда существует
гиперплоскость f(x) = <w, x> - b = 0 такая, что
все точки первого класса из нашей выборки лежат в верхнем полупространстве
<w, x_i> - b > 0, y_i > 0
y_i(<w, x_i> - b) > 0
а точки второго класса -- в нижнем полупространстве:
<w, x_i> - b < 0, y_i < 0
y_i(<w, x_i> - b) < 0
Формула для расстояния со знаком от точки x до гиперплоскости <w, x> - b = 0
h = (<w, p> - b)/|w|
Мы хотим, чтобы расстояние h1 от точек первого класса до гиперплоскости
равнялось расстоянию h2 от точек второго класса h = h1 = h2 и было бы
максимальным:
y_i (<w, x_i> - b)/|w| >= h для всех i=1, 2, ..., m
h --> max
Разделим обе части равенства на h
y_i (<w, x_i> - b)/(|w|*h) >= 1 для всех i
Обозначим w/(|w|*h) = w'
y_i (<w', x_i> - b) >= 1 для всех i
Выразим h через длину |w'|
w/(|w|*h) = w'
|w'| = |w|/(|w|*h) = 1/h
Получаем задачу максимизации
y_i (<w', x_i> - b) >= 1 для всех i
h --> max <==>
1/h --> min <==>
|w'| --> min
Давайте опустим штрих для упрощения. Получаем задачу:
y_i (<w, x_i> - b) >= 1 для всех i=1,...,m
|w|^2 --> min_w,b
Это задача квадратичного программирования: минимизировать квадратичную
функцию при заданном наборе ограничений типа равенств
или неравенств.
ПЕРЕРЫВ 10 минут до 13:42 Москва
Рассмотрим неразделимый случай, когда задача
y_i (<w, x_i> - b) >= 1 для всех i=1,...,m
|w|^2 --> min
неразрешима.
Разрешим неравенствам y_i (<w, x_i> - b) >= 1 выполняться не строго,
с некоторой ошибкой, и обозначим величину этой ошибки через xi_i >= 0
(переменные xi_i называются расслабляющими переменными slack variables):
y_i (<w, x_i> - b) >= 1 - xi_i
Но тогда нам надо минимизировать общую сумму ошибок Sum xi_i
В реузльтате получаем следующую задачу миннимизации:
y_i (<w, x_i> - b) >= 1 - xi_i, i=1,...,m
xi_i >= 0,
|w|^2 + C Sum xi_i --> min
Здесь C -- гиперпараметр нашего алгоритма. Чем C больше, тем большее
значение имеет минимизация суммы ошибок xi_i. При малых C большее
значение ширина разделяющей полосы между классами.
это также задача квадратичного программирования, которая решается
методами условной минимизации (например, по теореме Куна-Такера).
Попробуем избавится от переменных xi_i и свести эту задачу к задаче
безусловной минимизации.
Для каждой переменной xi_i мы имеем два неравенства:
xi_i >= 0
y_i (<w, x_i> - b) >= 1 - xi_i
Преобразуем второе неравенство:
xi_i >= 0
xi_i >= 1 - (y_i (<w, x_i> - b))
Поскольку xi_i -- это плата за ошибку, мы заинтересованы в том
чтобы сделать xi_i как можно меньше. То есть xi_i должно равняться
максимуму из правых частей этих двух неравенств:
xi_i = max(1 - (y_i (<w, x_i> - b)), 0)
Обычно рассмотривают функцию
{ 1 - x, x <= 1
HingeLoss(x) = {
{ 0, x > 1
Получаем xi_i = HingeLoss(y_i (<w, x_i> - b))
Таким образом, мы избавились от переменных xi_i и ограничений xi_i >= 0,
xi_i >= 1 - (y_i (<w, x_i> - b)). Наша задача переписывается следующим
образом:
|w|^2 + C Sum xi_i --> min
подставляем выражение xi_i = HingeLoss(y_i (<w, x_i> - b))
Получаем задачу безусловной минимизации:
|w|^2 + C*Sum HL(y_i*<w, x_i> - b) --> min_{w,b}
где через HL мы обозначили функцию HingeLoss.
Давайте продемонстрируем это на примере в Питоне.
В Питоне есть алгоритм минимизации чего угодно
from scipy import minimize
res = minimize(f, x0) # f = f(x) -- функция, минимум которой мы ищем
# x0 -- начальное приближение
======================================================
26 Apr 2022
В линейном методе опорных векторов далеко не всегда можно построить
гиперплоскость, удовлетворительно разделющая два класса точек в
n-мерном пространстве. Как быть???
Идея: дополним набор признаков (features) объектов еще несколькими
признаками, которые вычисляются как функции от исходных признаков.
Пример: точки первого класса расположены вблизи начала координат,
точки второго класса окружают точки первого класса и расположены
дальше от начала координат.
Добавим признак "квадрат расстояния до начала координат"
x == (x0, x1) --> x' == (x0, x1, z), где z = x0*x0 + x1*x1
x <- R^2 x' <- R^3
Тогда уже в трехмерном пространстве точки успешно разделяются гиперплоскостью,
поскольку координата z для точек второго класса больше, чем для первого.
В общем случае признаки объекта
x = (x0, x1, ..., x_n-1) можно дополнить всеми мономами от x_i\
степени <= d
Например, пусть x имеет 2 признака, x <- R^2
x = (x0, x1)
Пусть степень мономов d = 2. Тогда мы дополняем x следующими признаками:
x' = (1, x0, x1, x0^2, x0*x1, x1^2) <- R^6
=======================================================
27 Apr 2022
1. Как сгенерировать случайный набор точек на плоскости?
import random
x = random.normalvariate(mean, sigma)
Если нужно сгенерировать случайные точки по кругу, то
можно использовать полярную систему координат:
p = (r*cos(phi), r*sin(phi)),
где и r, и phi также распределены по нормальному закону.
p = ex * (r*cons(phi)) + ey * (r*sin(phi))
2. Пусть есть функция, определенная для точек плоскости
z = f(x, y)
Надо нарисоать линию уровня этой функции, заданную уравнением
f(x, y) = l = const
(l от слова level). Такие линии называются изолиниями или
изогипсами.
Воспользуемся модулем skimage и в нем подмодулем measure
Чтобы вычислить изолинии, соответствующие уравнению
f(x, y) = l,
надо заполнить матрицу
a_ij = f(x_j, y_i)
размера mxn (m строк, n столбцов),
где
x_j = xmin + dx*j, dx = (xmax - xmin)/n
y_i = ymin + dy*i, dy = (ymax - ymin)/m
Затем надо вызвать функцию measure.find_contours
contours = measure.find_contours(a, l)
Массив contours состоит из элементов типа массив точек,
где каждая точка -- это пара (номер строки, номер столбцаЗ)
где номера представлены вещественными числами (не обязательно целые!).
Чтобы получить координаты точек контура, надо проделать следующие
вычисления:
def xcoordinate(j):
dx = (xmax - xmin)/n
x = xmin + j*dx
return x
def ycoordinate(i):
dy = (ymax - ymin)/m
y = ymin + i*dy
return y
ПЕРЕРЫВ 10 минут до 13:48 (Moscow time)
15:48 Dushanbe
Построение классификатора.
Итак, мы имеем массив объектов x размера m, каждый объект описывается
точкой плоскости (число признаков = 2)
x = [ x[0], x[1], ..., x[m-1] ],
x[i] -- точка плоскости: ( x[i, 0], x[i, 1] )
Разметка нашей выборки:
y = [ y[0], y[1], ..., y[m-1] ],
y[i] = 1, если объект x[i] принадлежат 1-му классу,
= -1, если объект x[i] принадлежат 2-му классу.
Для построения классификатора w, b надо минимизировать
функцию потерь. Но мы вычисляем гиперплоскость, разделяющую 2 класса,
в расширенном пространстве: мы добавляем к признака (x0, x1) объекта x
всевозможные произведения признаков объекта x степени не выше degree.
Например, при degree = 2 мы получаем расширенный объект
xext = (1, x0, x1, x0^2, x0*x1, x1^2) <- R^6
При degree = 6 получаем
xext = (1, x0, x1, x0^2, x0*x1, x1^2, x0^3, x0^2*x1, x0*x1^2, x1^3) <- R^10
Обозначим через n размерность расширенного пространства.
В расширенном пространстве гиперплоскость задается двумя параметрами:
вектор нормали
w = (w[0], w[1], ..., w[n-1]) <- R^n
свободный член
b <- R
Гиперплоскость задается уравнением
<w, xext> - b = 0
Значение классификатора на векторе x исходного пространства:
xext = monomials(x, degree)
c = <w, xext> - b
Параметры классификатора находятся путем минимизации функции ошибки:
wb = (w, b)
lossFunction(wb) =
= |w|^2 + C*1/m*Sum HingeLoss( y[i]*(<w, xext[i]> - b) )
--> min_w,b
Как выполняется минимизация в Питоне?
res = minimize(lossFunction, wb0)
где wb0 -- начальное значение параметров wb0 = (w0, b0)
Точка минимума будет вычислена в переменной-члене x:
wb = res.x
============================================
28 Apr 2022
Ядровый метод опорных векторов (Kernel Support Vector Machine).
Теорема Карруша--Куна--Такера о нахождении условного минимума
функции от многих переменных при ограничениях типа равенств и
неравенств.
Вспомним, как мы формулировали задачу в методе опорных векторов.
В результате получаем следующую задачу миннимизации:
y_i (<w, x_i> - b) >= 1 - xi_i, i=1,...,m
xi_i >= 0,
|w|^2 + C Sum xi_i --> min
При добавлении новых признаков (features) объекта мы
как бы определяем отображение
phi: x --> x' <- R^k, k > n
В расширенном пространстве мы строим классификатор
f'(x') = <W, x'> - B
Получается классификатор для исходного пространства
f(x) = <W, phi(x)> - B
Обратите внимание, здесь мы вычисляем скалярное произведение в расширенном
пространстве R^k, k > n.
Определение (из функционального анализа):
функция K(x, y) = s <- R, x, y <- R^n называется ядром,
если существует отображение
phi: x --> x', x' <- H, где H -- гильбертово пространство
такое, что
K(x, y) = < phi(x), phi(y) >
В машинном обучении / искусственном интеллекте используются в основном
два ядра:
1) Гауссовское ядро:
rbf(x, y) = 1/(2 pi sigma) exp( -||x - y||^2 / (2*sigma))
Radial Base Function;
Здесь спрямляющее пространство бесконечномерно;
2) полиномиальное ядро степени d >= 1 (d -- целое число):
poly(x, y) = (<x, y> + rho)^d, rho <- R.
ПЕРЕРЫВ 14 минут до 13:35 (Miscow time)
15:35 (Dushanbe)
Для использования ядрового метода опорных векторов импортируем класс
SVC из модуля sklearn.svm:
from sklearn.svm import SVC
Затем создаем классификатор:
classifier = SVC(C=2., kernel="rbf")
Здесь мы создали классификатор с гиперпараметром C = 2.,
который будет использовать Гауссовское ядро rbf (Radial Basis Function):
K(x1, x2) = exp( -||x1 - x2||^2/(2*sigma^2) )
Тренируем классификатор на обучающей выборке:
classifier.fit(x, y)
Здесь x -- массив векторов, представляющих наши объекты,
y -- разметка объектов: y[i] == 1, если x[i] принадлежит первому классу,
y[i] == -1, если x[i] принадлежит второму классу.
Обычно x представляется как numpy-матрица. Например,
x = np.array([
[1., 2.],
[3., 4.],
[0., -1.],
[2.5, 0.3]
])
y = np.array([
1.,
-1.,
-1.,
1.
])
Здесь мы задали 4 объекта, первый и последний объекты принадлежат
первому классу, второй и третий -- второму классу.
После этого мы может применить классификатор к любому объекту xxx,
чтобу определить (предсказать) его класс:
с = classifier.predict(xxx)
Результатом будет с = +1 или c = -1.
Если мы хотим определить точное значение расстояния со знаком
до разделяющей гиперплоскости, то лучше использовать метод decision_function:
z = classifier.decision_function(xxx)
Число z принимает любые вещественные значения. В частности, z=0
на линии, разделяющей два наших класса.
Домашнее задание (deadline: вторник 3 мая)
Вариант 1.
Переделать программу svm.py, где реализован классический
линейный метод опорных векторов. Надо дополнить признаки объектов
до всех мономов (произведений) степени <= d от исходных признаков,
d должно задаваться слайдером. Например, при d = 2 мы
вместо точки x = (x0, x1) <- R^2 используем "расширенную" точку
x' = (1, x0, x1, x0^2, x0*x1, x1^2) <- R^6
при d=3 получится 10-мерное пространство и т.д.
Линейный классификатор вычисляется для расширенного пространства.
Программа должна нарисовать линию, разделяющую точки двух классов.
Вариант 2.
Переделать программу svm.py, где реализован классический
линейный метод опорных векторов. Модифицированная программа
должна использовать ядровый метод опорных векторов, реализованный
в стандартных модулях Питона в виде класса SVC:
from sklearn.svm import SVC
Программа должна применить классификатор SVC и нарисовать линию,
разделяющую точки двух классов.
=============================================================
30 Apr 2022
Вернемся к задаче линейной регрессии.
Напомним формулировку задачи
Надо построить линейную функцию, хорошо аппроксимирующую зависимость
между двумя величинами:
x = (x1, x2, ..., x_n) <- R^n -- независима величина.
y = y(x) <- R -- зависимая от x величина.
Априори известно, что зависимость описывается оинейной функцией
y = <w, x> + b
w <- R^n вектор весов, b <- R -- свободный член
(intercept)
У нас имеется тренировная выборка, для которой известно значение
зависимой величины y
x1, x2, x3, ..., x_m -- для этих значений x измерено значение y:
y1, y2, y3, ..., y_m
Отметим, что можно избавиться от свободного члена, добавив к векторам x
еще одну координату, тождественно равную 1:
x = (x1, x2, ..., x_n) --> x' = (x1, x2, ..., x_n, 1)
y = <w, x> + b = <w', x'>, где w' = (w1, w2, ..., w_n, b)
Можно считать, что зависимость выражается скалярным произведением
y = <w', x'>
В частности, можно рассмотреть задачу интерполяции многочленом степени d:
x1, x2, ..., x_m <- R
y1, y2, ..., y_m <- R
Требуется построить полином степени d, хорошо приближающий эту зависимость:
p(x) = w0*x^d + w1*x^d-1 + ... + w_d-1*x + w_d*1
w = (w0, w1, ..., w_d)
y = <w, z>, где z = (x^d, x^d-1, ..., x^2, x, 1) -- набор степеней x
по убыванию.
В методе наименьших квадратов мы минимизируем среднюю квадратичную ошибку:
Выборка:
x1, x2, x3, ..., x_m -- для этих значений x измерено значение y:
y1, y2, y3, ..., y_m
Ищем зависимость
y = <w, x>
Вектор w находим путем минимизации средней ошибки на тренировочной выборке:
1/m Sum (y_i - <w, x_i>)^2 --> min_w
Для метода наименьших существует точное решение задачи с помощью
псевдообратной матрицы Мура--Пенроуза:
[ x11 x12 ... x1n ] [ w1 ] [ y1 ]
A = [ x21 x22 ... x2n ] W = [ w2 ] y = [ y2 ] m >= n
[ . . . ] [ ...] [ y3 ]
[ . . . ] [ w_n] [ ...]
[ x_m1 x_m2 ... x_mn] [ y_m]
Мы ищем приблизительное "решение" системы
A W = Y
(система, вообще говоря, несовместна). Вектор W находится с помощью
псевдообратной матрицы:
W = pinv(A) Y
где, если столбцы матрицы A линейно независимы, то
pinv(A) = (A* A)^(-1) A*
Всё это очень красиво, за исключением единственной проблемы:
метод наименьших кадратов очень чувствителен к выбросам.
Выбросы -- это измерения, сделанные с большой ошибкой или
не соответствующие общей закономерности.
Пример:
берем статистику по использованию такси.
у -- время обслуживания одного клиента в зависимости от времени суток.
В среднем это времена порядки от 0.5 часа до 2 часов.
Но с небольшой вероятностью найдется клиент, который наймет машину на целый день
(10 часов). В методе наименьших квадратов подобные "нестандартные" клиенты
оказывают чрезмерно сильное влияние на результат.
Причина: в сумме квадратов отклонений большие отклонения приводят к
очень сильному увеличению средней ошибки => метод наименьших квадратов
имеет тенденцию минимизировать скорее отклонения для выбросов, чем для
регулярных точек.
Как быть???
1. Давайте попробуем минимизировать сумму модулей отклонений
вместо суммы квадратов отклонений.
Функция сумма квадратов отклонений называется MSE (Mean Squared Error).
MSE = 1/m Sum (y_i - <w, x_i>)^2 --> min_w
Функция сумма абсолютных величин отклонений называется MAE (Mean Absolute Err)
MAE = 1/m Sum |y_i - <w, x_i>| --> min_w
Проблема: функция модуля 1) не дифференцируема в нуле;
2) хоть и выпукла вниз, но не строго выпукла.
Большинство методов поиска минимума требует дифференцируемости
минимизируемой функции.
Есть два решения этой проблемы.
Первое решение -- использовать функцию Хубера вместо функций x^2 и |x|
Пусть мы имеем ошибку в одной точке
e = y_i - <w, x_i>
Квадрат ошибки вычисляется с помощью функции e^2.
Абсолютня величина ошибки вычисляется с помощью функции |e|.
Функция Хубера -- это компромис между функцией x^2 и |x|
{ x^2/2, при |x| <= delta
H_delta(x) = {
{ delta*(|x| - delta/2) при |x| > delta
Итак, метод Хубера использует минимизацию функции ошибки
HuberError(w) = 1/m Sum H_delta(y_i - <w, x_i>) --> min_w
Чаще всего delta = 1:
{ x^2/2, при |x| <= 1
H(x) = {
{ |x| - 1/2 при |x| > 1
В Питоне вычисление линейной регрессии с помощью функции Хубера
реализовано в модуле sklearn.linear_model:
from sklearn.linear_model import HuberRegressor
Пусть мы имеем выборку x (с добавленным признаком 1, решаем задачу без
свободного члена)
x = [ x1, x2, ..., x_m ], x_i <- R^n
y = [ y1, y2, ..., y_m ] y_i <- R
Создаем объект типа HuberRegressor для решения нашей задачи:
huber_reg = HuberRegressor(fit_intercept=False)
huber_reg.fit(x, y)
w = huber_reg.coef_
ПЕРЕРЫВ 15 минут до 13:55 (Moscow)
15:55 (Dishanbe)
==========================================
Минимизация функции средней абсолютной ошибки MAE
методами линейного программирования
Линейное программирование
У нас есть
1) набор переменных;
2) набор ограничений в виде равенств или неравенств;
3) целевая функция, которую надо минимизировать или максимизировать.
При этом все выражения линейно зависят от переменных.
Найти точку максимума или минимума.
Это называется задачей линейного программирования.
В Питоне для решения задач линейного программировани
используетс модуль pulp
from pulp import *
Рассмотрим простую задачу линейного программирования.
Есть
1) две переменные x, y (на плоскости)
2) ограничения в виде неравеств:
y >= 2 - 2*x
y <= 3 + 0.5*x
y <= 3*x
y <= 2 - 2*(x - 3)
y >= 2*(x - 2)
3) целевая функция:
x + y --> max
======================================
3 May 2022
Решение задачи линейной регрессии методами линейного
программирования
Задачи линейной регрессии
Надо найти линейную зависимость между независимой переменной
x <- R^n x = (x1, x2, ..., x_n)
и зависимой переменной
y = y(x), y <- R
При этом мы априори предполагаем, что зависимость линейная:
y = <w, x> + b, w <- R^n, b <- R -- свободный член
Часто избавляются от свободного члена b, добавляя к вектору x
еще одну координату, тождественно равную 1:
x --> x' = (x1, x2, ..., x_n, 1) <- R^n+1
Тогда линейная зависимость выражается просто скалярным произведением:
y = <w', x'>
где w' = (w1, w2, ..., w_n, b)
Такие вещи очень любят бизнесмены, экономисты, банковские работники:
вектор x -- признаки объекта (features)
вектор w -- это веса, которые мы придаем признакам.
Частный случай задачи:
пусть зависимость выражается многочленом степени d:
y = w0*x^d + w1*x^d-1 + ... + w_d-1*x + w_d (запись многочлена
по убыванию степеней).
Дана выборка
x1, x2, x3, ..., x_m, x_i -- это просто числа x_i <- R
y1, y2, y3, ..., y_m y_i <- R -- измеренные значения многочлена
Значение многочлена
f(x) = <w, z>, где z -- вектор степеней переменной x
z = (x^d, x^d-1, ..., x, 1) <- R^d+1
В методе наименьших квадратов мы минимизируем квадратичную ошибку:
MSE = 1/m Sum (f(x_i) - y_i)^2 --> min_w
Как было показано, этот метод очень чувствителен к выбросам.
Один из способов решения проблемы выбросов -- это использование функции
Хубера вместо квадрата разности. Функция Хубера для малых x ведет себя
как квадрат, а для больших x -- как модуль:
H_delta(e) = { 1/2 e^2 при |e| <= delta
{ delta(|e| - delta/2) при |e| > delta
Можно взять delta = 1.
В методе Хубера минимизируется штрафная функция
Huber = 1/m Sum H(f(x_i) - y_i) --> min_w
Довольно часто функцию потерь Хубера обозначают как MAE
(Mean Absolute Error), хотя это не то же самое в точности.
Давайте попробуем минимизировать функцию MAE:
MSE = 1/m Sum |f(x_i) - y_i| --> min_w
Никаких квадратов здесь нет, но мы не можем просто так применить
методы линейного программирования, потому что в симплекс-методе
минимизируемая (целевая) функция должна быть линейной, а здесь
у нас применяется модуль.
Как быть???
Есть красивый трюк, который позволяет избавиться от модуля.
Вводим для каждого объекта x_i две неотрицательные переменные
u_i >= 0, v_i >= 0, u_i, v_i <- R
Мы хотим избавиться от |f(x_i) - y_i|. Обозначим ошибку на i-м
объекте через e_i = f(x_i) - y_i
Представим e_i как разность двух неотрицательных переменных:
e_i = u_i - v_i
Тогда выполняется неравенство
|e_i| <= u_i + v_i
Поэтому нам достаточно минимизировать сумму переменных u_i + v_i
u_i + v_i --> min
В точке минимума либо v_i = 0 (когда e_i >= 0), либо u_i = 0
(когда e_i <= 0).
Например, e_i = 7
Если e_i = 10 - 3, где u_i = 10, v_i = 3, то эту пару можно заменить на
u_i = 7, v_i = 0,
разность u_i - v_i останется такой же, а сумма u_i + v_i уменьшится.
Значит, в точке минимума справедливо точное равенство
|e_i| = u_i + v_i
Теперь мы можем точно сформулировать нашу задачу как задачу линейного
программирования.
Есть тренировочная выборка
x1, x2, x3, ..., x_m <- R^n
y1, y2, y3, ..., y_m <- R
Найти w <- R^n такой, что выполняются органичения:
<w, x_i> - y_i = u_i - v_i для всех i = 1, 2, ..., m
u_i >= 0,
v_i >= 0,
при этом надо минимизировать целевую функцию
Sum (u_i + v_i) --> min_w
Попробуем ее решить на Питоне.
Домашнее задание: реализовать оконное приложение, где мы ищем
наилучший многочлен
f(x) = w0*x^d + w1*x^d-1 + ... + w_d-1*x + w_d
степени d по заданному набору узлов
(x_i, y_i)
Надо минимизировать функцию ошибки
MAE = 1/m Sum |f(x_i) - y_i|
ПЕРЕРЫВ 10 минут до 13:35 Moscow
15:35 Dushanbe
4 May 2022
Алгоритмы машинного обучения:
метод главных компонент (PCA --
Principal Component Analysis)
Этот метод относится к предварительной обработке объектов
еще до применения методов машинного обучения
(определение параметров линейной регрессии, построение линейного
классификатора и т.п.).
Какие могут быть проблемы с данными?
1. Всегда хорошо бы данные нормировать -- например, чтобы среднее
значение векторов x_i было бы равно нулю (центрировать данные).
2. Данные могут иметь очень большую размерность не по существу --
например, точки в n-мерном пространстве, представляющие наши объекты,
могут быть расположены на какой-то гиперплоскости (или вблизи
этой гиперплоскости), т.е. реальная размерность наших данных может
быть понижена.
Пример: рассматриваем распознавание рукописных цифр. Каждая цифра --
это матрица чисел размера 8x8, т.е. размерность нашего пространства
равна 64. Однако ясно, что по существу размерность должна быть меньше
хотя бы потому, что пиксели вблизи границ изображения (например, в углах
картинки) почти всегда черные (если цифры изображаются
белым по черному) и не несут никакой информации.
Метод главных компонент
1) центрирует данные, т.е. находит среднюю точку pca.mean_
2) строит ортонормированный базис pca.components_ (это массив из n элементов),
каждый элемент базиса -- это вектор длины 1, все векторы базиса
ортогональны друг другу и
-- дисперсии проекций векторов x_i на векторы базиса e1, e2, ..., en
убывают (нестрого).
Дисперсия -- это мера разброса точек
variance_i = Sum_j (<x_j, e_i> - E_k(<x_k, e_i>)^2
variance[1] >= variance[2] >= ... >= variance[n]
-- ковариация по разным направлениям равна нулю (мера зависимость
двух случайных величин), т.е. точки разбросаны независимо
по каждому из направлений e_i
Идея этого метода в том, что, если дисперсия по направлению k
близка к нулю, то то же самое будет и по всем направлениям l >= k,
т.е. реально все координаты, начиная с k-й координаты, можно
отбросить. Таки образом, мы сократим размерность нашей задачи.
Итак, имеем описания наших объектов в исходной системе координат.
В методе главных компонент выбирается новая ортогональная система координат,
в которой координаты, начиная с номера k, будут очень малы, и их можно
вообще не учитывать. Оставшиеся первый k-1 координата называются
главными компонентами. Каждая из главных компонент -- это некоторая
линейная комбинация исходных признаков. Причем все главные компоненты уже
не зависят друг от друга (ортогональны).
Пример:
в военном деле давно известно, что при стрельбе снаряды рассеиваются
по нормальному закону (закону Гаусса), по так называемым эллипсам
рассеивания. Полуоси этого эллипса соответствуют главным компонентам.
Если мы рассматриваем случайные точки в n-мерном пространстве,
распределенные, скорее всего, по закону, близкому к нормальному распределению,
то они тоже находятся внутри n-мерного эллипсоида; полуоси этого эллипсоида
и будут главными компонентами.
Если только первый k полуосей этого эллипсоида имеют значительный размер
(остальные очень малы), то мы можем взять проекции исходных точек
на гиперплоскость, натянутую на первые k полуосей эллипсоида.
Таким образом, мы потеряем очень мало информации, поскольку разброс
вдоль остальных направлений минимален.
Главные компоненты вычисляются путем приведения матрицы ковариаций
к диагональному виду.
В Питоне все это реализовано в модуле
sklearn.decomposition
(предварительная обработка данных, обучение еще без учителя).
Из этого модуля мы импортируем объект PCA (Principal Component Analysis).
Для того, чтобы обработать исходные данные x (массив векторов,
описывающих наши объекты x[i] -- n-мерный вектор, содержащий
признаки (features) i-го объекта), мы вызываем метод fit:
from sklearn.decomposition import PCA
. . .
pca = PCA() # Не указываем количество компонент, их число будет равно
# числу признаков
pca = PCA(n_components = 3) # Выделяются всего 3 первых главных компоненты
pca.fit(x) # переходим к новой системе координат, связанной
# с главными компонентами
x_tranformed = pca.transform(x)
Размерность векторов из x_transformed равна числу главных компонент,
указанных в конструкторе объекта pca.
Давайте напишем простую программу в Jupyter-Lab, которая иллюстрирует
метод главных компонент в 2-мерном случае.
Попробуем сгенерировать случайный набор данных (точки на плоскости),
применить для него метод главных компонент и нарисовать результат.
ПЕРЕРЫВ 10 минут до 11:18 (Moscow)
13:18 (Dushanbe)
Попробуем применить метод главных компонент к
базе данных рукописных цифр
from sklearn.datasets import load_digits
digits = load_digits()
Мы загрузили изображения рукописных цифр в виде матриц 8x8.
Сохранять только первые k компонент нужно, чтобы
1) уменьшить размерность задачи (число признаков) --
оставляем только существенные признаки
(причем независимые!), отбрасываем второстепенные признаки;
2) убрать шум из исходного множества объектов.
Дома:
загрузить базу данных рукописных цифр,
оставить только две главные компоненты,
нарисовать каждую цифру в виде цветной точки на плоскости.
===============================================
7 May 2022
Классификация рукописных цифр с помощью метода опорных векторов
Метод опорных векторов позволяет отделить друг от друга 2 класса точек
(binary classification). Строится разделяющая гиперплоскость в
(расширенном) пространстве и линейный классификатор
f(x) = <w, x> - b линейный классификатор
F(x) = K(w, x) - b ядровый классификатор
Классификатор принимает значения > 0 на точках первого класса и
значения < 0 на точках второго класса.
Как быть, когда у нас много классов?
Пусть мы имеем l классов. Тогда можно применить одну из двух стратегий:
1) one-versus-rest ("ovr") "каждый против остальных".
Имеем класс c, разбиваем точки всех классов на два подмножества
-- точки класса c;
-- точки всех остальных классов.
Строим линейный классификатор f_c(x)
Всего мы построим l классификаторов f1(x), f2(x), ..., f_l(x)
Возмем произвольный вектор x, надо определить, какому классу он
принадлежит.
predict(x) = maxarg f_i(x), i=1,2,...,l
2) one-versus-one ("ovo") "каждый против каждого".
Для каждой пары классов i < j строится классификатор f_i,j(x)
разделяющий точки только двух классов i и j. Всего будет построено
(l-1) + (l-2) + ...+ 1 = (l-1)*l/2
Для рукописных цифр классов 10, значит будет построено
45 классификаторов.
Рассмотрим конкретный вектор. Применяем все 45 классификаторов.
Применение классификатора f_i,j(x) -- это как бы матч двух команд.
В зависимости от результата сравнения
f_i,j(x) >= 0 -- выиграла команда i, i-й команде прибавляется
очко к ее результату,
f_i,j(x) < 0 -- выиграла команда j, j-й команде прибавляется
очко.
Какая-то из команд s выигрывает круговой турнир, т.е. набирает
максимальное количество очков -- это и будет значение
классификатора
predict(x) = s
В методе опорных векторов в стандартной библиотеке Питона sklearn.svm
для многоклассовой классификации применяется метод "one-versus-one".
Применим метод опорных векторов для классификации рукописных цифр
из базы данных NIST (National Institute of Standards and Technology),
созданной в 1998 г.
import numy as np
from sklearn.datasets import load_digits
digits = load_digits()
Что мы имеем в результате?
digits.DESCR -- описание базы данных:
======================================================
Optical recognition of handwritten digits dataset
--------------------------------------------------
**Data Set Characteristics:**
:Number of Instances: 1797
:Number of Attributes: 64
:Attribute Information: 8x8 image of integer pixels in the range 0..16.
:Missing Attribute Values: None
:Creator: E. Alpaydin (alpaydin '@' boun.edu.tr)
:Date: July; 1998
This is a copy of the test set of the UCI ML hand-written digits datasets
https://archive.ics.uci.edu/ml/datasets/Optical+Recognition+of+Handwritten+Digits
The data set contains images of hand-written digits: 10 classes where
each class refers to a digit.
Preprocessing programs made available by NIST were used to extract
normalized bitmaps of handwritten digits from a preprinted form. From a
total of 43 people, 30 contributed to the training set and different 13
to the test set. 32x32 bitmaps are divided into nonoverlapping blocks of
4x4 and the number of on pixels are counted in each block. This generates
an input matrix of 8x8 where each element is an integer in the range
0..16. This reduces dimensionality and gives invariance to small
distortions.
For info on NIST preprocessing routines, see M. D. Garris, J. L. Blue, G.
T. Candela, D. L. Dimmick, J. Geist, P. J. Grother, S. A. Janet, and C.
L. Wilson, NIST Form-Based Handprint Recognition System, NISTIR 5469,
1994.
.. topic:: References
- C. Kaynak (1995) Methods of Combining Multiple Classifiers and Their
Applications to Handwritten Digit Recognition, MSc Thesis, Institute of
Graduate Studies in Science and Engineering, Bogazici University.
- E. Alpaydin, C. Kaynak (1998) Cascading Classifiers, Kybernetika.
- Ken Tang and Ponnuthurai N. Suganthan and Xi Yao and A. Kai Qin.
Linear dimensionalityreduction using relevance weighted LDA. School of
Electrical and Electronic Engineering Nanyang Technological University.
2005.
- Claudio Gentile. A New Approximate Maximal Margin Classification
Algorithm. NIPS. 2000.
======================================================================
>>> dir(digits)
['DESCR', 'data', 'feature_names', 'frame', 'images', 'target', 'target_names']
>>> len(digits.data)
1797
digits.data -- матрица, каждая ее строка -- массив из 64 чисел 0..16.
Этот массив можно рассматривать как матрицу 8x8,
представляющую изображение цифры.
digits.images -- массив матриц 8x8, представляющих изображения цифр.
digits.target -- номер класса, т.е. номер цифры, изображение которой
в массивах data и images
Всего изображений 1797.
Наша цель:
применить классификаторы, построенные по методу опорных векторов
с различными параметрами, и определить, при каком наборе параметров получается
наилучший результат.
Для этого база данных разделяется на две части:
-- тренировочная часть;
-- тестовая часть.
Классификаторы создаются по тренировочной части и
затем проверяются на тестовой части.
Результат -- это процент правильных предсказаний на тестовой части.
Пусть размер тренировочной части составляет 4/5 = 80% от всей базы,
размер тестовой части составляет 1/5 = 20% от базы.
Какие параметры метода SVM мы используем?
kernel = [ "linear", "poly", "rbf" ]
Для ядра poly используется степень полинома
degree = [ 2, 3, 4, 5 ]
Гиперпараметр C
C = [ 0.5, 1., 2., 4., 10. ]
Всего число комбинаций 3*5 + 3 = 18 комбинаций.
Для какой комбинации будет наилучший результат?
ПЕРЕРЫВ 10 минут до 11:15 (Moscow)
13:15 (Dushanbe)
import numpy as np
from sklearn.datasets import load_digits
from sklearn.svm import SVC
import matplotlib.pyplot as plt
%matplotlib inline
digits = load_digits()
m = len(digits.data)
# Divide the data into train and test sets
m_train = int(m*0.8)
m_test = m - m_train
# Shuffle data before separation
x = digits.data.copy()
y = digits.target.copy()
p = np.random.permutation(m)
x = x[p]
y = y[p]
# plt.imshow(x[234].reshape(8, 8), cmap="gray_r")
# print(y[234])
x_train = x[:m_train]
x_test = x[m_train:]
y_train = y[:m_train]
y_test = y[m_train:]
len(x_train), len(x_test)
# C = 2.
C = 5.
# kernel = "linear"
kernel = "poly"
# kernel = "rbf"
degree = 3
classifier = SVC(C=C, kernel=kernel, degree=degree)
classifier.fit(x_train, y_train)
y_predicted = classifier.predict(x_test)
# Define the precision of predictions
correct_predictions = 0
for i in range(len(y_test)):
if y_predicted[i] == y_test[i]:
correct_predictions += 1
score = correct_predictions / len(x_test)
print("score =", score)
# Find an error prediction and show a problematic image
idx = (-1)
for i in range(len(y_test)):
if y_predicted[i] != y_test[i]:
idx = i
break
print("predicted:", y_predicted[idx], " real:", y_test[idx])
plt.imshow(x_test[idx].reshape(8, 8), cmap="gray_r")
Домашнее задание
Для базы данных digits определить, какой набор параметров дает
наиболее точные предсказания при использовании метода опорных
векторов.
Гиперпараметр C:
C = [ 0.5, 1., 2., 4., 10. ]
Параметр "Ядро":
kernel = [ "linear", "poly", "rbf" ]
Для ядра poly используется параметр "степень полинома":
degree = [ 2, 3, 4, 5 ]
Всего получаем 18 наборов параметров:
5*3 + 3 = 18
(поскольку параметр degree используется только для полиномиального ядра).
Вся база цифр разбивается на две части:
80% -- тренировочная часть
20% -- тестовая часть
Классификатор тренируется на тренировочной части и проверяется
на тестовой части. Про проверке на тестовой части мы определяем
процент правильных предсказаний.