import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn import preprocessing
%matplotlib inline
!wget --no-check-certificate -O space_ga.csv https://www.dropbox.com/s/t1webrwixdlaacv/space_ga.txt
Линейная регрессия с батч-оптимизацией(2 балла)¶
Рассмотрим случай, когда данных в выборке много. В таких случаях используется стохастическая или батч-оптимизация.
Загрузите данные из файла space_ga.csv и нормализуйте их. Мы будем предсказывать первый столбец по шести остальным. Эти данные получены с выборов в США в 1980 году.
df = pd.read_csv('space_ga.csv', header=None)
y = df[0].values
scaler = preprocessing.StandardScaler().fit(df)
df = pd.DataFrame(data=scaler.transform(df))
X = df[[1, 2, 3, 4, 5, 6]].values
X = np.hstack((np.ones(shape=(X.shape[0])).reshape((X.shape[0], 1)), X))
X
y
Как вы могли заметить, датасет больше предыдущего. На нём мы попробуем батч-оптимизацию.
Измените функцию для минимизации написанную на семинаре так, чтобы на вход они принимала дополнительный параметр — размер батча. Запустите функцию при разных размерах батча. Прокомментируйте результаты.
def loss(X:np.ndarray, y:np.ndarray, w:np.ndarray, batch_size=None):
if batch_size is None or batch_size > y.shape[0]:
X_i, y_i = X, y
else:
index = np.random.choice(y.shape[0], batch_size, replace=False)
X_i, y_i = X[index], y[index]
return np.mean((X_i @ w - y_i)**2) / 2
def grad(X:np.ndarray, y:np.ndarray, w:np.ndarray, batch_size=None):
if batch_size is None or batch_size > y.shape[0]:
X_i, y_i = X, y
else:
index = np.random.choice(y.shape[0], batch_size, replace=False)
X_i, y_i = X[index], y[index]
return np.mean((X_i @ w - y_i)[..., None] * X_i, axis=0)
w = np.ones((X.shape[1]))
w
def log_int_iterator(start, end, step):
i = start
last = None
while i <= end:
if int(i) != last:
last = int(i)
yield last
i *= step
N = 3000
plt.figure(figsize=[20,5])
plt.xscale("log")
plt.xlim([0.9, X.shape[0]])
plt.title("Loss в зависимости от batch_size")
plt.xlabel("Batch size")
plt.ylabel("Loss")
history_y = []
history_x = []
np.random.seed(42)
for i in log_int_iterator(1, X.shape[0], 1.3):
losses = [loss(X,y,w,i) for _ in range(N)]
print(f"Batch_size:{i}\tloss_min:{np.min(losses)}\tmax:{np.max(losses)}\tmean:{np.mean(losses)}\tstd:{np.std(losses)}")
plt.scatter([i] * len(losses), losses)
history_x.append(i)
history_y.append(np.mean(losses))
plt.plot(history_x, history_y)
Усредненный результат функции для минимизации практически не зависит от размера выборки, так как при каждом вычислении высчитывается средний результат для одной записи, который затем усредняется от количества экспериментов.
Loss для отдельного эксперимента сильно зависит от размера батча, при большем размере разброс значений невелик, однако при малых - сильно растет. Это можно объяснить природой данных - некоторые хорошо описывают обобщенное значение выборки, некоторые - плохо и, возможно, явяляются выбросами.
Возможный идеальный результат:¶
from sklearn.linear_model import LinearRegression
reg = LinearRegression(fit_intercept=False).fit(X, y)
print("Weights:", *[np.round(it, 5) for it in reg.coef_], sep='\t')
print("Loss: ", loss(X, y, reg.coef_))
import time
def super_optimizer_and_plotter_steps(learning_rate, X, y,
batch_size=None,
w=None,
eps = 0.0001,
max_iter=10000,
learning_rate_mul=0.99):
if w is None:
w = np.zeros(shape=(X.shape[1]))
# eps *= X.shape[1]
np.random.seed(42)
start = time.time()
w_old = w + eps
history = [loss(X, y, w)]
i = 0
while np.linalg.norm(w_old - w) > eps:
learning_rate *= learning_rate_mul
w_old = w
w = w - learning_rate * grad(X, y, w, batch_size)
history.append(loss(X, y, w))
i += 1
if i > max_iter:
break
print("Batch: {}\tSteps:{}\tLoss: {}\tTime: {}\tW: {weights}".format(batch_size,
i,
np.round(history[-1], 5),
np.round(time.time() - start, 3),
weights=[np.round(it, 3) for it in w]))
# print("Batch:", batch_size, "\tКоличество шагов:", len(history), "\tloss: ", history[-1], "\ttime:", time.time() - start)
#
plt.grid(True)
plt.plot(history)
plt.title("Descent trajectory. Batch: {} Steps: {}\nLoss: {}\nTime: {}".format(batch_size, i, history[-1], time.time() - start))
plt.xlabel("# iteration")
plt.ylabel("Loss")
plt.figure(figsize=[20,20])
i = 0
max_step = None
min_loss = None
max_loss = None
for batch in [None, 3000, 2000, 1000, 500, 200, 100, 50, 10, 7, 5, 2]:
plt.subplot(4, 3, i + 1)
super_optimizer_and_plotter_steps(0.4, X, y, batch)
_, r = plt.xlim()
if max_step is None or r > max_step:
max_step = r
l, r = plt.ylim()
if min_loss is None or l < min_loss:
min_loss = l
if max_loss is None or r > max_loss:
max_loss = r
i += 1
for i in range(i):
plt.subplot(4, 3, i + 1)
# plt.subplot2grid((4,3),(i//3, i%3))
# plt.subplot(431 + j)
plt.xlim([0, max_step])
plt.ylim([0.008, 0.05])
plt.tight_layout()
Масштаб по обоим осям на графиках привен к одному виду.
Итерации с $batch\_size < 50$ сильно колеблются, у них наблюдается увеличение количества шагов обучения, так как шаг в сторону антиградиента сильно зависит от тех элементов, которые были выбраны случайным образом.
В целом, все выборки с $batch\_size > 2$ привели к примерно одинаковым результатам.
Все выборки с $batch\_size >= 1000$ затратили примерно одинаковое количество шагов.
Все итерации с батч оптимизацией затратили больше времени, чем итерации без нее. Это можно объяснить тем, что тратится много времени для генерации случайных чисел и взятии элементов по данным индексам. Возможно, на больших выборках генерация случайных чисел и взятие индекса будет занимать меньше времени, чем перемножение больших матриц.
Двумерная классификация(1 балл)¶
Решим задачу 2D классификации синтетических данных.
with open('train.npy', 'rb') as fin:
X = np.load(fin)
with open('target.npy', 'rb') as fin:
y = np.load(fin)
plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Paired, s=20)
plt.show()
Как можно заметить, данные сверху линейно неразделимы. Поэтому мы должны добавить дополнительные признаки(или использовать нелинейную модель). Можно заметить, что гиперплоскость разделяющая два класса принимает форму круга, поэтому мы можем добавить квадратичные признаки чтобы сделать классы линейно разделимыми.
def expand(X):
"""
Adds quadratic features.
This expansion allows your linear model to make non-linear separation.
For each sample (row in matrix), compute an expanded row:
[feature0, feature1, feature0^2, feature1^2, feature0*feature1, 1]
:param X: matrix of features, shape [n_samples,2]
:returns: expanded features of shape [n_samples,6]
"""
X_expanded = np.zeros((X.shape[0], 6))
X_expanded[:,0] = X[:, 0]
X_expanded[:,1] = X[:, 1]
X_expanded[:,2] = X[:, 0]**2
X_expanded[:,3] = X[:, 1]**2
X_expanded[:,4] = X[:, 0] * X[:,1]
X_expanded[:,5] = np.ones(X.shape[0])
return X_expanded
X_expanded = expand(X)
# simple test on random numbers
dummy_X = np.array([
[0,0],
[1,0],
[2.61,-1.28],
[-0.59,2.1]
])
# call your expand function
dummy_expanded = expand(dummy_X)
# what it should have returned: x0 x1 x0^2 x1^2 x0*x1 1
dummy_expanded_ans = np.array([[ 0. , 0. , 0. , 0. , 0. , 1. ],
[ 1. , 0. , 1. , 0. , 0. , 1. ],
[ 2.61 , -1.28 , 6.8121, 1.6384, -3.3408, 1. ],
[-0.59 , 2.1 , 0.3481, 4.41 , -1.239 , 1. ]])
#tests
assert isinstance(dummy_expanded,np.ndarray), "please make sure you return numpy array"
assert dummy_expanded.shape == dummy_expanded_ans.shape, "please make sure your shape is correct"
assert np.allclose(dummy_expanded,dummy_expanded_ans,1e-3), "Something's out of order with features"
print("Seems legit!")
Логистическая регрессия(3 балла)¶
Для классификации объектов мы будем получать вероятность того что объект принадлежит к классу '1'. Чтобы предсказывать вероятность мы будем использовать вывод линейной модели и логистической функции:
$$ a(x; w) = \langle w, x \rangle $$$$ P( y=1 \; | \; x, \, w) = \dfrac{1}{1 + \exp(- \langle w, x \rangle)} = \sigma(\langle w, x \rangle)$$def probability(X, w):
"""
Given input features and weights
return predicted probabilities of y==1 given x, P(y=1|x), see description above
Don't forget to use expand(X) function (where necessary) in this and subsequent functions.
:param X: feature matrix X of shape [n_samples,6] (expanded)
:param w: weight vector w of shape [6] for each of the expanded features
:returns: an array of predicted probabilities in [0,1] interval.
"""
return 1/(1 + np.exp(-X.dot(w)))
dummy_weights = np.linspace(-1, 1, 6)
ans_part1 = probability(X_expanded[:1, :], dummy_weights)[0]
Для логистической регрессии оптимальное значение весов $w$ находится с помощью минимизации кросс-энтропии:
Loss для одного сэмпла: $$ l(x_i, y_i, w) = - \left[ {y_i \cdot log P(y_i = 1 \, | \, x_i,w) + (1-y_i) \cdot log (1-P(y_i = 1\, | \, x_i,w))}\right] $$
Loss для нескольких сэмплов: $$ L(X, \vec{y}, w) = {1 \over \ell} \sum_{i=1}^\ell l(x_i, y_i, w) $$
def compute_loss(X, y, w):
"""
Given feature matrix X [n_samples,6], target vector [n_samples] of 1/0,
and weight vector w [6], compute scalar loss function L using formula above.
Keep in mind that our loss is averaged over all samples (rows) in X.
"""
n = X.shape[0]
prob = probability(X, w)
return -np.sum(y * np.log(prob) + (1-y) * np.log(1 - prob))/n
Т.к мы обучаем нашу модель с помощью градиентного спуска мы должны считать градиенты. Для этого нам нужны производные функции потерь по каждому из весов.
$$ \nabla_w L = {1 \over \ell} \sum_{i=1}^\ell \nabla_w l(x_i, y_i, w) $$Выведите формулу для подсчета градиента.
def compute_grad(X, y, w):
"""
Given feature matrix X [n_samples,6], target vector [n_samples] of 1/0,
and weight vector w [6], compute vector [6] of derivatives of L over each weights.
Keep in mind that our loss is averaged over all samples (rows) in X.
"""
prob = probability(X, w) - y
return prob.dot(X)/X.shape[0]
Вспомогательная функция для визуализации предсказаний:
from IPython import display
h = 0.01
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
def visualize(X, y, w, history):
"""draws classifier prediction with matplotlib magic"""
Z = probability(expand(np.c_[xx.ravel(), yy.ravel()]), w)
Z = Z.reshape(xx.shape)
plt.subplot(1, 2, 1)
plt.contourf(xx, yy, Z, alpha=0.8)
plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Paired)
plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.subplot(1, 2, 2)
plt.plot(history)
plt.grid()
ymin, ymax = plt.ylim()
plt.ylim(0, ymax)
display.clear_output(wait=True)
plt.show()
visualize(X, y, dummy_weights, [0.5, 0.5, 0.25])
Обучение¶
В данной секции мы будем использовать функции, написанные вами, чтобы обучить наш классификатор с помощью стохастического градиентного спуска.
Mini-batch SGD(1 балл)¶
Стохастический градиентный спуск берет рандомный батч из $m$ сэмплов на каждой итерации, подсчитывает градиент функции потерь на этом батче и делает шаг градиентного спуска:
$$ w_t = w_{t-1} - \eta \dfrac{1}{m} \sum_{j=1}^m \nabla_w l(x_{i_j}, y_{i_j}, w_t) $$np.random.seed(42)
w = np.array([0, 0, 0, 0, 0, 1])
eta= 0.1 # learning rate
n_iter = 100
batch_size = 4
loss = np.zeros(n_iter)
plt.figure(figsize=(12, 5))
for i in range(n_iter):
ind = np.random.choice(X_expanded.shape[0], batch_size)
loss[i] = compute_loss(X_expanded, y, w)
if i % 10 == 0:
visualize(X_expanded[ind, :], y[ind], w, loss)
# Keep in mind that compute_grad already does averaging over batch for you!
# TODO:<your code here>
w = w - eta * compute_grad(X_expanded[ind, :], y[ind], w)
visualize(X, y, w, loss)
plt.clf()
SGD with momentum(1 балл)¶
Momentum это метод позволяющий корректировать шаг SGD в нужное направление и уменьшать осцилляции как показано на рисунке. Данный эффект достигается с помощью добавления предыдущих шагов с коэффициентом $\alpha$ к текущему градиенту для каждого шага с обновлением весов.
np.random.seed(42)
w = np.array([0, 0, 0, 0, 0, 1])
eta = 0.05 # learning rate
alpha = 0.9 # momentum
nu = np.zeros_like(w)
n_iter = 100
batch_size = 4
loss = np.zeros(n_iter)
plt.figure(figsize=(12, 5))
for i in range(n_iter):
ind = np.random.choice(X_expanded.shape[0], batch_size)
loss[i] = compute_loss(X_expanded, y, w)
if i % 10 == 0:
visualize(X_expanded[ind, :], y[ind], w, loss)
nu = nu * alpha + eta * compute_grad(X_expanded[ind, :], y[ind], w)
w = w - nu
visualize(X, y, w, loss)
plt.clf()
ADAM(2 балла)¶
Реализуйте метод ADAM, использующий градиенты и квадраты градиентов сглаженные экспоненциальным скользящим средним:
\begin{eqnarray} m_t &=& \beta_1 m_{t-1} + (1-\beta_1) g_t\\ s_t &=& \beta_2 s_{t-1} + (1-\beta_2) g_t^2 \\ w_t &=& w_{t-1} - \eta \times \frac{\sqrt{ 1 - \beta_2^t}}{ 1 - \beta_1^t} \times \frac{ m_t }{ \sqrt{s_t+eps}} \end{eqnarray}np.random.seed(42)
w = np.array([0, 0, 0, 0, 0, 1.])
m = np.zeros(w.shape)
s = np.zeros(w.shape)
eta = 0.1 # learning rate
beta_1 = 0.9 # moving average of gradient
beta_2 = 0.999 # moving average of gradient norm squared
g2 = None # we start with None so that you can update this value correctly on the first iteration
eps = 1e-8
n_iter = 100
batch_size = 4
loss = np.zeros(n_iter)
plt.figure(figsize=(12,5))
for i in range(n_iter):
ind = np.random.choice(X_expanded.shape[0], batch_size)
loss[i] = compute_loss(X_expanded, y, w)
if i % 10 == 0:
visualize(X_expanded[ind, :], y[ind], w, loss)
# TODO:<your code here>
g = compute_grad(X_expanded[ind, :], y[ind], w)
m = beta_1 * m + (1 - beta_1) * g
s = beta_2 * s + (1 - beta_2) * g**2
w = w - eta * np.sqrt(1 - beta_2 ** (i + 1)) / (1 - beta_1 ** (i + 1)) * m / np.sqrt(s + eps)
visualize(X, y, w, loss)
plt.clf()
