Самостоятельная работа состоит в решении задач по темам, которые обсуждаются в лекциях. По каждой теме нужно решить одну задачу из предлагаемого списка. Номер задачи равен номеру студента в журнале по модулю n, где n — число задач по данной теме. Решения задач присылайте мне на электронную почту в виде одного или нескольких файлов, присоединенных к письму (если файлов много, то лучше присылать zip-архив с файлами в поддиректории). Адрес моей электронной почты:
vladibor266 (собака) gmail.com
Тема письма должна начинаться со слов "Душанбе 2022".

Задачи для зачета

Содержание

1. Занятие 16.04.2022
   Работа с матрицами в Python'е с использованием модуля numpy. Метод Гаусса, обращение матрицы, решение системы линейных уравнений. Особенности реализации при использовании модуля numpy.
2. Занятие 19.04.2022
   Метод наименьших квадратов. Задача линейной регрессии как пример задачи машинного обучения. Решение задачи регрессии методом наименьших квадратов.
3. Занятие 20.04.2022
   Метод сопряженных градиентов для минимизации квадратичной функции. Построенные на его основе методы Хестенса-Штифеля (минимизация квадратичной функции) и Флетчера-Ривса (минимизация произвольной выпуклой функции).

   Создание оконных и графических приложений на Python'е с использованием модуля tkinter. Оконная программа, иллюстрирующая метод наименьших квадратов.

4. Занятие 23.04.2022
   Метод опорных векторов, классический линейный вариант.
5. Занятие 26.04.2022
   Нелинейный метод опорных векторов: добавление некоторого набора базовых функций от исходных признаков объекта и расширение описания объекта с их помощью. Построение линейного классификатора в пространстве более высокой размерности.
6. Занятие 27.04.2022
   Реализация нелинейного метода опорных векторов путем добавления всех мономов степени не выше d к исходному набору признаков. Построение линий уровня функции от двух переменных
   f(x, y) = level = const
7. Занятие 28.04.2022
   Ядровые методы и их использование в нелинейном методе опорных векторов.
8. Занятие 30.04.2022
   Задача линейной регрессии: проблема выбросов. Функции оценки ошибки MSE (Mean Squared Error) и MAE (Mean Absolute Error). Решение задачи регрессии с использованием функции ошибки Хубера. Линейное программирование.
9. Занятие 3.05.2022
   Решение задачи линейной регрессии методами линейного программирования.
10. Занятие 4.05.2022
    Метод главных компонент.
11. Занятие 7.05.2022
    Классификация рукописных цифр с помощью метода опорных векторов

Домашние задания

1. Работа с матрицами в модуле numpy
2. Линейная регрессия и метод наименьших квадратов
3. Метод опорных векторов
4. Решение задачи линейной регрессии путем минимизации фунции средней абсолютной величины ошибки (MAE — Mean Absolute Error) методами линейного программирования
5. Классификация рукописных цифр с помощью метода опорных векторов

Журнал группы ПМИ-4, весенний семестр 2022

Тексты лекций (записи на "виртуальной доске")

import numpy as np

>>> a = np.array([2, 3, 5, 7])

>>> a = np.array([[1., 2., 3., 2.], [4., 5., 6., 3.], [7., 8., 0., 5.]])
>>> a
array([[1., 2., 3., 2.],
       [4., 5., 6., 3.],
       [7., 8., 0., 5.]])
>>> a.shape
(3, 4)

>>> a = np.array([1., 2., 3., 4.])
>>> b = np.array([10., 20., -10., -20.])
>>> a+b
array([ 11.,  22.,  -7., -16.])
>>> a*b
array([ 10.,  40., -30., -80.])
>>> a*7.
array([ 7., 14., 21., 28.])
>>> np.sin(a)
array([ 0.84147098,  0.90929743,  0.14112001, -0.7568025 ])

>>> a.dot(b)
-60.0
>>> np.dot(a, b)
-60.0

>>> a @ b
-60.0

>>> a = np.array([[1, 2], [3, 4]])
>>> a
array([[1, 2],
       [3, 4]])
>>> b = np.array([[2, 1], [-1, 2]])
>>> b
array([[ 2,  1],
       [-1,  2]])
>>> a.dot(b)
array([[ 0,  5],
       [ 2, 11]])
>>> a @ b
array([[ 0,  5],
       [ 2, 11]])

