Элементы работы с матрицами. Меры расстояний и сходства

Шестаков А.В. Майнор по анализу данных - 26/01/2016

Работа с матрицами

Начнем работу с матрицами пока безотносительно анализа данных. Как я заявил на самом первом семинаре, класс numpy.matrix достаточно глючный. В своё время мне порекомендовали использовать вместо него простые numpy.array.

Пока я буду следовать этому совету, но иногда рядом писать код, если бы что-то было проделано с numpy.matrix

Итак, обзор основных функций:

Арифметические операции, свойства и трансформации

In [ ]:
import numpy as np
import scipy as sp
In [ ]:
# Задаём свою матрицу
A = np.array([[3, 1, 4], [1, 5, 9], [2, 6, 5]])
# A = np.mat([[3, 1, 4], [1, 5, 9], [2, 6, 5]])

print A
In [ ]:
# Транспонируем
A.T
In [ ]:
# A.A # Представляем numpy.matrix в виде numpy.array

# Устраняем одну размерность
A.flatten()
# A.A1
In [ ]:
A + np.eye(A.shape[0])
In [ ]:
A.dot(np.linalg.inv(A))
# A.I * A
In [ ]:
np.linalg.matrix_power(A, 3)
# A**3 # A*A*A
In [ ]:
b = np.array([1,2,3])
b
In [ ]:
A.dot(b)
# A*b
In [ ]:
b.dot(A)
# b.T * A
In [ ]:
A.reshape(9,1)
# A.resize(9,1)
In [ ]:
np.vstack((A,b))
In [ ]:
# np.hstack((A,b))
np.hstack((A,b[:,np.newaxis]))
In [ ]:
e, u = np.linalg.eig(A)
print A.dot(u)
print u.dot(np.diag(e))

Решение линейных систем

Если мы имеет линейную систему $$ Ax = b, $$ то в общем виде $x$ можно попытаться найти так: $$x = A^{-1}b$$

In [ ]:
print A
print b
In [ ]:
x = np.linalg.inv(A).dot(b)
x
In [ ]:
# Проверка:
A.dot(x) - b

В случае матриц не самого приятного вида ине лучших размеров, лучше воспользоваться другими способами, например:

In [ ]:
np.linalg.solve(A, b)
In [ ]:
from scipy import linalg as lin

LU = lin.lu_factor(A)
lin.lu_solve(LU, b)

Расстояние (сходство) между наблюдениями

Для простоты рассмотрим эту тему на простом датасете - Iris.

In [ ]:
from sklearn import datasets

iris = datasets.load_iris()
In [ ]:
iris.keys()
In [ ]:
iris['feature_names']
In [ ]:
X = iris.data.copy()
X[:5, :]
In [ ]:
y = iris.target
print y[:5]
print np.unique(y)

Расчет расстояний (схожести) между наблюдениями призват ответить на вопрос: "насколько близки (схожи) эти наблюдения". Рассмотрим наиболее известные расстояния:

  • Euclidean Distance $$ d(a, b) = \sqrt{\sum_i (a_i - b_i)^2} $$
  • Manhattan Distance $$ d(a, b) = \sum_i |a_i - b_i| $$
  • Cosine Distance $$ d(a, b) = 1 - \frac{a^\top b}{||a||_2\cdot||b||_2}$$
  • Correlation Distance $$ d(a, b) = 1 - \frac{(a-\bar{a})^\top (b - \bar{b})}{||(a-\bar{a})||_2\cdot||(b - \bar{b})||_2}$$
  • ...

Многие из вариантов можно посмотреть здесь

In [ ]:
from scipy.spatial import distance
import matplotlib.pyplot as plt

plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = (10.0, 8.0)

%matplotlib inline
In [ ]:
d = distance.pdist(X, metric='euclidean')
In [ ]:
D = distance.squareform(d)
In [ ]:
D.shape
In [ ]:
plt.figure(figsize=(10,10))
plt.imshow(D, interpolation=None)
plt.colorbar()

Подводные камни:

Работа с разными шкалами признаков

  • Трансформация признаков (dummy-переменные)
  • Как вариант, можно придумать собственную меру, которая в зависимости от шкалы сответствующего признака будет считать то расстояние, которое вам нужно.
  • Взвешивание признаков.
In [ ]:
# Представим, что класс растений - тоже является признаком, причем номинальным.
# Считать любые из вышеперечисленных расстояний не имеет смысла
#
# Но можно сделать так:
dummy = (y[:, None] == np.unique(y)).astype(float)
In [ ]:
dummy[:5,:]

Разные физические размерности признаков

Представьте, что в Iris признаки 2-4 стали измеряться в метрах

In [ ]:
X[:, 1:] = X[:, 1:]/100
X[:5, :]
In [ ]:
d = distance.pdist(X, metric='euclidean',)
D = distance.squareform(d)
In [ ]:
plt.figure(figsize=(10,10))
plt.imshow(D, interpolation=None)
plt.colorbar()

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

Z-score

In [ ]:
X = (X - X.mean(axis=0))/X.std(axis=0)
In [ ]:
print X.mean(axis=0)
print X.std(axis=0)
print X.min()
print X.max()

Min-max

In [ ]:
X = (X - X.min(axis=0))/(X.max(axis=0) - X.min(axis=0))
In [ ]:
print X.mean(axis=0)
print X.std(axis=0)
print X.min()
print X.max()

Проклятье размерности

Надо понижать размерность признакового пространства. Но как?

Метод главных компонент

Метод главных компонент - метод уменьшения размерности данных. Он позволяет отобрать самые "изменчивые" переменные. Но это не просто отбор.

Метод главных компонент производит замену исходных переменных на новые, т.ч.:

1. Новые переменные равны линейной комбинацией старых переменных
2. Новые переменные не коррелируют между собой
3. Новые переменные подбираются так, чтобы максимировать собственную выборочную дисперсию

Для любой матрицы $A$ размера $n \times m$ можно найти разложение вида: $$ A = U \Sigma V^\top ,$$ где

  • $U$ - унитарная матрица, состоящая из собственных векторов $AA^\top$
  • $V$ - унитарная матрица, состоящая из собственных векторов $A^\top A$
  • $\Sigma$ - Диагональная матрица с сингулярными числами

Матрицы $U$ и $V$ ортогональны и могут быть использованы для перехода к ортогональному базису: $$ AV = U\Sigma $$

Сокращение размерности заключается в том, что вместо того, чтобы умножать $A$ на всю матрицу $V$, а лишь на первые $k<m$ её столбцов - матрицу $V'$

Квадраты сингулярных чисел в $\Sigma$ содержат долю дисперсии, которая содержится в соответствующих главных компонентах

In [ ]:
# Вновь возьмем нажу матрицу из Iris
X = iris.data.copy()

# Стандартизируем её
X = (X - X.mean(axis=0))/X.std(axis=0)
In [ ]:
u.shape
In [ ]:
u,s,v = np.linalg.svd(X)
In [ ]:
PC = X.dot(v.T)
In [ ]:
v_expl = s*s/(s*s).sum()
In [ ]:
# Доля дисперсии в главных компонентах
plt.bar(xrange(4), v_expl)
plt.ylabel('\% of variance')
In [ ]:
palette = plt.cm.flag
plt.scatter(PC[:,0], PC[:, 1])

А теперь, если мы успеем, давайте посмотрим на эти данные