Метод главных компонент (PCA)#
Ковариационная матрица#
Рассмотрим подробнее матрицу объекты-признаки.
А точнее следующее её преобразование:
Видим, что это буквально выглядит как ковариация! Только в ковариации надо, чтобы случайные величины (а наши признаки - это по сути случайные величины) имели нулевое среднее. Так давайте тогда перейдём к таким признакам, чтобы их среднее на выборке было ноль.
Внимание! Здесь и далее подразумеваем, что средние значения столбцов в матрице Х равны нулю. Столбец Intercept, естественно, надо будет отбросить совсем - для этого также вычтем среднее из столбца лэйблов \(y_i\). Такая штука называется центрированием данных:
\[ f_k(x) \stackrel{\text{замена}}{\longmapsto} f_k(x) - \frac{1}{N}\sum^N_{i=1} f_k(x_i) \]\[ y_i \stackrel{\text{замена}}{\longmapsto} y_i - \frac{1}{N}\sum^N_{i=1} y_i \]
Тогда в чистом виде, видим ковариационную матрицу
А тогда, из свойств матрицы \(X^\top X\) понимаем, что существует такое преобразование признаков, что в новых «координатах» ковариационная матрица имеет диагональный вид (это спектральное разложение \(X^\top X\)):
То есть, мы всегда можем перейти от текущего набора признаков к некоррелированным друг к другу признакам, что существенно повысит обусловленность задачи (вспомните смысл числа обусловленности).
Более того, собственные числа \(\Lambda\) матрицы \(X^\top X\) (или сингулярные числа \(\sqrt{\Lambda}\) матрицы \(X\)) - будут скоррелированы с важностью того или иного признака (Чем больше дисперсия признака, тем он «важнее»). Таким образом мы можем даже отсекать «неважные» признаки.
В этом и есть суть PCA.
Метод главных компонент (PCA - principal component analysis)#
Как работает преобразование признаков?
По сути, заметим, что:
Определяет поворот в пространстве признаков.
Из нашей цели диагонализации \(X^\top X\), мы хотим видим:
То есть \(Z = X \cdot U\) - действительно определяет нужное нам преобразование к максимально ортогональным друг другу признакам.
На практике это выглядит как
\[X \cdot \mathbf{w} \stackrel{\text{want}}{=} Y\]\[Z \cdot U^\top \cdot \mathbf{w} \stackrel{\text{want}}{=} Y\]\[Z \cdot \widetilde{\mathbf{w}} \stackrel{\text{want}}{=} Y\]
Обрезка признаков
Вспоминая сингулярное разложение,
Возьмём и модифицируем центральную матрицу, убрав самые наименьшие сингулярные числа, оставив \(k\) наибольших. Составим таким образом новую матрицу объекты-признаки:
Теорема:: Эта \(\widetilde{Z}\) наилучшим образом решает задачу о малоранговом приближении
Таким образом, остаётся просто решать задачу с новой матрицей объекты-признаки \(\widetilde{Z} = \widetilde{X} \cdot U\).
Пример PCA в борьбе с лишними признаками#
import numpy as np
from matplotlib import pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error
# Пример с обрезкой лишних признаков
MAX_ORDER = 20 # Возьмём многочлены 19 степени и ниже (первые 20, считая от нуля)
xs = 2*np.random.random(200).reshape((-1, 1)) - 1
noise = 0.1*np.random.normal(size = xs.shape)
ys = -4 + 3*xs - 2*xs**2 + 5*xs**3 + noise
xs_train, xs_test, ys_train, ys_test = train_test_split(xs, ys, test_size=0.5, random_state=13)
poly_regresssor = PolynomialFeatures(degree=MAX_ORDER, include_bias=False) # А вот так просто можно фитить многочленами
X_train = poly_regresssor.fit_transform(xs_train.reshape(-1, 1)) # Или генерировать матрицу X
X_test = poly_regresssor.fit_transform(xs_test.reshape(-1, 1))
shift_x = X_train.mean(axis=0).reshape((1, -1)) # Сдвигаем на нулевое среднее
X_train = X_train - shift_x
X_test = X_test - shift_x
shift_y = ys_train.mean()
ys_train -= shift_y
ys_test -= shift_y
w = np.linalg.pinv(X_train) @ ys_train
y_train_predicted = X_train @ w
print('Train MSE', mean_squared_error(y_train_predicted, ys_train))
y_test_predicted = X_test @ w
print('Test MSE', mean_squared_error(y_test_predicted, ys_test))
Train MSE 0.008937804137053496
Test MSE 0.018415062974522566
import seaborn as sns
sns.barplot(x = np.arange(len(w)), y = w.flatten()) # Отобразим веса на барплоте для ясности
<Axes: >
Видим огромные значения весов - один из явных признаков переобучения.
Это можно «на пальцах» объяснить следующим образом: какой-то признак имеет очень маленькую дисперсию, поэтому соответствующий ей вес должен быть очень большим для «искуственного» учитывания этого признака.
Сделаем PCA и посмотрим, какие признаки с маленькой дисперсией можно обрезать.
from sklearn.decomposition import PCA # За нас уже всё реализовано
pca = PCA(n_components = 5) # Оставим только 10 наиважнейших признака
pca.fit(X_train)
X_train_reduced = pca.transform(X_train) # Преобразование координат
X_test_reduced = pca.transform(X_test)
w = np.linalg.pinv(X_train_reduced) @ ys_train
y_train_predicted = X_train_reduced @ w
print('Train MSE', mean_squared_error(y_train_predicted, ys_train))
y_test_predicted = X_test_reduced @ w
print('Test MSE', mean_squared_error(y_test_predicted, ys_test))
Train MSE 0.01686699427417666
Test MSE 0.01962030834633463
sing_values = pca.singular_values_
sns.barplot(x = np.arange(len(sing_values)), y = sing_values, log = True)
print("Видим, что 5 признак в 10 раз меньше самого важного - можно и обрезать признаки далее")
Видим, что 5 признак в 10 раз меньше самого важного - можно и обрезать признаки далее
plt.scatter(xs, ys)
plt.scatter(xs_test, y_test_predicted + shift_y, c = 'r')
print('Срезка ненужных признаков помогла получить хорошую аппроксимацию')
Срезка ненужных признаков помогла получить хорошую аппроксимацию
PCA в сжатии изображений#
Для плотных матриц требуется хранить \(N^2\) элементов. Аппроксимация ранга r (метод главных компонент) может уменьшить это число O(\(N_r\)) . В качестве примера рассмотрим сжатие изображений.
#!pip install scikit-image
from skimage import data
from skimage.color import rgb2gray
from numpy.linalg import svd
from skimage import img_as_float
cat = rgb2gray(img_as_float(data.chelsea()))
cat.shape
(300, 451)
plt.imshow(cat)
plt.tight_layout()
U, s, Vc = svd(cat, full_matrices=False)
def compress_and_show(k):
os = cat.shape
cat_compressed = U[:,:k] @ np.diag(s[:k]) @ Vc[:k,:]
compression_ratio = 100.0* (k*(os[0] + os[1])+k)/(os[0]*os[1])
plt.title("compression ratio={:.2f}".format(compression_ratio)+"%")
plt.imshow(cat_compressed)
plt.axis('off')
plt.tight_layout()
return cat_compressed
cat_compressed = compress_and_show(100)
cat_compressed.shape
(300, 451)
from ipywidgets import interact, widgets
def on_value_change(b):
compress_and_show(b)
b_slider = widgets.IntSlider(value=90, min=1, max=100, step=1, description='b')
#interact(on_value_change, b=b_slider)
<function __main__.on_value_change(b)>