В файле "npmatrix.py" для матриц, представленных numpy-массивами, реализован метод Гаусса приведения матрицы к ступенчатому виду (функция echelonize), а также обратный проход в методе Гаусса (функция reversePass). Добавить в тот же файл реализацию одной из следующих двух функций.

1. Вычисление обратной матрицы:

   def invMatrix(a, eps = 1e-8):
       """Return an inverse matrix of a matrix a"""
   

2. Решение системы линейных уравнений с невырожденной квадратной матрицей:

   def solveLinearSystem(a, b, eps = 1e-8):
       """Solve a linear system
           a x = b
       where a is non-singular square matrix
       and b is array of free terms. Return an array x"""
   

Описание метода наименьших квадратов

Создание оконных приложений в Python'е с помощью модуля tkinter

Модуль tkinter предоставляет, наверно, самый простой из всех возможных интерфейсов для создания интерактивных оконных приложений на языке Python. Рассмотрим пример простейшей оконной программы, которую можно использовать как заготовку для создания собственных приложений. Программа создает окно-форму верхнего уровня, в верхней части которого расположены управляющие элементы, в нашем случае две кнопки "Draw" и "Clear". Также в верхней части окна расположен слайдер, с помощью которого можно изменять значение целочисленной переменной (мы будем ее использовать для задания степени многочлена, аппроксимирующего зависимость между переменными x и y в методе наименьших квадратов). В нижней части окна верхнего уровня расположено дочернее окно, предназначенное для рисования. В нем изображена координатная сетка. По нажатию кнопки "Draw" в нем рисуется график функции, заданной непосредственно в тексте программы. Также можно в этом окне кликать мышкой, программа в ответ рисует крестики в этих точках; цвет крестика соответствует номеру нажатой клавиши мыши. Нажатие кнопки "Clear" в верхней части стирает график функции и крестики в отмеченных точках.

Приложение использует два файла: файл "plotFunc.py" содержит непосредственно саму оконную программу; файл "R2Graph.py" содержит определения классов Вектор и Точка на плоскости R2Vector и R2Point, которые используются в программе. Файл R2Graph.py импортируется в оконной программе как вспомогательный модуль. Архив всех файлов: "plotFunc.zip".

Запускается программа командой

    python3 plotFunc.py

Реализовать на языке Python3 с использованием модуля tkinter интерактивное приложение, иллюстрирующее метод наименьших квадратов.

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

Содержание: постановка задачи нахождения линейного классификатора для разделения двух классов объектов, которые описываются векторами n-мерного проостранства. Метод опорных векторов, который вычисляет гиперплоскость коразмерности 1, разделяющую два класса точек; различные функции потери. Графическая программа на Python+tkinter, иллюстрирующая метод опорных векторов: "svm.py" (программа использует также модуль R2Graph.py, архив всех файлов: "svm.zip").

Описание метода опорных векторов (Support Vector Machine)

В классическом методе опорных векторов мы имеем два класса объектов, каждый объект описывается n признаками, представляемыми вещественными числами. Таким образом, объекты x₁, x₂, … x_k представляются точками n-мерного пространства; координаты такой точки равны значениям признаков объекта. В тренировочной выборке объекты x_i первого класса отмечены значением y_i = 1, объекты x_j второго класса — значением y_j = -1. В методе опорных векторов SVM (Support Vector Machine) строится линейный классификатор f(x), который определяется как гиперплоскость, разделяющая два класса точек:

В классическом методе опорных векторов далеко не всегда можно разделить гиперплоскостью два класса точек в n-мерном пространстве:

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

Итак, пусть отображение

В простейшем случае в качестве функций φ₁(x₁, … x_n), φ₂(x₁, … x_n), …, φ_m(x₁, … x_n) можно взять всевозможные мономы степени не выше d от переменных x₁, x₂, … x_n (то есть от исходных признаков).

Пусть, например, d = 2. Пусть число исходных признаков n = 2, т.е. объекты представляются точками плоскости R². Тогда расширенный набор признаков представляется всеми мономами степени не выше 2 от переменных x₁, x₂:

Программа "svmNonlinear.py" иллюстрирует использование нелинейного метода опорных векторов. Исходные объекты представляются точками плоскости двух цветов, отмечаемые кликами левой и правой клавиш мыши. В качестве расширенных признаков используются все мономы степени не выше d от координат x и y очередной точки; степень d задается с помощью слайдера в пределах от 1 до 5. Программа после построения нелинейного классификатора f(x) изображает линию нулевого уровня функции

Классический случай: степень 1

Квадратичные мономы: степень 2

