Математика для машинного обучения, часть 2

Вот Часть-2 - Математика для машинного обучения. Часть-1. Если вы еще не прошли через это, сделайте это прямо сейчас!

СОДЕРЖАНИЕ

• Собственные значения и собственные векторы
• Симметричные матрицы
• Диагонализация
• Полиномиальная интерполяция
• Линейная регрессия методом наименьших квадратов

Собственные значения и собственные векторы

import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la

Определение

Пусть A - квадратная матрица. ненулевой вектор v - это собственный вектор для A с собственным значением λ, если -

Преобразуя уравнение, мы видим, что v является решением однородной системы уравнений (где I - единичная матрица размера n)

Нетривиальные решения существуют, только если матрица A − λI сингулярна, что означает det (A − λI) = 0 . Следовательно, собственные значения матрицы A являются корнями характеристического многочлена.

scipy.linalg.eig

Функция scipy.linalg.eig вычисляет собственные значения и собственные векторы квадратной матрицы A.

Рассмотрим простой пример с диагональной матрицей:

A = np.array([[1,0],[0,-2]])
print(A)
Out:
[[ 1  0]
 [ 0 -2]]

Функция la.eig возвращает кортеж (eigvals,eigvecs), где eigvals - это 1D Numpy-массив из комплексных чисел, дающий собственные значения A, а eigvecs - 2D-массив Numpy с соответствующими собственными векторами в столбцах:

results = la.eig(A)

Собственные значения A:

print(results[0])
Out:
[ 1.+0.j -2.+0.j]

Соответствующие собственные векторы:

print(results[1])
Out:
[[1. 0.]
 [0. 1.]]

Мы можем распаковать кортеж:

eigvals, eigvecs = la.eig(A)
print(eigvals)
Out:
[ 1.+0.j -2.+0.j]
print(eigvecs)
Out:
[[1. 0.]
 [0. 1.]]

Если мы знаем, что собственные значения являются действительными числами (т. Е. Если A является симметричным), то мы можем использовать метод массива Numpy .real для преобразования массив собственных значений действительных чисел:

eigvals = eigvals.real
print(eigvals)
Out:
[ 1. -2.]

Обратите внимание, что положение собственного значения в массиве eigvals соответствует столбцу в eigvecs с его собственным вектором:

lambda1 = eigvals[1]
print(lambda1)
Out:
-2.0
v1 = eigvecs[:,1].reshape(2,1)
print(v1)
Out:
[[0.]
 [1.]]
A @ v1
Out:
array([[ 0.],
       [-2.]])
lambda1 * v1
Out:
array([[-0.],
       [-2.]])

Примеры

Симметричные матрицы

собственные значения симметричной матрицы всегда действительны, а собственные векторы - всегда ортогонально!

Давайте проверим эти факты с помощью случайных матриц:

n = 4
P = np.random.randint(0,10,(n,n))
print(P)
Out:
[[7 0 6 2]
 [9 5 1 3]
 [0 2 2 5]
 [6 8 8 6]]

Создайте симметричную матрицу S = P.PT:

S = P @ P.T
print(S)
Out:
[[ 89  75  22 102]
 [ 75 116  27 120]
 [ 22  27  33  62]
 [102 120  62 200]]

Давайте распаковываем собственные значения и собственные векторы матрицы S:

evals, evecs = la.eig(S)
print(evals)
Out:
[361.75382302+0.j  42.74593101+0.j  26.33718907+0.j   7.16305691+0.j]

Все собственные значения имеют нулевую мнимую часть, поэтому они действительно являются действительными числами:

evals = evals.real
print(evals)
Out:
[361.75382302  42.74593101  26.33718907   7.16305691]

Соответствующие собственные векторы матрицы A:

print(evecs)
Out:
[[-0.42552429 -0.42476765  0.76464379 -0.23199439]
 [-0.50507589 -0.54267519 -0.64193252 -0.19576676]
 [-0.20612674  0.54869183 -0.05515612 -0.80833585]
 [-0.72203822  0.4733005   0.01415338  0.50442752]]

Давайте проверим, ортогональны ли собственные векторы друг другу:

