В итоге вы должны получить функцию $f(x) \to [0,1]^{10}$, которая на вход получает картинку с написанной от руки цифрой, а на выход дает 10 вероятностей от 0 до 1 принадлежности к каждому из классов (цифре). Картинка это вектор из 8x8 = 64 чисел. Мы будем рассматривать параметрическое семейство функций $F(c)$, таких, что если $f_с \in F$, то она удовлетворяет нашим требованиям. Кроме того, для каждой функции $f_c \in F$ мы можем посчитать, насколько она хорошо работает на некотором наборе картинок - это будет функционал качества этой функции $\text{loss}(f_c, \text{images})$. Чем он меньше, тем лучше: в идеале $\text{loss}$ будет давать $0$ в том случае, если наша функция на всех картинках, на которых нарисована цифра $i$ выдала вектор, где все числа отличны от $1$ и только на $i$-м месте стоит $1$. Иначе это будет некоторое положительное число, которое тем больше, чем хуже работает классификатор (потеря).
Итак, возьмем функцию $g(c) = \text{loss}(f_c, \text{images})$, и будем ее минимизировать. Если мы найдем глобальный минимум, то научимся максимально качественно решать задачу классификации с помощью того семейства функций, которое мы выбрали. Глобальный минимум мы, конечно, не сможем аналитически найти, поэтому будем решать задачу минимизации методом градиентного спуска.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import math
from sklearn.datasets import load_digits
from sklearn.cross_validation import train_test_split
import tqdm
digits = load_digits()
In [2]:
digits.data.shape
Out[2]:
In [3]:
plt.figure(figsize=(20, 4))
for index, (image, label) in enumerate(zip(digits.data[0:5], digits.target[0:5])):
plt.subplot(1, 5, index + 1)
plt.imshow(np.reshape(image, (8, 8)), cmap=plt.cm.gray)
plt.title('Training: %i\n' % label, fontsize=20)
Это необходимо для того, чтобы понимать, насколько наша функция реально умеет решать поставленную задачу: понимать, где какая цифра. У нас ограниченная выборка - всего 1797 картинок. Но в реальности нарисованных цифр может быть значительно больше! Если даже наша функция безошибочно работает на всех 1797 картинках, но ошибается вне - это плохо. Обычно график обучения должен выглядит примерно так, если зеленое - обучающая выборка, а красное - тестовая.
In [4]:
plt.plot(
[1/math.log(i+2) for i in range(100)],
color='green',
label='train'
)
plt.plot(
[1/math.log(i+2)*1.1+( (i/50.0-1)**2./5. if i>50 else 0.) for i in range(100)],
color='red',
label='test')
plt.xlabel('Gradient descent iteration')
plt.ylabel('Loss')
plt.legend()
plt.show()
То есть с каждым шагом мы уменьшаем наш $\text{loss}$ на обучающей выборке за счет градиентного спуска, но в какой-то момент функция может начать хорошо работать только на обучающей выборке. Этот эффект (в данном примере около 55 итерации по оси $x$) назвается переобучением.
In [5]:
def one_hot(y, n_classes):
# Делаем вектор из 10 координат с 0 везде, кроме правильного ответа.
tmp = np.zeros(
(len(y), n_classes),
dtype=np.uint8
)
tmp[range(len(tmp)), y] = 1
return tmp
x_train, x_test, y_train, y_test = train_test_split(
(digits.data - 8) / 16,
one_hot(digits.target, 10),
test_size=0.33,
random_state=0
)
Для каждой картинки мы хотим найти вектор $(p_0,\ldots,p_{10})$, вероятностей, такой, что $p_i$ - вероятность того, что на картинка цифра $i$.
Реализуя логистическую регрессию, мы хотим приближать вероятности к их настоящему распределению.
Выражение выдает ответ вида $$ W x + b, $$
где x - наш вектор картинки, а результат - числовой вектор размерности 10 с какими-то числами. Для того, чтобы эти числа стали вероятностями от 0 до 1, реализуем функцию $$ \text{softmax}(W, b, x) = \frac{e^{Wx+b}}{\sum(e^{Wx+b})}, $$ и полученные значения будут как раз давать в сумме 1, и ими мы будем приближать вероятности.
Оценивать качество нашей модели будем с помощью кросс-энтропии, см. https://en.wikipedia.org/wiki/Cross_entropy.
Сначала поймем, что $x$ - вектор размерности 64 (у картинки 64 точки), $W$ - матрица 10 на 64, $b$ - вектор размерности 10.
Положим $x'_i = x_i$ для $ i \leqslant 64 $ и $x'_{65} = 1$. Получили вектор $x'$ размерности 65. Положим $W'_{i,j} = W_{i,j}$ для $ i \leqslant 10, j \leqslant 64$ и $W'_{i,65}=b_i$ для $ i \leqslant 10 $.
Таким образом, к вектору $x$ просто дописали 1, а к матрице $W$ просто приписали вектор $b$ справа.
Заметим теперь, что в точности верно равенство: $Wx+b=W'x'$. Теперь забьем на вектор $b$ и будем считать, что у нас есть матрица 10 на 65, элементы которой надо оценить. Далее везде считаем $W' = W$ и $x' = x$.
Градиент считается по формуле: $W_{k+1} = W_k - \eta_k \nabla L(W_k)$, где $\eta_k$ - шаг, а $L$ - функция $\text{loss}$. Значит, нам надо посчитать градиент функции $L$, то есь найти ее частные производные по всем 650 переменным.
Вспомним, как определяется $L$. Обозначим через $y$ вектор вида $(0, \ldots, 0, 1, 0, \ldots, 0)^T$, где 1 на $k$-м месте, где $k - 1$ - цифра, изображенная на исследуемой картинке. Размерность $y$ равна 10.
Тогда $$ L(W) = -y_1 \ln \frac{e^{(Wx)_1}}{\sum_{i=1}^{10} e^{(Wx)_i}} - \ldots -y_{10} \ln \frac{e^{(Wx)_{10}}}{\sum_{i=1}^{10} e^{(Wx)_i}}. $$
Теперь найдем производную по $W_{i,j}$: $$ \frac{dL(W)}{dW_{i,j}} = -y1 \frac{\sum{k=1}^{10} e^{(Wx)_k}}{e^{(Wx)_1}} \cdot \frac{-e^{(Wx)_1} e^{(Wx)_i} xj} {(\sum{i=k}^{10} e^{(Wx)_k})^2}
Упростим немного: $$ \frac{dL(W)}{dW_{i,j}} = \frac{ x_j e^{(Wx)_i} \sum_{k=1}^{10} y_k } {\sum_{i=k}^{10}} -y_i x_j. $$
Упрощая еще сильнее, приходим к окончательному ответу: $$ \frac{dL(W)}{dW_{i,j}} =\left( \frac{e^{(Wx)_i}}{\sum_{i=k}^{10} e^{(Wx)_k}} -y_i \right) x_j . $$ Соответственно, если $i = 65$, то есть дифференцируем по переменным $ W_{1, 65} = b_1, \ldots, W_{10, 65} = b_{10} $, то коэффициент перед скобкой просто 1.
In [6]:
def softmax(W, x):
# В исходной версии было плохо, что экспонента может переполниться.
# Эту проблему решает функция logaddexp.reduce.
p = np.dot(x, W.T)
return np.exp(p - np.logaddexp.reduce(p, axis=1).reshape(-1, 1))
def loss(y, softmax):
# Формула из Википедии по ссылке выше.
return np.mean(-np.sum(y * np.log(softmax), axis=1))
In [7]:
# Считаем средний по всем картинкам градиент.
def gradients(W, x, y):
p = softmax(W, x)
grads = (p - y).T @ x
return grads / x.shape[0]
In [8]:
# Шаг, который мы выбираем для градиентного спуска. Как сказано, он должен быть постоянным.
get_constant = lambda step_size: lambda W, x, y: step_size
In [9]:
def classify(x_train, x_test, y_train, y_test, step_size=get_constant(1), iters=5000):
# Как было замечено выше, W Размера 10 на 64, а b размера 10, но мы приписываем b к W.
W = np.zeros((10, 64))
b = np.zeros(10)
# Для приписывания запишем b как вектор столбец и воспользуемся функцией hstack.
b = b.reshape(b.size, 1)
W = np.hstack([W, b])
# Соответственно, нужно поменять x_train и x_test.
fictious = np.ones((x_train.shape[0], 1))
x_train = np.hstack([x_train, fictious])
fictious = np.ones((x_test.shape[0], 1))
x_test = np.hstack([x_test, fictious])
# Будем записывать потери на каждом шаге спуска.
losses_train = [loss(y_train, softmax(W, x_train))]
losses_test = [loss(y_test, softmax(W, x_test))]
# Собственно, сам спуск.
for i in tqdm.tqdm(range(iters)):
eta = step_size(W, x_train, y_train)
W = W - eta * gradients(W, x_train, y_train)
losses_train.append(loss(y_train, softmax(W, x_train)))
losses_test.append(loss(y_test, softmax(W, x_test)))
# На выходе имеется оптимальное значение W и массивы потерь.
return W, losses_train, losses_test
In [10]:
W, losses_train, losses_test = classify(x_train, x_test, y_train, y_test, iters=5000)
plt.plot(losses_train, color='green', label='train')
plt.plot(losses_test, color='red', label='test')
plt.xlabel('Gradient descent iteration')
plt.ylabel('Loss')
plt.ylim((0, 0.5))
plt.legend()
plt.show()
In [11]:
print(np.argmin(losses_test))
Чтобы все было лучше видно, ограничили по оси $y$. Как и следовало ожидать, на обучающей выборки ошибка только уменьшается и достигает минимума на последнем шаге. При этом ошибка на тестовой выборке почти стабилизируется и остается больше, чем на тестовой. Но на самом деле она немного растет и достигает минимума на 3840-й итерации, просто растет еле заметно.
In [12]:
def armijo(W, x, y, alpha=0.5, beta=0.5):
s = 100
grad = gradients(W, x, y)
dW = -grad # Направление спуска.
loss_1 = loss(y_train, softmax(W + s * dW, x))
loss_0 = loss(y_train, softmax(W, x))
while loss_1 > loss_0 + alpha * s * (grad * dW).sum():
s = beta * s
loss_1 = loss(y_train, softmax(W + s * dW, x))
loss_0 = loss(y_train, softmax(W, x))
return s
In [13]:
W, losses_train_const, losses_test_const = classify(x_train, x_test,
y_train, y_test,
get_constant(100))
W, losses_train_armijo, losses_test_armijo = classify(x_train, x_test,
y_train, y_test,
armijo)
plt.plot(losses_train_const, color='green', label='train, const step')
plt.plot(losses_test_const, color='red', label='tets, const step')
plt.plot(losses_train_armijo, color='blue', label='train, armijo step')
plt.plot(losses_test_armijo, color='black', label='test, armijo step')
plt.xlabel('Gradient descent iteration')
plt.ylabel('Loss')
plt.ylim((0, 1))
plt.legend()
plt.show()
In [14]:
print(np.argmin(losses_test_const))
In [15]:
iters = np.argmin(losses_test_armijo)
print(iters)
Отсюда видно, что метод Армихо не ускорил достижение минимума на обучающей выборке, однако ускорил на тестовой. Причем на тестовой выборке алгоритм работает намного лучше. Минимум ошибки на тестовой выборке достигается на 160-й итерации.
Рассмотрим метод Армихо и 160-ю итерацию, то есть ту итерацию, где достигается минимум потерь. Получим "оптимальную" в каком-то смысле матрицу $W$. Затем посчитаем $\text{softmax}(\text{x_test})$ и возьмем аргмакс по всем картинкам. Это и будет наиболее вероятная цифра для каждой картинки. Эти данные сравним с истинным ответом.
In [16]:
W, losses_train_armijo, losses_test_armijo = classify(x_train, x_test,
y_train, y_test,
armijo,
iters)
In [17]:
# Не забудем к x_test добавить единицу на 65-ю компоненту.
nx_test = np.hstack([x_test, np.ones(x_test.shape[0]).reshape(x_test.shape[0], 1)])
probabilities = softmax(W, nx_test)
recognized = np.argmax(probabilities, axis=1)
answers = np.argmax(y_test, axis=1)
Посчитаем теперь долю ошибок.
In [18]:
print(100 * (recognized != answers).mean())
Получается, алгоритм ошибается только в 4.38% картинок тестовой выборки!
Посмотрим на примеры картинок, где алгоритм работает правильно.
In [19]:
k = 4
plt.figure(figsize=(k * k, k * k))
index = 0
for (image, r, a) in zip(x_test, recognized, answers):
if r != a:
continue
plt.subplot(k, k, index + 1)
plt.imshow(np.reshape(image, (8, 8)), cmap=plt.cm.gray)
plt.title('Digit: %i - Recognized: %i' % (a, r), fontsize=14)
index += 1
if index == k * k:
break
Теперь посмотрим на примеры картинок, где алгоритм ошибся.
In [20]:
k = 4
plt.figure(figsize=(k * k, k * k))
index = 0
for (image, r, a) in zip(x_test, recognized, answers):
if r == a:
continue
plt.subplot(k, k, index + 1)
plt.imshow(np.reshape(image, (8, 8)), cmap=plt.cm.gray)
plt.title('Digit: %i - Recognized: %i' % (a, r), fontsize=14)
index += 1
if index == k * k:
break