Кубический многочлен: степень 3

Многочлен степени 4

Многочлен степени 5

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

# Non-linear Support Vector Machine method

from scipy.optimize import minimize
import numpy as np
import itertools
import random
from math import sin, cos
import matplotlib.pyplot as plt
from skimage import measure

get_ipython().run_line_magic('matplotlib', 'inline')
plt.axes(aspect="equal")

def monomials(x, degree=1):
    m = [1.]
    for d in range(1, degree+1):
        for c in itertools.combinations_with_replacement(x, d):
            m.append(np.prod(c))
    return np.array(m)


m1 = 100
class1 = [ np.array([
    random.normalvariate(1., 1.5), random.normalvariate(1., 0.5)
]) for i in range(m1) ]

m2 = m1
angle1 = -np.pi*2/3
angle2 = np.pi*2/3
ex = np.array([1., 0.]); ey = np.array([0., 1.])
class2 = []
for i in range(m2):
    phi = random.normalvariate(0., 1.)
    r = random.normalvariate(5., 0.5)
    p = ex*cos(phi)*r + ey*sin(phi)*r
    class2.append(p.copy())
    
# plt.scatter([c[0] for c in class1], [c[1] for c in class1])
# plt.scatter([c[0] for c in class2], [c[1] for c in class2])

x = class1 + class2
x = np.array(x)
y = [1.]*len(class1) + [-1.]*len(class2)
y = np.array(y)
m = len(x)
perm = np.random.permutation(m)
x = x[perm]
y = y[perm]

plt.scatter(
    [c[0] for c in x], [c[1] for c in x],
    color=[ "r" if yy > 0. else "g" for yy in y ]
)

degree = 3
xext = np.array([monomials(xx, degree) for xx in x])
n = len(xext[0])
print("n =", n)
C = 10.    # Hyperparameter
w = np.array([0.]*n)
b = 0.

def hingeLoss(c):
    return max(1. - c, 0.)

def errorFunction(wb):
    w = wb[:-1]
    b = wb[-1]
    w2 = w @ w
    s = 0.
    m = len(xext)
    for i in range(m):
        c = y[i]*(w @ xext[i] - b)
        s += hingeLoss(c)
    return w2 + (C/m)*s

def classifier(p):
    pext = monomials(p, degree)
    return w @ pext - b

res = minimize(errorFunction, np.zeros(n + 1))
w = res.x[:-1]
b = res.x[-1]
print("w =", w)
print(" b =", b)

# Draw the separating line
xmin = min([xx[0] for xx in x]) - 1.
xmax = max([xx[0] for xx in x]) + 1.
ymin = min([xx[1] for xx in x]) - 1.
ymax = max([xx[1] for xx in x]) + 1.
rows = 100
cols = rows

def xcoord(j):
    dx = (xmax - xmin)/cols
    return xmin + j*dx
def ycoord(i):
    dy = (ymax - ymin)/rows
    return ymin + i*dy

# a is a matrix of dimension (rows, cols)
a = np.array([[0.]*steps for iy in range(rows)])
for i in range(rows):
    for j in range(cols):
        xx = np.array([xcoord(j), ycoord(i)])
        a[i, j] = classifier(xx)
nullContours = measure.find_contours(a, 0.)
for c in nullContours:
    plt.plot(
        [xcoord(cc[1]) for cc in c], 
        [ycoord(cc[0]) for cc in c],
        color="blue"
    )

Вот результат работы этой программы (запущенной в среде Jupyter-Lab):

Исходный код программы: "svmNonlinear.ipynb".

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

# Drawing Isolines
import numpy as np
import matplotlib.pyplot as plt
from skimage import measure

get_ipython().run_line_magic('matplotlib', 'inline')
plt.axes(aspect="equal")

def f(x, y):
    return ((x - 2) + (y - 1))**2/4 + ((x - 2) - (y - 1))**2

xmin = -5.; xmax = 10.
ymin = -5.; ymax = 10.
m = 100; n = 100

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

a = np.array([[0.]*n for i in range(m)])
for i in range(m):
    y = ycoordinate(i)
    for j in range(n):
        x = xcoordinate(j)
        a[i, j] = f(x, y)
        
for l in range(20):
    level = 0.1*(2**l)
    contours = measure.find_contours(a, level)
    for c in contours:
        plt.plot(
            [ xcoordinate(cc[1]) for cc in c ],
            [ ycoordinate(cc[0]) for cc in c ]
            #, color="blue"
        )