v1 = evecs[:,0] # First column is the first eigenvector
print(v1)
Out:
[-0.42552429 -0.50507589 -0.20612674 -0.72203822]
v2 = evecs[:,1] # Second column is the second eigenvector
print(v2)
Out:
[-0.42476765 -0.54267519  0.54869183  0.4733005]
v1 @ v2
Out:
-1.1102230246251565e-16

Скалярное произведение собственных векторов v1 и v2 равно нулю (указанное выше число очень близко к нулю и связано с ошибками округления в вычислениях), поэтому они ортогональны!

Диагонализация

Квадратная матрица M является диагонализуемой, если она подобна диагональной матрице. Другими словами, M диагонализуема, если существует обратимая матрица P такая, что D - диагональная матрица.

Прекрасный результат в линейной алгебре состоит в том, что квадратная матрица M размера n является диагонализуемой тогда и только тогда, когда M имеет n независимых собственных векторов. Кроме того,

Если столбцы P являются собственными векторами M и D, имеет соответствующие собственные значения по диагонали.

Давайте воспользуемся этим, чтобы построить матрицу с заданными собственными значениями λ1 = 3, λ2 = 1 и собственными векторами.

P = np.array([[1,1],[1,-1]])
print(P)
Out:
[[ 1  1]
 [ 1 -1]]
D = np.diag((3,1))
print(D)
Out:
[[3 0]
 [0 1]]
M = P @ D @ la.inv(P)
print(M)
Out:
[[2. 1.]
 [1. 2.]]

Давайте проверим, что собственные значения M равны 3 и 1:

evals, evecs = la.eig(M)
print(evals)
Out:
[3.+0.j 1.+0.j]

Проверьте собственные векторы:

print(evecs)
Out:
[[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]

Матричные силы

Пусть M - квадратная матрица. Вычислительная мощность M по умножению матриц -

Это дорого с точки зрения вычислений. Вместо этого давайте воспользуемся диагонализацией для более эффективного вычисления M мощности k.

Давайте посчитаем M power 20 в обоих направлениях и сравним время выполнения:

Pinv = la.inv(P)
k = 20
%%timeit
result = M.copy()
for _ in range(1,k):
    result = result @ M
Out:
42.1 µs ± 11.4 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Давайте воспользуемся диагонализацией для тех же вычислений:

%%timeit
P @ D**k @ Pinv
Out:
6.42 µs ± 1.36 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Диагонализация вычисляет M power k намного быстрее!

Приложения

Полиномиальная интерполяция

Полиномиальная интерполяция находит уникальный полином степени n, который проходит через n + 1 точку в плоскости xy .

Например, две точки в плоскости xy определяют линию, а три точки определяют параболу.

Формулировка

Предположим, у нас есть n + 1 точка в плоскости xy -

такие, что все значения x разные (xi ≠ xj для i ≠ j). Общая форма полинома степени n -

Если p (x) является единственным полиномом степени n, который интерполирует все точки, то коэффициенты a0, a1,…, an удовлетворяют следующим уравнениям:

Следовательно вектор коэффициентов:

является уникальным решением линейной системы уравнений -

где X - это матрица Вандермонда, а y - вектор значений y

Примеры

Простая парабола

Приведем простой пример. Мы знаем, что y = x² - единственный полином степени 2, который интерполирует точки (−1,1), (0,0) и (1,1 ). Давайте вычислим полиномиальную интерполяцию этих точек и проверим ожидаемый результат a0 = 0, a1 = 0 и a2 = 1.

Создайте матрицу Вандермонда X с массивом значений x:

x = np.array([-1,0,1])
X = np.column_stack([[1,1,1],x,x**2])
print(X)
Out:
[[ 1 -1  1]
 [ 1  0  0]
 [ 1  1  1]]

Создайте вектор y значений y:

y = np.array([1,0,1]).reshape(3,1)
print(y)
Out:
[[1]
 [0]
 [1]]

Мы ожидаем решения a = [0,0,1] T:

a = la.solve(X,y)
print(a)
Out:
[[0.]
 [0.]
 [1.]]

Успех!

Другая парабола

Полиномиальная интерполяция трех точек (x0, y0), (x1, y1) и (x2, y2) - это парабола p (x) = a0 + a1x + a2x2 такие, что коэффициенты удовлетворяют

Найдем полиномиальную интерполяцию точек (0,6), (3,1) и (8,2).

Создайте матрицу Вандермонда X:

x = np.array([0,3,8])
X = np.column_stack([[1,1,1],x,x**2])
print(X)
Out:
[[ 1  0  0]
 [ 1  3  9]
 [ 1  8 64]]

И вектор значений y:

y = np.array([6,1,2]).reshape(3,1)
print(y)
Out:
[[6]
 [1]
 [2]]

Вычислить вектор a коэффициентов:

a = la.solve(X,y)
print(a)
Out:
[[ 6.        ]
 [-2.36666667]
 [ 0.23333333]]

И нанесите на график результат:

xs = np.linspace(0,8,20)
ys = a[0] + a[1]*xs + a[2]*xs**2
plt.plot(xs,ys,x,y,'b.',ms=20)
plt.show()
Out:

Подгонка 10 случайных точек

Теперь давайте интерполируем точки с помощью xi = i, i = 0,…, 9 и 10 случайных целых чисел, выбранных из [0,10) как значения y:

N = 10
x = np.arange(0,N)
y = np.random.randint(0,10,N)
plt.plot(x,y,'r.')
plt.show()
Out:

Создайте матрицу Вандермонда и проверьте первые 5 строк и столбцов:

X = np.column_stack([x**k for k in range(0,N)])
print(X[:5,:5])
Out:
[[  1   0   0   0   0]
 [  1   1   1   1   1]
 [  1   2   4   8  16]
 [  1   3   9  27  81]
 [  1   4  16  64 256]]

Мы также можем использовать функцию Numpy numpy.vander.

Мы указываем параметр increasing=True так, чтобы степени xi увеличивались слева направо:

X = np.vander(x,increasing=True)
print(X[:5,:5])
Out:
[[  1   0   0   0   0]
 [  1   1   1   1   1]
 [  1   2   4   8  16]
 [  1   3   9  27  81]
 [  1   4  16  64 256]]

Решите линейную систему:

a = la.solve(X,y)

Постройте интерполяцию:

xs = np.linspace(0,N-1,200)
ys = sum([a[k]*xs**k for k in range(0,N)])
plt.plot(x,y,'r.',xs,ys)
plt.show()
Out:

Успех! Но обратите внимание, насколько нестабильна кривая.

Вот почему лучше использовать кубический сплайн для интерполяции большого количества точек.

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

Линейная регрессия методом наименьших квадратов

Предположим, у нас есть n + 1 балл

в плоскости xy, и мы хотим подогнать линию

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

В линейной регрессии методом наименьших квадратов мы хотим минимизировать сумму квадратов ошибок.

Формулировка

Если мы сформируем матрицы

тогда сумма квадратов ошибок может быть выражена как

Теорема. (Линейная регрессия методом наименьших квадратов). Учитывайте n + 1 балл

в плоскости xy. Коэффициенты a = [a0, a1] T , которые минимизируют сумма квадратов ошибок -

это уникальное решение системы

Эскиз подтверждения. Товар Xa находится в пространстве столбца X. Линия, соединяющая y с ближайшей точкой в ​​пространстве столбца X, перпендикулярна пространству столбца X. Следовательно,

и так

Примеры

Поддельные линейные данные с шумом

Приведем пример с поддельными данными. Давайте построим набор случайных точек на основе модели

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

Мы можем генерировать случайные числа, взятые из стандартного нормального распределения, используя функцию Numpy numpy.random.rand.

Цель состоит в том, чтобы продемонстрировать, что мы можем использовать линейную регрессию для извлечения коэффициентов a0 и a1 из расчета линейной регрессии.

a0 = 2
a1 = 3
N = 100
x = np.random.rand(100)
noise = 0.1*np.random.randn(100)
y = a0 + a1*x + noise
plt.scatter(x,y);
plt.show()
Out:

Давайте воспользуемся линейной регрессией, чтобы получить коэффициенты a0 и a1.

Построить матрицу X:

X = np.column_stack([np.ones(N),x])
print(X.shape)
Out:
(100, 2)

Давайте посмотрим на первые 5 строк X, чтобы убедиться, что они имеют правильную форму:

X[:5,:]
Out:
array([[1.        , 0.92365627],
       [1.        , 0.78757973],
       [1.        , 0.51506055],
       [1.        , 0.51540875],
       [1.        , 0.86563343]])

Используйте scipy.linalg.solve , чтобы найти a -

a = la.solve(X.T @ X, X.T @ y)
print(a)
Out:
[2.02783873 2.95308228]

Мы почти точно получили коэффициенты модели!

Давайте изобразим случайные точки данных с помощью только что вычисленной линейной регрессии.

xs = np.linspace(0,1,10)
ys = a[0] + a[1]*xs
plt.plot(xs,ys,'r',linewidth=4)
plt.scatter(x,y);
plt.show()
Out:

Настоящие данные Коби Брайанта

Для всех любителей баскетбола - поработаем с реальными данными. Коби Брайант завершил карьеру в 2016 году, набрав 33643 очков, что является третьим по величине общим количеством очков в истории НБА. Сколько еще лет нужно было сыграть Коби Брайанту, чтобы передать рекорд Карим Абдул-Джаббар в 38387 очков?

Покойся с миром, Коби Брайант: ’(Вы были лучшими! ❤

Пиком Коби Брайанта стал сезон НБА 2005–2006 годов. Давайте посмотрим на общее количество сыгранных игр Коби Брайанта и количество очков за игру с 2006 по 2016 годы.

years = np.array([2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016])
games = [80,77,82,82,73,82,58,78,6,35,66]
points = np.array([35.4,31.6,28.3,26.8,27,25.3,27.9,27.3,13.8,22.3,17.6])
fig = plt.figure(figsize=(12,10))
axs = fig.subplots(2,1,sharex=True)
axs[0].plot(years,points,'b.',ms=15)
axs[0].set_title('Kobe Bryant, Points per Game')
axs[0].set_ylim([0,40])
axs[0].grid(True)
axs[1].bar(years,games)
axs[1].set_title('Kobe Bryant, Games Played')
axs[1].set_ylim([0,100])
axs[1].grid(True)
plt.show()
Out:

Коби был травмирован большую часть сезона НБА 2013–2014 годов и сыграл всего 6 игр.

Это выброс, поэтому мы можем отбросить эту точку данных:

years = np.array([2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2015, 2016])
games = np.array([80,77,82,82,73,82,58,78,35,66])
points = np.array([35.4,31.6,28.3,26.8,27,25.3,27.9,27.3,22.3,17.6])

Давайте посчитаем среднее количество игр, сыгранных за сезон за этот период:

avg_games_per_year = np.mean(games)
print(avg_games_per_year)
Out:
71.3

Вычислите линейную модель для очков за игру:

X = np.column_stack([np.ones(len(years)),years])
a = la.solve(X.T @ X, X.T @ points)
model = a[0] + a[1]*years
plt.plot(years,model,years,points,'b.',ms=15)
plt.title('Kobe Bryant, Points per Game')
plt.ylim([0,40])
plt.grid(True)
plt.show()
Out:

Теперь мы можем экстраполировать на будущие годы и умножить очки за игры на игры за сезон и вычислить совокупную сумму, чтобы увидеть общее количество очков Коби:

future_years = np.array([2017,2018,2019,2020,2021])
future_points = (a[0] + a[1]*future_years)*avg_games_per_year
total_points = 33643 + np.cumsum(future_points)
kareem = 38387*np.ones(len(future_years))
plt.plot(future_years,total_points,future_years,kareem)
plt.grid(True)
plt.xticks(future_years)
plt.title('Kobe Bryant Total Points Prediction')
plt.show()
Out:

Осталось всего 4 года!

Скучаю по тебе, Коби ❤

Это все, что касается Части-2! Я знаю, что сразу нужно принять во внимание многое! Но ты дожил до конца! Престижность на этом! Не забудьте ознакомиться с другими частями этой статьи - Часть 3, Часть 4 и Часть 5!

Получите доступ к экспертному обзору - Подпишитесь на DDI Intel

Дополнительные ресурсы и ссылки

Есть много других хороших ресурсов, если вы все еще заинтересованы в получении максимальной отдачи от этой темы -

Для полной реализации ознакомьтесь с моим репозиторием GitHub -