В данной заметке подйет речь о преобразования над самими признаками, т.е. над матрицей X. Значения каждого из прецендента нам не важно. Постараемся из X выжать как можно больше данных.
Это предварительная версия! Любые замечания приветсвуются.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Формирование данных на приямой, которая под углом.
# Сгенерируем данные
x = np.random.uniform(-10, 10, 30) # Раномерная выборка на отрезке [-10, 10] в количестве 30 элементов.
x.shape
X = np.row_stack( (x, np.zeros_like(x) ) ) # Две строчки -- к каждому числу добавлена вторая координата (равная 0).
X # Выводим матрицу.
plt.plot(X[0,:], X[1,:], ".") # Нарисуем их
# Они находятся на прямой.
# Теперь подкрутим данные. Разместим их на под уголом.
theta = np.radians( 30 ) # Градусы в радианы.
c, s = np.cos(theta), np.sin(theta) # Поэлементвное присвоение,
# Фактически эквивалентно:
# c = np.cos(theta)
# s = np.sin(theta)
XX = R.dot(X) + T # Применяем к нашим данным вращение и трансляцию.
XXt = XX.transpose()
XXt
plt.plot(XXt[:,0], XXt[:,1], "*", 0, 0, ".")
# Данные теперь под углом. Ордината (красная точка) находится сбоку.
Метод главных компонент (PCA), данные на прямой.
Как восстановить их? Как устранить факт того, что они под уголом? Для этого есть метод главных компонент.
Он позволяет снизить размерность данных, т.е. найти подпространство, которое содержит данные. Самое пространство и преобразование линейны.
# Загружаем метод главных компонент (PCA -- Principal component analysis):
from sklearn.decomposition import PCA
# PCA это конструктор объектов. Для выполнения вычислений нужно объект создать.
# Самое главное задать требуемое размер пространства.
pca = PCA( n_components = 1 ) # Хотим преобзовать данные к прямой, т.е. сделать их одномерными (как были исходные).
pca.fit( XXt ) # Запускаем метод pca на данных XXt, т.е. ишим подплоскость (смещение и поворот) содержащую данные.
# Выводятся используемые значения параметров метода.
XXtd = pca.transform( XXt ) # Применяем найденное (сохраненное в объекте pca) преобразование к данным.
XXtd
pca.explained_variance_ratio_ # Хорошо-ли данные коррелируют с вычисленным преобзованием?
# Фактически говорит о том насколько хорошо удалось снизить размерность,
# много ли данных оказалось искаженными (спроецированными с сильным смещением)
# 1 -- 100% по данных по первой оси (она же единственная).
# Теперь нужно сравнить восстановленные данные с исходными. Для этого их нужно привести к нормальному виду:
from sklearn import preprocessing # preprocessing -- первичные преобразования данных.
# Будем применять preprocessing.scale -- делает выборочное среднее равно нулю, а дисперсию равной 1.
# Иногда, для корретной работы примера, требуется:
# XXtd = -XXtd # см ниже.
XXtds = preprocessing.scale( XXtd[:, 0] ) # Нормируем. Используем функию scale из модуля sklearn.preprocessing.
XXtds
xs = preprocessing.scale( x ) # x -- это исходные даные. Их тоже нормируем.
xs
# Вычислим расстояние между векторами, т.е. отнормированными данными.
np.linalg.norm((xs - XXtds)) # Вызываем функцию из модуля, без from...
# Длина фактически равна нулю. Значит данные были восстановлены полностью.
# Если норма не равна нулю, то требуется массив отразить. см выше: XXtd = -XXtd. Нужно убрать комментарий.
# Что будет если на две оси спроецировать?
pca = PCA( n_components = 2 ) # Настраиваем объект на двумерную плоскость.
XXtd = pca.fit_transform( XXt ) # Настроим параметры модели и сразу преобразуем данные к двум компонентам.
XXtd # Видим, что даные состоят из двумерных точек. Вторая компонента фактически равна нулю (машинное эпсилон).
# Построим преобращованные (с целью снижения размерность) точки.
plt.plot(XXtd[:,0], XXtd[:,1], "*" ) # Но график будет обманчевым. Ось по y отмасштабирована (разные масштабы)!
plt.plot(XXtd[:,0], XXtd[:,1], "*" )
plt.axis('equal') # Деламем оси одинакового масштаба.
# Правильнее так. Теперь видно, что данные расположены вдоль оси x.
pca.explained_variance_ratio_ # Теперь тоже проверяем насколько хорошо данные коррелируют с осями.
# Видно, что первая ось все содержит. Для второй оси коеффициент равен 0.
Видно, что увеличиение количества осей мы не теряем информацию о том, что данные на самом деле одномерные.
Метод главных компонент (PCA), данные на прямой плюс шум
# Добавим шума (к y, раньше он был раен 0)
X[1,:] = np.random.randn( (30) )/4 # вторая компонента содержит номальный шум, раньше там был 0.
XX = R.dot(X) + T # Как и раньше применяем вращение и трансляцию.
XXt = XX.transpose() # Для удобства. Например, при вызове plot, fit_transform.
plt.plot(XXt[:,0], XXt[:,1], "*", 0, 0, ".")
pca = PCA( n_components = 1 ) # Настраиваем объект на одну компоненту.
XXtd = pca.fit_transform( XXt ) # Ищим параметры преобразования и преобразуем.
pca.explained_variance_ratio_ # Теперь важность первой оси немного упала. Но совсем немного!
# Теперь спроецируем на две оси.
pca = PCA( n_components = 2 )
XXtd = pca.fit_transform( XXt )
XXtd # Видно, что вторая компонента фактически равна шуму.
plt.plot(XXtd[:,0], XXtd[:,1], "*" ) # Но график как и раньше будет обманный. Ось по y отмасштабирована!
plt.plot(XXtd[:,0], XXtd[:,1], "*" ) # Правильнее так. Теперь видно, что данные расположены вдоль оси x.
plt.axis('equal') # Деламем оси одинакового масштаба. Вдоль y расположен шум.
pca.explained_variance_ratio_ # Теперь тоже проверяем насколько хорошо данные коррелируют с осями.
# Видно, что первая ось содержит почти все. Для второй оси коеффициент почти равен 0, но не 0.
XXtd[:,1] # Это шум.
import scipy.stats as models
models.norm.fit( XXtd[:,1] )[1], models.norm.fit( X[1,:] )[1]
# Почти как исходное 1/4 = 0.25. У исходного и преобразованного оценка почти совпадает.
Аналогичный пример в 3d
данные находятся в плоскости.
from sklearn.decomposition import PCA
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
e = np.exp(1)
np.random.seed(4)
def pdf(x):
return 0.5 * (stats.norm(scale=0.25 / e).pdf(x)
+ stats.norm(scale=4 / e).pdf(x))
y = np.random.normal(scale=0.5, size=(30000))
x = np.random.normal(scale=0.5, size=(30000))
z = np.random.normal(scale=0.1, size=len(x))
density = pdf(x) * pdf(y)
pdf_z = pdf(5 * z)
density *= pdf_z
a = x + y
b = 2 * y
c = a - b + z
norm = np.sqrt(a.var() + b.var())
a /= norm
b /= norm
def plot_figs(fig_num, elev, azim):
fig = plt.figure(fig_num, figsize=(8, 6))
plt.clf()
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=elev, azim=azim)
ax.scatter(a[::10], b[::10], c[::10], c=density[::10], marker='+', alpha=.4)
Y = np.c_[a, b, c]
# Using SciPy's SVD, this would be:
# _, pca_score, V = scipy.linalg.svd(Y, full_matrices=False)
pca = PCA(n_components=3)
pca.fit(Y)
pca_score = pca.explained_variance_ratio_
V = pca.components_
x_pca_axis, y_pca_axis, z_pca_axis = V.T * pca_score / pca_score.min()
x_pca_axis, y_pca_axis, z_pca_axis = 3 * V.T
x_pca_plane = np.r_[x_pca_axis[:2], - x_pca_axis[1::-1]]
y_pca_plane = np.r_[y_pca_axis[:2], - y_pca_axis[1::-1]]
z_pca_plane = np.r_[z_pca_axis[:2], - z_pca_axis[1::-1]]
x_pca_plane.shape = (2, 2)
y_pca_plane.shape = (2, 2)
z_pca_plane.shape = (2, 2)
ax.plot_surface(x_pca_plane, y_pca_plane, z_pca_plane)
ax.w_xaxis.set_ticklabels([])
ax.w_yaxis.set_ticklabels([])
ax.w_zaxis.set_ticklabels([])
elev = -40
azim = -80
plot_figs(1, elev, azim)
elev = 30
azim = 20
plot_figs(2, elev, azim)
Другие примеры с PCA
X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
X
plt.plot(X[:,0], X[:,1], "*")
plt.axis('equal')
pca = PCA(n_components=2)
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.singular_values_)
X2 = pca.transform(X)
X2
plt.plot(X2[:,0], X2[:,1], "*")
#plt.axis('equal')
pca_full = PCA(n_components=2, svd_solver='full')
pca_full.fit( X )
pca_arpack = PCA(n_components=1, svd_solver='arpack') #долюно быть n_components < n_features
pca_arpack.fit( X )
X_arpack=pca_arpack.transform(X)
X_arpack
plt.plot(np.sort(X_arpack[:,0]), "*")
from sklearn.decomposition import KernelPCA
from sklearn.datasets import make_circles
# Для обеспечения повтаряемости эксперимента
np.random.seed(0)
# Создадим "кривые" данные, две концентрические окружности с центром в точке (0, 0)
X, y = make_circles(n_samples=400, factor=.3, noise=.05)
X.shape, y.shape
np.unique(y)
plt.title("Исходное пространство")
reds = y == 0
blues = y == 1
plt.scatter(X[reds, 0], X[reds, 1], c="red",
s=20, edgecolor='k')
plt.scatter(X[blues, 0], X[blues, 1], c="blue",
s=20, edgecolor='k')
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
plt.scatter(X_pca[reds, 0], X_pca[reds, 1], c="red",
s=20, edgecolor='k')
plt.scatter(X_pca[blues, 0], X_pca[blues, 1], c="blue",
s=20, edgecolor='k')
plt.title("После МГК преобразования")
plt.xlabel("1ая главная компонента")
plt.ylabel("2ая главная компонента")
В общем ничего не поменялось. Теперь пробуем ядровый метод главных компонент.
kpca = KernelPCA(n_components=2, kernel="rbf", fit_inverse_transform=True, gamma=10) # Ядро rbf, radial basis function
X_kpca = kpca.fit_transform(X)
X_back = kpca.inverse_transform(X_kpca)
plt.scatter(X_kpca[reds, 0], X_kpca[reds, 1], c="red",
s=20, edgecolor='orange')
plt.scatter(X_kpca[blues, 0], X_kpca[blues, 1], c="blue",
s=20, edgecolor='c')
plt.title("После Ядровый МГК преобразования")
plt.xlabel("1ая главная компонента")
plt.ylabel("2ая главная компонента")
Исходные два множества четко разделились.
Начертим кривые координаты.
l = np.linspace(-1.5, 1.5, 5)
l
grid = np.meshgrid( l, l )
grid
grid[0][1,3],grid[1][1,3]
type(grid), len(grid) #! Это список, а не массив !
grid[0]
grid[1]
np.ravel([[1, 2],[3,4]]) # выпрямляет массив в одномерный, схоже с flatten
plt.title("Кривые координаты в исходном пространстве")
reds = y == 0
blues = y == 1
plt.scatter(X[reds, 0], X[reds, 1], c="red",
s=20, edgecolor='k')
plt.scatter(X[blues, 0], X[blues, 1], c="blue",
s=20, edgecolor='k')
X1, X2 = np.meshgrid(np.linspace(-1.0, 1.0, 50), np.linspace(-1.0, 1.0, 50))
X_grid = np.array([np.ravel(X1), np.ravel(X2)]).T
# projection on the first principal component (in the phi space)
Z_grid = kpca.transform(X_grid)[:, 0].reshape(X1.shape)
plt.contour(X1, X2, Z_grid, colors='grey', linewidths=1, origin='lower')
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.scatter(X_back[reds, 0], X_back[reds, 1], c="red",
s=20, edgecolor='k')
plt.scatter(X_back[blues, 0], X_back[blues, 1], c="blue",
s=20, edgecolor='k')
plt.title("Обратно в исходном пространстве")
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
# Загрузка данных из самой бибилотеки.
from sklearn import datasets # В dataset есть нексолько традиционых наборов.
iris = datasets.load_iris()
X = iris.data # Сами данные
y = iris.target # Класс, в данной заметке он нужен нам только для визуализации.
target_names = iris.target_names # Имена классов.
X.shape, y.shape, target_names.shape # Покажем чему раны размеры.
Инача говоря данные содержат 150 примеров. Для каждого примера заданы 4 признака, а также разметку (чем пример является). Для числовых обозначений классов даны соответствующие текстовые названия.
X # Данные, четыре признака. 150 примера (прецендента).
y # Класса. Фактически сначала идут преценденты относящиеся к первому классу, потом второму, а уже потом третьему.
target_names # Имена классов, т.е. 0 классу соответствует имя setosa, ..., 2 -- virginica
# Да, тип элементов может быть текстом!
np.unique(y) # Найдем номера различающихся классов. Функция ищет уникальные элементы массива.
# Напомню.
y * 2 + 1 # Выполняет данное вычисление поэлементно.
# По аналогии.
y==1 # Применили проверку на равенство 1 для каждого из элементов массива. Истина -- True, ложь -- False.
# Напечатам размер каждого из классов, т.е. количество прецендентов в каждом классе.
# Далее идет так называемый цикл, конкретно цикл for. Он выполняет некий код для каждого из элементов по одному.
for i in np.unique(y): # i пробегает список элементов [0,1,2], т.е. принимает каждое из данных значений по отдельности.
# Тело начинается с отступа. Замечу, в предыдущей строчке в конце был символ :!
print( "размер класса ", i, sum( y==i ) ) # Для каждого из жлементов выполянется данный код. Идет печать.
# В арифметических операциях истина (True) равна 1, а лож (False).
# Видим, что классы имеют один и тот же размер, равный 50.
nums = np.unique(y) # Фактически это [0, 1, 2]
list( zip( nums, target_names ) ) # Списки сшиваются по-элементно вместе. Как молния с бегунком у куртки.
# nums[0] и target_names[0] в первый элемент нового списка... (0, 'setosa') ... и так далее
# по аналогии все далается при большем количестве списоков. Рзамеры их должны сопадать!
X_1 = X[ y == 1 ] # Выбираем только те элементы для которых условие истинно, т.е. относящиеся к классу 1.
X_1
X_1.shape, sum(y==1) # Перепроверяем размеры.
X[y == 1, 0] # Берем массив нулевых признаков отсносящихся к классу 1. Фактически 0ой вектор-столбец из X_1.
# Иначе говоря, идет проекция всего признакового пространства (там было 4 признака) на 0 ось.
colors = ['navy', 'turquoise', 'darkorange']
lw = 2
for color, i, target_name in zip(colors, [0, 1, 2], target_names): # Хитрым образом записаный цикл for.
plt.scatter(X[y == i, 0], X[y == i, 1], color=color, alpha=.8, lw=lw, # Отрисовываем каждый класс по отдельности.
label=target_name) # Названия наборов данных.
plt.legend(loc='best', shadow=False, scatterpoints = 1 ) # Отображение легенды.
plt.title('Исходные IRIS данные') # Название графика.
Синий класс можно почти отделить, а вот оранжевый и берюзовый скорее нет. Как быть? Приобразуем систему координат так, что бы было видно. Обращаю внимание: исходные даные 4х мерные, а картинка 2х мерная.
Упр. Изучить другие проекции.
pca = PCA(n_components=2, svd_solver='full')
X_r = pca.fit_transform(X) # Сразу строим модель и преобразуем.
for color, i, target_name in zip(colors, [0, 1, 2], target_names):
plt.scatter(X_r[y == i, 0], X_r[y == i, 1], color=color, alpha=.8, lw=lw,
label=target_name)
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('PCA of IRIS dataset')
Теперь данные гораздо проще отделить данные друг отдруга. По крайней мере можно провести четкую границу/ломанную между классами.
Посмотрим на все это дело в 3d
from mpl_toolkits.mplot3d import Axes3D
from sklearn import decomposition
from sklearn import datasets
np.random.seed(5)
iris = datasets.load_iris()
X = iris.data
y = iris.target
np.choose([0, 2, 1, 2], ['a','b','c']) # Один список используется для индексации элементов из второго списка
fig = plt.figure(figsize=(8, 6))
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=20, azim=60)
for name, label in [('Setosa', 0), ('Versicolour', 1), ('Virginica', 2)]:
ax.text3D(X[y == label, 0].mean(),
X[y == label, 1].mean() + 1.0,
X[y == label, 2].mean(), name,
horizontalalignment='center')
# Отрисовываем данные.
yy = np.choose(y, [1, 2, 0]).astype(np.float)
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=yy, cmap=plt.cm.spectral, edgecolor='k')
pca = decomposition.PCA(n_components=3)
pca.fit(X)
Xpca = pca.transform(X)
fig = plt.figure(figsize=(8, 6))
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=20, azim=90)
for name, label in [('Setosa', 0), ('Versicolour', 1), ('Virginica', 2)]:
ax.text3D(Xpca[y == label, 0].mean(),
Xpca[y == label, 1].mean() + 1.0,
Xpca[y == label, 2].mean(), name,
horizontalalignment='center')
# Отрисовываем данные.
yy = np.choose(y, [1, 2, 0]).astype(np.float)
ax.scatter(Xpca[:, 0], Xpca[:, 1], Xpca[:, 2], c=yy, cmap=plt.cm.spectral, edgecolor='k')
# В Методе главных компонент считалось, что разброс энергии (шума) постоянен по всем осям.
# В данном методе каждая из ортогональных осей независима.
from sklearn.decomposition import FactorAnalysis
fa = FactorAnalysis( n_components=2 )
X_fa = fa.fit_transform(X)
for color, i, target_name in zip(colors, [0, 1, 2], target_names):
plt.scatter(X_fa[y == i, 0], X_fa[y == i, 1], color=color, alpha=.8, lw=lw,
label=target_name)
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('FA of IRIS dataset')
Когда будет разница видна?
В методе гоавных компонентов оси ортонормированы. Это жесткое требование. от него можно отказатся применив следующий подход.
Еще более общий случай метода главных компонент
# Теперь оси необязательно ортогональны
from sklearn.decomposition import FastICA
# ?rng
# Создадим данные
rng = np.random.RandomState(42)
S = rng.standard_t(1.5, size=(20000, 2)) #standard_t 1.5, size=(20000, 2) randn
S[:, 0] *= 2.
S, S.shape
SS = S / S.std()
plt.scatter(SS[:, 0], SS[:, 1], s=2, marker='o', zorder=10,
color='steelblue', alpha=0.5)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
# Смешиваем данные распределенные по разным осям
A = np.array([[1, 1], [0, 2]]) # Матрица смешивания
X = np.dot(S, A.T) # Формируем наблюдения
XX = X / X.std()
plt.scatter(XX[:, 0], XX[:, 1], s=2, marker='o', zorder=10,
color='steelblue', alpha=0.5)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
pca = FactorAnalysis() #PCA()
S_pca_ = pca.fit(X).transform(X)
ica = FastICA(random_state=rng)
S_ica_ = ica.fit(X).transform(X) # Estimate the sources
S_ica_ /= S_ica_.std(axis=0)
def plot_samples(S, axis_list=None):
plt.scatter(S[:, 0], S[:, 1], s=2, marker='o', zorder=10,
color='steelblue', alpha=0.5)
if axis_list is not None:
colors = ['orange', 'red']
for color, axis in zip(colors, axis_list):
axis /= axis.std()
x_axis, y_axis = axis
# Trick to get legend to work
plt.plot(0.1 * x_axis, 0.1 * y_axis, linewidth=2, color=color)
plt.quiver(0, 0, x_axis, y_axis, zorder=11, width=0.01, scale=6,
color=color)
plt.hlines(0, -3, 3)
plt.vlines(0, -3, 3)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.xlabel('x')
plt.ylabel('y')
plt.figure(figsize=(8, 6))
plt.subplot(2, 2, 1)
plot_samples(S / S.std())
plt.title('True Independent Sources')
axis_list = [pca.components_.T, ica.mixing_]
plt.subplot(2, 2, 2)
plot_samples(X / np.std(X), axis_list=axis_list)
legend = plt.legend(['PCA', 'ICA'], loc='upper right')
legend.set_zorder(100)
plt.title('Observations')
plt.subplot(2, 2, 3)
plot_samples(S_pca_ / np.std(S_pca_, axis=0))
plt.title('PCA recovered signals')
plt.subplot(2, 2, 4)
plot_samples(S_ica_ / np.std(S_ica_))
plt.title('ICA recovered signals')
from scipy import signal
from sklearn.decomposition import FastICA, PCA
# #############################################################################
# Generate sample data
np.random.seed(0)
n_samples = 2000
time = np.linspace(0, 8, n_samples)
s1 = np.sin(2 * time) # Signal 1 : sinusoidal signal
s2 = np.sign(np.sin(3 * time)) # Signal 2 : square signal
s3 = signal.sawtooth(2 * np.pi * time) # Signal 3: saw tooth signal
S = np.c_[s1, s2, s3]
S += 0.2 * np.random.normal(size=S.shape) # Add noise
S /= S.std(axis=0) # Standardize data
# Mix data
A = np.array([[1, 1, 1], [0.5, 2, 1.0], [1.5, 1.0, 2.0]]) # Mixing matrix
X = np.dot(S, A.T) # Generate observations
# Compute ICA
ica = FastICA(n_components=3)
S_ = ica.fit_transform(X) # Reconstruct signals
A_ = ica.mixing_ # Get estimated mixing matrix
# We can `prove` that the ICA model applies by reverting the unmixing.
assert np.allclose(X, np.dot(S_, A_.T) + ica.mean_)
# For comparison, compute PCA
pca = PCA(n_components=3)
H = pca.fit_transform(X) # Reconstruct signals based on orthogonal components
# #############################################################################
# Plot results
plt.figure(figsize=(12,12))
models = [X, S, S_, H]
names = ['Observations (mixed signal)',
'True Sources',
'ICA recovered signals',
'PCA recovered signals']
colors = ['red', 'steelblue', 'orange']
for ii, (model, name) in enumerate(zip(models, names), 1):
plt.subplot(4, 1, ii)
plt.title(name)
for sig, color in zip(model.T, colors):
plt.plot(sig, color=color)
X.shape
A_/(-46)
from time import time
import scipy as sp
from sklearn.decomposition import MiniBatchDictionaryLearning
from sklearn.feature_extraction.image import extract_patches_2d
from sklearn.feature_extraction.image import reconstruct_from_patches_2d
try: # SciPy >= 0.16 have face in misc
from scipy.misc import face
face = face(gray=True)
except ImportError:
face = sp.face(gray=True)
# Convert from uint8 representation with values between 0 and 255 to
# a floating point representation with values between 0 and 1.
face = face / 255.
# downsample for higher speed
face = face[::2, ::2] + face[1::2, ::2] + face[::2, 1::2] + face[1::2, 1::2]
face /= 4.0
height, width = face.shape
# Distort the right half of the image
print('Distorting image...')
distorted = face.copy()
distorted[:, width // 2:] += 0.075 * np.random.randn(height, width // 2)
# Extract all reference patches from the left half of the image
print('Extracting reference patches...')
t0 = time()
patch_size = (7, 7)
data = extract_patches_2d(distorted[:, :width // 2], patch_size)
data = data.reshape(data.shape[0], -1)
data -= np.mean(data, axis=0)
data /= np.std(data, axis=0)
print('done in %.2fs.' % (time() - t0))
# #############################################################################
# Learn the dictionary from reference patches
print('Learning the dictionary...')
t0 = time()
dico = MiniBatchDictionaryLearning(n_components=100, alpha=1, n_iter=500)
V = dico.fit(data).components_
dt = time() - t0
print('done in %.2fs.' % dt)
plt.figure(figsize=(4.2, 4))
for i, comp in enumerate(V[:100]):
plt.subplot(10, 10, i + 1)
plt.imshow(comp.reshape(patch_size), cmap=plt.cm.gray_r,
interpolation='nearest')
plt.xticks(())
plt.yticks(())
plt.suptitle('Dictionary learned from face patches\n' +
'Train time %.1fs on %d patches' % (dt, len(data)),
fontsize=16)
plt.subplots_adjust(0.08, 0.02, 0.92, 0.85, 0.08, 0.23)
# #############################################################################
# Display the distorted image
def show_with_diff(image, reference, title):
"""Helper function to display denoising"""
plt.figure(figsize=(5, 3.3))
plt.subplot(1, 2, 1)
plt.title('Image')
plt.imshow(image, vmin=0, vmax=1, cmap=plt.cm.gray,
interpolation='nearest')
plt.xticks(())
plt.yticks(())
plt.subplot(1, 2, 2)
difference = image - reference
plt.title('Difference (norm: %.2f)' % np.sqrt(np.sum(difference ** 2)))
plt.imshow(difference, vmin=-0.5, vmax=0.5, cmap=plt.cm.PuOr,
interpolation='nearest')
plt.xticks(())
plt.yticks(())
plt.suptitle(title, size=16)
plt.subplots_adjust(0.02, 0.02, 0.98, 0.79, 0.02, 0.2)
show_with_diff(distorted, face, 'Distorted image')
# #############################################################################
# Extract noisy patches and reconstruct them using the dictionary
print('Extracting noisy patches... ')
t0 = time()
data = extract_patches_2d(distorted[:, width // 2:], patch_size)
data = data.reshape(data.shape[0], -1)
intercept = np.mean(data, axis=0)
data -= intercept
print('done in %.2fs.' % (time() - t0))
transform_algorithms = [
('Orthogonal Matching Pursuit\n1 atom', 'omp',
{'transform_n_nonzero_coefs': 1}),
('Orthogonal Matching Pursuit\n2 atoms', 'omp',
{'transform_n_nonzero_coefs': 2}),
('Least-angle regression\n5 atoms', 'lars',
{'transform_n_nonzero_coefs': 5}),
('Thresholding\n alpha=0.1', 'threshold', {'transform_alpha': .1})]
reconstructions = {}
for title, transform_algorithm, kwargs in transform_algorithms:
print(title + '...')
reconstructions[title] = face.copy()
t0 = time()
dico.set_params(transform_algorithm=transform_algorithm, **kwargs)
code = dico.transform(data)
patches = np.dot(code, V)
patches += intercept
patches = patches.reshape(len(data), *patch_size)
if transform_algorithm == 'threshold':
patches -= patches.min()
patches /= patches.max()
reconstructions[title][:, width // 2:] = reconstruct_from_patches_2d(
patches, (height, width // 2))
dt = time() - t0
print('done in %.2fs.' % dt)
show_with_diff(reconstructions[title], face,
title + ' (time: %.1fs)' % dt)
V.shape