Вот результат работы этой программы (запущенной в среде Jupyter-Lab):

Описание ядрового варианта метода опорных векторов

Вот несколько скриншотов оконной программы на Python+tkinter, иллюстрирующей применение ядрового метода опорных векторов:
Гауссовское ядро (RBF)

Полиномиальное ядро, степень 2

Полиномиальное ядро, степень 3

В этой программе используется модуль sklearn и его подмодуль svm, в котором реализован ядровый метод опорных векторов. В программе кликами мыши задается массив точек points, при этом точки, отмеченные кликом левой клавиши мыши, относятся к первому классу, другой клавишей — ко второму, номера клавиш мыши содержится в массиве mouseButtons. Вот как выглядит код, вызывающий метод опорных векторов.

Описания и глобальные переменные:

import numpy as np
from sklearn.svm import SVC

classifier = SVC(kernel="rbf")

        global classifier
        . . .
        if kernel == "rbf":
            classifier = SVC(C=c, kernel="rbf")
        else:
            classifier = SVC(C=c, kernel="poly", degree=polynomialDegree)

        data = [
            [points[i].x, points[i].y]
            for i in range(len(points))
        ]
        data = np.array(data)

        y = [
            (1. if mouseButtons[i] == 1 else (-1.))
            for i in range(len(points))
        ]
        y = np.array(y)

        classifier.fit(data, y)

def classifierValue(p):
    x = np.array([[p.x, p.y]])
    # return classifier.predict(x)
    return classifier.decision_function(x)

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

Пример применения ядрового метода опорных векторов в среде Jupyter-Lab

В среде Jupyter-Lab мы сначала генерируем на плоскости случайные точки двух классов так, чтобы их невозможно было разделить прямой линией. Затем мы применяем ядровый метод опорных векторов с Гауссовским ядром "rbf" (Radial Basis Function) и рисуем линию, которая разделяет два класса. Вот код программы (в среде Jupyter-Lab):

# Kernel Support Vector Machine Method

from scipy.optimize import minimize
import numpy as np
import random
from math import sin, cos
import matplotlib.pyplot as plt
from skimage import measure

from sklearn.svm import SVC

get_ipython().run_line_magic('matplotlib', 'inline')
plt.axes(aspect="equal")

m1 = 100
class1 = [ np.array([
    random.normalvariate(1., 1.5), random.normalvariate(1., 0.5)
]) for i in range(m1) ]

m2 = m1
angle1 = -np.pi*2/3
angle2 = np.pi*2/3
ex = np.array([1., 0.]); ey = np.array([0., 1.])
class2 = []
for i in range(m2):
    phi = random.normalvariate(0., 1.)
    r = random.normalvariate(5., 0.5)
    p = ex*cos(phi)*r + ey*sin(phi)*r
    class2.append(p.copy())
    
x = class1 + class2
x = np.array(x)
y = [1.]*len(class1) + [-1.]*len(class2)
y = np.array(y)
m = len(x)
perm = np.random.permutation(m)
x = x[perm]
y = y[perm]

plt.scatter(
    [c[0] for c in x], [c[1] for c in x],
    color=[ "r" if yy > 0. else "g" for yy in y ]
)

degree = 2
C = 2.    # Hyperparameter

# classifier = SVC(C=C, kernel="poly", degree=degree)
classifier = SVC(C=C, kernel="rbf")
classifier.fit(x, y)

def classifierValue(x):
    xx = np.array([  [x[0], x[1]]  ])
    # return classifier.predict(xx)
    return classifier.decision_function(xx)

# Draw the separating line
xmin = min([xx[0] for xx in x]) - 1.
xmax = max([xx[0] for xx in x]) + 1.
ymin = min([xx[1] for xx in x]) - 1.
ymax = max([xx[1] for xx in x]) + 1.
steps = 100

def xcoord(j):
    dx = (xmax - xmin)/steps
    return xmin + j*dx
def ycoord(i):
    dy = (ymax - ymin)/steps
    return ymin + i*dy

# a is a matrix of dimension (steps, steps)
a = np.array([[0.]*steps for iy in range(steps)])
for i in range(steps):
    for j in range(steps):
        xx = np.array([xcoord(j), ycoord(i)])
        a[i, j] = classifierValue(xx)
nullContours = measure.find_contours(a, 0.)
for c in nullContours:
    plt.plot(
        [xcoord(cc[1]) for cc in c], 
        [ycoord(cc[0]) for cc in c],
        color="blue"
    )

Вот результат работы этой программы (запущенной в среде Jupyter-Lab):

Исходный код программы: "svmKernel.ipynb".

В задании 2 варианта. Студенты с нечетными номерами по журналу делают вариант 1, с четными — вариант 2. В обоих вариантах надо изменить оконную программу "svm.py" (полный архив всех исходных файлов "svm.zip"), в которой реализован классический линейный метод опорных векторов:

1. Переделать программу "svm.py", в которой реализован линейный метод опорных векторов. Надо дополнить признаки объектов всеми мономами (произведениями) степени ≤ d от исходных признаков, число d должно задаваться слайдером. Например, при d=2 мы вместо точки
   x = (x₀, x₁) ∈ R²
   используем "расширенную" точку
   x' = (1, x₀, x₁, x₀², x₀x₁, x₁²) ∈ R⁶,
   при d=3 получится 10-мерное пространство и т.д. Линейный классификатор вычисляется для расширенного пространства. Программа должна нарисовать линию, разделяющую точки двух классов.
2. Переделать программу "svm.py", в которой реализован классический линейный метод опорных векторов. Модифицированная программа должна использовать ядровый метод опорных векторов, реализованный в стандартных модулях Питона в виде класса SVC:

   from sklearn.svm import SVC
   

   Программа должна применить классификатор SVC и нарисовать линию, разделяющую точки двух классов. В программе должна быть реализована возможность выбора ядра из двух вариантов: "rbf" (Гауссовское ядро), "poly" (полиномиальное ядро). В случае полиномиального ядра степень полинома должна задаваться с помощью слайдера.

Занятие 30.04.2022
Проблема выбросов в задаче линейной регрессии и способы ее решения

В задаче линейной регрессии мы строим зависимость между независимой переменной x∈Rⁿ и зависимой переменной y∈R. Зависимость описывается вектором w∈Rⁿ и свободным членом b∈R:

Если вместо средней величины квадрата разности для каждого измерения минимизировать абсолютную величину разности, то построенная модель будет в значительно меньшей степени зависеть от выбросов. Среднее значение абсолютной величины ошибки обозначается через MAE (Mean Average Error):

Хорошим компромиссом является использование функции потерь Хубера, зависящей от параметра δ:

На следующей картинке изображен график функции Хубера для δ=1 (синий цвет), а также графики функций y = |x| (красный) и y = x²/2 (зеленый):

Видно, что функция потерь Хубера объединяет в себе лучшие свойства квадратичной функции ошибки и функции абсолютной величины ошибки: она всюду дифференцируема, так же как квадратичная функция, но при этом не растет чрезмерно быстро, как и абсолютная величина. В задаче регрессии модель строится путем минимизации функции средней ошибки, вычисление которой использует функцию потери Хубера H_δ(x):

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

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

В использованной программе вычисляется полином заданной степени, построенный по набору узлов, которые отмечаются кликами мыши. При вычислении наилучшего полинома используется либо метод наименьших квадратов, при этом минимизируется функция ошибки MSE; либо модель строится путем минимизации функции средней ошибки с использованием функции Хубера (это близко к функции ошибки MAE). Вот исходный код программы на Python'е: "huberRegression.py". Программа использует также модуль "R2Graph.py", архив всех файлов: "huberRegression.zip".

Задача линейного программирования

Теперь мы попробуем минимизировать функцию ошибки MAE (Mean Absolute Error) методами линейного программирования (симплекс-метод и т.п.). Рассмотрим сначала общую постановку задачи линейного программирования и средства Питона для ее решения.

Формулировка задачи линейного программирования

В задаче линейного программирования заданы

В Питоне для решения задач линейного программировани используетс модуль pulp:

from pulp import *

Способ записи формулировки задачи в модуле pulp очень простой. Сначала мы создаем объект, отвечающий за решение нашей задачи:

    problem = LpProblem("Regression", LpMinimize)

Затем надо определить 3 вещи.

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

       x = LpVariable("x")         # Свободная вещественная переменная x
       y = LpVariable("y", 0)      # Неотрицательная вещ. переменная y >= 0
       v = LpVariable("v", lowBound=0) # Неотрицательная вещ. переменная v >= 0
       z = LpVariable("z", upBound=0)  # Вещ. переменная z <= 0
       t = LpVariable("t", 0, 3)   # Вещ. переменная 0 <= t <= 3
       i = LpVariable("i", 0, cat="Integer") # Целочисленная переменная i >= 0
   

   Можно описать сразу целый список (или словарь) переменных:

       n = 3
       x = LpVariable.dicts("x", range(n), lowBound=0)
   

   Описаны неотрицательные вещественные переменные x[0], x[1], x[2] с именами x_0, x_1, x_2.
2. Ограничения в форме равенств или неравенств. Ограничения добавляются к нашему объекту, описывающему задачу, при помощи операции +=. Примеры:

       problem += 2*x + 3 <= 0
       problem += x - 2*y == 1
   

   Возможны только три операции сравнения: <=, ==, >= (равенство или нестрогое неравенство).
   Для суммирования списка выражений можно использовать функцию lpSum:

       a = [0.5, -2.7, 1.9]
       x = LpVariable.dicts("x", range(3), lowBound=0)
       problem += lpSum([a[i]*x[i] for i in range(3)]) <= 10
   

   Здесь мы добавили ограничение

       0.5*x_0 - 2.7*x_1 + 1.9*x_2 <= 10.
   

3. Целевая функция, которую мы минимизируем или максимизируем. В отличие от ограничений, целевая функция только одна. Она также добавляются к объекту, описывающему задачу, при помощи операции +=. Запись целевой функции отличается от записи ограничений тем, что в ней не присутствует операция сравнения. Пример:

       problem += x + y - 2*z
   

   В записи целевой функции тоже можно использовать функцию суммирования lpSum:

       problem += 0.5*w + lpSum([x[i] for i in range(3)])
   

Пример решения простой задачи линейного программирования средствами Питона

Работаем в среде Jupyter-Lab:

from pulp import *
import matplotlib.pyplot as plt
import numpy as np

plt.axes(aspect="equal")
xmin = -1.
xmax = 5.
ymin = -1.
ymax = 5.
plt.xlim(xmin, xmax)
plt.ylim(ymin, ymax)
plt.grid()

xx = np.linspace(xmin, xmax, 50)
plt.plot(xx, 2. - 2.*xx)
plt.plot(xx, 3. + 0.5*xx)
plt.plot(xx, 3.*xx)
plt.plot(xx, 2. - 2.*(xx - 3.))
plt.plot(xx, 2.*(xx - 2.))
plt.plot(xx, 0.*xx)
plt.plot([0., 0.], [ymin, ymax])
plt.fill(
    [1, 2, 3, 2, 6/5, 2/5], 
    [0, 0, 2, 4, 18/5, 6/5], 
    color="lightGray", alpha=0.5
)  

x = LpVariable("x", 0)    # x >= 0
y = LpVariable("y", 0)    # y >= 0

myProblem = LpProblem("FirstProblem", LpMaximize)
# myProblem = LpProblem("FirstProblem", LpMinimize)

myProblem += y >= 2. - 2.*x
myProblem += y <= 3. + 0.5*x
myProblem += y <= 3.*x
myProblem += y <= 2. - 2.*(x - 3.)
myProblem += y >= 2.*(x - 2.)

myProblem += x + y

myProblem.solve()
print("Status: ", LpStatus[myProblem.status])

print("Solution: x =", value(x), " y =", value(y))

Status:  Optimal
Solution: x = 2.0  y = 4.0

В задаче линейной регрессии мы строим зависимость между независимой переменной x∈Rⁿ и зависимой переменной y∈R. Зависимость описывается вектором w∈Rⁿ и свободным членом b∈R:

Если вместо средней величины квадрата разности минимизировать абсолютную величину разности, то построенная модель будет в значительно меньшей степени зависеть от выбросов. Среднее значение абсолютной величины ошибки обозначается через MAE (Mean Average Error):

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

Трюк с избавлением от модуля

Приведенная формулировка задачи минимизации не позволяет сразу применить методы линейного программирования из-за использования функции модуля. Однако можно применить следующий трюк. Пусть мы минимизируем сумму модулей некоторых величин |e_i|, i=1,…,m. Каждую величину e_i, которая может быть положительной и отрицательной, можно представить в виде разности двух положительных величин:

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

Опустим штрихи в обозначениях x', w', чтобы не загромождать запись. Итак, нам надо минимизировать следующую сумму:

Пример решения задачи регрессии методом линейного программирования

from pulp import *
x = [-6, -4, -3, -1, 1, 2, 3]
y = [-2, -1, 1, 2, 2, 3, 4]
m = len(x)
problem = LpProblem("Regression", LpMinimize)

# Variables:
w = LpVariable.dicts("w", range(2))
u = LpVariable.dicts("u", range(m), lowBound=0)
v = LpVariable.dicts("v", range(m), lowBound=0)

# Constraints:
for i in range(m):
    problem += w[0]*x[i] + w[1]*1 - y[i] == u[i] - v[i]

# Objective function:
problem += lpSum([u[i] + v[i] for i in range(m)])

# Solve the problem using the default solver
problem.solve()

# Print the status of solution
print("Status:", LpStatus[problem.status])

# Print the result: w = (w_0, w_1)
print("Regression coefficients:", value(w[0]), value(w[1]))

Status: Optimal
Regression coefficients: 0.625 1.75

Реализовать на языке Python3 с использованием модуля tkinter интерактивное оконное приложение, где мы ищем наилучший многочлен

Пользователь отмечает кликами мыши множество точек на координатной плоскости, точки отмечаются крестиками. С помощью слайдера в оконной форме задается степень аппроксимирующего многочлена. По нажатию на кнопку "Draw" программа вычисляет коэффициенты аппроксимирующего многочлена, применяя методы линейного программироания модуля pulp, и рисует график многочлена в окне:

В качестве образца использовать программу "huberRegression.py". Программа использует также модуль "R2Graph.py", архив всех файлов: "huberRegression.zip".

Этот метод относится к предварительной обработке объектов еще до применения методов машинного обучения (определения параметров линейной регрессии, построения классификатора и т.п.).

Какие могут быть проблемы с данными?

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

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

Метод главных компонент
1) центрирует данные, т.е. находит среднюю точку pca.mean_
2) строит ортонормированный базис pca.components_ (это массив из n элементов), каждый элемент базиса — это вектор нормы 1, все векторы базиса ортогональны друг другу, и

• дисперсии проекций векторов x_i на векторы базиса e₁, e₂, ..., e_n убывают (нестрого).
  Дисперсия — это мера разброса точек
  var_i = ∑_j ( ⟨x_j, e_i⟩ − (1/n) ∑_k⟨x_k, e_i⟩ )²,
  var₁ ≥ var₂ ≥ ... ≥ var_n.
• ковариация по разным направлениям равна нулю (мера зависимости двух случайных величин), т.е. точки разбросаны независимо по каждому из направлений e_i.

Итак, имеем описания наших объектов в исходной системе координат. В методе главных компонент выбирается новая ортогональная система координат, в которой координаты, начиная с номера k+1, будут очень малы, и их можно вообще не учитывать. Оставшиеся первые k координат называются главными компонентами. Каждая из главных компонент — это некоторая линейная комбинация исходных признаков. Причем все главные компоненты уже не зависят друг от друга (ортогональны).

Пример: в военном деле давно известно, что при стрельбе снаряды рассеиваются по нормальному закону (закону Гаусса), по так называемым эллипсам рассеивания. Полуоси этого эллипса соответствуют главным компонентам.

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

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

Главные компоненты вычисляются путем приведения матрицы ковариаций к диагональному виду.

    sklearn.decomposition

from sklearn.decomposition import PCA
. . .
pca = PCA()     # Не указываем количество компонент, их число будет равно
                # числу признаков
pca = PCA(n_components = 3) # Выделяются всего 3 первых главных компоненты

pca.fit(x)  # Вычисляем новую системе координат, связанную
            #           с главными компонентами
x_tranformed = pca.transform(x) # Вычисляем координаты объектов в
                                # новой системе координат, связанной с
                                # заданным числом главных компонент.

Давайте напишем простую программу в Jupyter-Lab, которая иллюстрирует метод главных компонент в 2-мерном случае.

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

import numpy as np
import random
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
%matplotlib inline

# Generate a random data
origin = np.array([0., 1.])
v = np.array([2., 1.])
v /= np.linalg.norm(v)       # make norm(v) = 1.
n = np.array([-v[1], v[0]]) 
m = 50

x = [
    origin  + v * random.normalvariate(3., 4.) +
    n * random.normalvariate(0., 2.)
    for i in range(m)
]
x = np.array([[xx[0], xx[1]] for xx in x])

plt.axes(aspect="equal")
# plt.xlim(min(x[:, 0])-1., max(x[:, 0])+1)
# plt.ylim(min(x[:, 1])-1., max(x[:, 1])+1)
# plt.grid()

plt.scatter(x[:, 0], x[:, 1], color="blue")

pca = PCA()
pca.fit(x)
e0 = pca.components_[0]
variance0 = pca.explained_variance_[0]
e1 = pca.components_[1]
variance1 = pca.explained_variance_[1]
center = pca.mean_

pca = PCA(n_components = 1)
pca.fit(x)
x_tranformed = pca.transform(x)
x_reduced = pca.inverse_transform(x_tranformed)

plt.scatter(x_reduced[:, 0], x_reduced[:, 1], color="red", alpha=0.3)

plt.arrow(
    center[0], center[1],
    e0[0] * variance0**0.5, e0[1] * variance0**0.5,
    color="black",
    width=0.1
)
plt.arrow(
    center[0], center[1],
    e1[0] * variance1**0.5, e1[1] * variance1**0.5,
    color="black",
    # length_includes_head=True,
    width=0.1
)

В следующем примере мы рассматриваем базу NIST рукописных цифр (это матрицы размера 8x8, размерность объектов 64). Применяя метод главных компонент, мы оставляем сначала только две первых компоненты, затем 3. Полученные точки изображаются разными цветами, соответствующими десяти разным цифрам, в первом случае на плоскости, во втором случае — в трехмерной проекции.

import numpy as np
import random
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
%matplotlib inline

from sklearn.datasets import load_digits
digits = load_digits()
print(digits.data.shape)
# plt.imshow(digits.data[1001].reshape((-1, 8)), cmap='gray')

pca = PCA()
x = digits.data
pca.fit(x)
# print("Variance:")
# print(pca.explained_variance_)
# print("Variance ratio:")
# print(pca.explained_variance_ratio_)
# xx = np.linspace(0, 64, 64)
# plt.plot(xx, np.cumsum(pca.explained_variance_ratio_))

variance_ratio = 0.8
pca = PCA(variance_ratio)   # 80% of total variance must be preserve
pca.fit(x)

print(
    "for variance ratio ", variance_ratio, " we need", 
    len(pca.components_), "components"
)
pca = PCA(n_components = 2)
x = digits.data
pca.fit(x)
x_transformed = pca.transform(x)
plt.scatter(
    x_transformed[:, 0],
    x_transformed[:, 1],
    c = digits.target,
    alpha = 0.5,
    cmap=plt.cm.get_cmap('tab10')
)
plt.xlabel('component 1')
plt.ylabel('component 2')
plt.colorbar();
# print(digits.target[111])
# plt.imshow(digits.data[111].reshape((-1, 8)), cmap='gray')

pca = PCA(n_components = 3)
x = digits.data
pca.fit(x)
x_transformed = pca.transform(x)
print("Variance_ration:", pca.explained_variance_ratio_)
print("Total variance ratio:", sum(pca.explained_variance_ratio_))

fig = plt.figure()
# ax = fig.add_subplot(111, projection='3d')
ax = plt.axes(projection='3d')
p = ax.scatter3D(
    x_transformed[:, 0],
    x_transformed[:, 1],
    x_transformed[:, 2],
    c = digits.target,
    alpha = 0.5,
    cmap=plt.cm.get_cmap('tab10')
);
fig.colorbar(p)

Занятие 7.05.2022
Классификация рукописных цифр с помощью метода опорных векторов

Метод опорных векторов позволяет отделить друг от друга 2 класса точек (binary classification). Строится разделяющая гиперплоскость в (расширенном) пространстве и линейный классификатор

f(x) = ⟨w, x⟩ − b      линейный классификатор
F(x) = K(w, x) − b    ядровый классификатор

Классификатор принимает положительные значения на точках первого класса и отрицательные значения на точках второго класса.

Как быть, когда у нас много классов?

Пусть мы имеем l классов. Тогда можно применить одну из двух стратегий.

1. one-versus-rest ("ovr") — "каждый против остальных".
   Имеем класс c, разбиваем точки всех классов на два подмножества:
   — точки класса c;
   — точки всех остальных классов.
   Для этих двух множеств строим линейный классификатор f_c(x). Всего мы построим l классификаторов
   f₁(x), f₂(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 классификаторов.

   Рассмотрим конкретный вектор x. Применяем все 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 комбинаций.

Для какой комбинации будет наилучший результат?

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")

Результат выполнения этой программы:

score = 0.9805555555555555
predicted: 6  must be: 5
problematic image:

Для базы данных digits определить, какой набор параметров дает наиболее точные предсказания при использовании метода опорных векторов.

Гиперпараметр C:

Всего получаем 18 наборов параметров:

Вся база цифр разбивается на две части: