Сегодня я хочу показать вам возможности анализа главных компонентов (АГК). Это метод уменьшения размерности данных, повышения интерпретируемости, но в то же время минимизации потери информации.

При этом давайте посмотрим, как происходит это волшебство! Я продемонстрирую код Python для реализации PCA с нуля. Это поможет вам понять концепцию более подробно, а также облегчит практическое обучение. Чтобы визуально увидеть силу PCA и понять ее, я выбрал данные изображения. Независимо от того, данные состоят из различных видов аниме-вайфу! в наборе данных все еще есть некоторые статистические закономерности, которые PCA может изучать и интерпретировать изображения с использованием меньшего количества функций.

Итак, начнем! Я создаю код, который я разработал в качестве ассистента преподавателя в своей магистратуре с помощью и под руководством проф. Вишванат Гопалакришнан, чтобы учащиеся поняли и внедрили код для PCA по предмету Математика для машинного обучения. Я использую подмножество набора данных, доступного на Kaggle. Набор данных можно найти здесь: Набор данных. Я предполагаю, что вы используете Google Colab. Однако это также можно сделать в мощной системе в автономном режиме в Jupyter Notebook, используя при необходимости сюжет в автономном режиме. Я бы также рекомендовал использовать темный режим для лучшего опыта.

Первым шагом будет подключение вашего диска к Google Colab. Это легко сделать с помощью следующего кода:

#code to mount drive
from google.colab import drive
drive.mount('/content/drive')

Просто следуйте инструкциям, предоставьте разрешение, и вы будете готовы к следующему шагу — загрузке набора данных. Обратите внимание, что я создал папку с именем «PCA» в MyDrive и загрузил набор данных в эту папку.

Теперь мы импортируем необходимые библиотеки. Просто используйте приведенный ниже фрагмент. Здесь мы также устанавливаем значение рендерера по умолчанию для Google Colab, чтобы наши графики были видны в блокноте.

#importing libraries
import os
import sys 
import cv2 as cv
import numpy as np
import plotly.io as pio
import plotly.graph_objs as go


from PIL import Image
from skimage import color
from plotly import subplots
from sklearn.model_selection import train_test_split


#setting the rederer as colab
pio.renderers.default = "colab"

Теперь мы загрузим набор данных. Но прежде чем я продолжу, я хотел бы сделать одно предостережение. Я выбрал значения количества изображений до 500 и размеры изображения 64x64, чтобы вы могли безупречно запускать код в пределах предоставленного ограничения ресурсов в Google Colab. Однако, если вы запускаете его локально на мощной системе с ОЗУ 32 ГБ или выше, вы можете изменить эти значения, чтобы получить лучшие результаты. Например, количество изображений, предоставленных в наборе данных, равно 1000. Или вы также можете увеличить размеры изображения до 128x128.

При всем при этом нигде нет кода для загрузки изображений в систему. Я объясню код подробно шаг за шагом.

#Loading data and reducing size to 64 x 64 pixels
IMG_DIR = '/content/drive/MyDrive/PCA/Dataset'
X = []
X_flat = []
count = 1
size = 64
total = 500
print("Loading...")
for img in os.listdir(IMG_DIR):
    if count == total + 1:
        break
    sys.stdout.write("\r" + str(count) + " / " + str(total))
    sys.stdout.flush()
    img_array = cv.imread(os.path.join(IMG_DIR, img), cv.IMREAD_GRAYSCALE)
    img_pil = Image.fromarray(img_array)
    img_64x64 = np.array(img_pil.resize((size, size), Image.ANTIALIAS))
    X.append(img_64x64)
    img_array = img_64x64.flatten()
    X_flat.append(img_array)
    count += 1
print()
print("Done!")
  • Расположение набора данных загружается в переменную «IMG_DIR».
  • Для хранения изображений создаются два списка «X» и «X_flat».
  • Счетчик count инициализируется равным 1, чтобы отслеживать, какое изображение загружается.
  • «size» инициализирован равным 64, чтобы сохранить размеры изображения 64x64.
  • «Всего» инициализируется значением 500, чтобы загрузить только 500 изображений из набора данных.
  • Теперь доступ к каждому изображению в наборе данных осуществляется с помощью цикла for.
  • Условие «если» используется для отслеживания количества уже загруженных изображений. Если «count» становится больше, чем «total», цикл завершается.
  • Постоянная обратная связь о том, какое изображение загружается, предоставляется с помощью двух операторов sys.
  • Каждое изображение одно за другим считывается в переменную img_array с помощью команды openCV imread(). Также обратите внимание, что каждое изображение преобразуется в оттенки серого.
  • Изображение (в настоящее время массив numpy) преобразуется в изображение PIL с помощью команды Image.fromarray(). Это сделано для облегчения изменения размера изображения.
  • Затем размер изображения изменяется до 64x64 с использованием метода Anti Aliasing и снова преобразуется обратно в массив numpy.
  • Затем он добавляется в список «X».
  • Плоская версия изображений также хранится в списке «X_flat», где изображение 64x64 формируется в массив (4096) numpy. Мы будем использовать «X_flat» для дальнейшей обработки. «X» просто для нас, чтобы увидеть, как выглядят изображения сейчас.

Я надеюсь, что приведенные выше пункты сделали код очень понятным. Теперь давайте посмотрим, загрузились изображения или нет, а также получим представление о том, как они выглядят. Ниже приведен фрагмент, показывающий первые 25 изображений, загруженных в «X».

#visualizing some images
size = 5
count = 0
fig = subplots.make_subplots(rows = size, cols = size, 
                 vertical_spacing = 0.06, horizontal_spacing = 0.02)
for row in range(size):
  for col in range(size):
    fig.add_trace(go.Image(z = color.gray2rgb(X[count])), 
                  row = row + 1, col = col + 1)
    count += 1
fig["layout"].update(title = "Anime Images", template = "plotly_dark", height = 900)
fig.show()

Как только вы запустите код, вы увидите 25 изображений очень красивых вайфусов! Пожалуйста, обратите внимание, что каждое из них представляет собой изображение в градациях серого и не отображает изображения в градациях серого. Поэтому мне пришлось использовать утилиту gray2rgb() в skimage, чтобы изображения в градациях серого имели 3 канала, как в стандартном изображении RGB. Цветовая модель по умолчанию в go.Image() — только RGB. Если вы не понимаете этого построения трасс в plotly, я настоятельно рекомендую вам заглянуть в другую мою статью о plotly. Это поможет вам пройти весь путь от новичка до продвинутого графолога в кратчайшие сроки! Доступ к статье можно получить здесь: Полное введение в сюжет, от новичка до продвинутого. Вот ожидаемый результат.

Хорошо! давайте не будем отвлекаться на эти симпатичные лица, и быстро перейдем к основной повестке дня, которая у нас есть сегодня. Теперь я конвертирую «X_flat» в массив numpy в приведенном ниже фрагменте. Вы увидите вывод (500, 4096). Пожалуйста, помните, что одно сплющенное изображение было (4096), и я объединил все 500 изображений, которые мы загрузили вместе. Размеры набора данных могут ввести вас в заблуждение, если вы поверите, что существует 4096 функций, но, пожалуйста, обратите внимание, что каждое изображение представляет собой функцию. Следовательно, существует 500 функций (которое будет уменьшено до 250 из-за разделения тестового поезда дальше по линии). ). [Вопрос 1: Можете ли вы сказать, что произойдет, если мы возьмем 4096 в качестве признаков после того, как вы проработаете всю концепцию?]

#converting X_flat to numpy array
X_flat = np.asarray(X_flat)
X_flat.shape

Вывод: (500, 4096)

Теперь давайте разделим тест-поезд. Обратите внимание, что разделение должно быть 50–50. Вы можете видеть, что для атрибута test_size установлено значение 0,5. Кроме того, я не использовал транспонирование, чтобы данные больше походили на естественное понимание того, что столбцы являются функциями, а строки - точками данных. Используйте атрибут «random_state», чтобы получить такое же разделение, как у меня. [Вопрос 2: Можете ли вы ответить, почему нам нужно разделить вещи на равные части, после того как вы прошли через всю концепцию?]

#Test-Train split
X_train, X_test = train_test_split(X_flat, 
                           test_size = 0.5, random_state = 69)
X_train = X_train.T
X_test = X_test.T
print(X_train.shape)
print(X_test.shape)

Вывод: (4096, 250) ‘\n’ (4096, 250)

Итак, теперь у нас есть данные, готовые для обучения PCA. Итак, приступим к реализации PCA! Но сначала позвольте мне дать вам небольшую информацию о том, как работает PCA. PCA можно разбить на 5 шагов. Я пройдусь по ним один за другим последовательно. Если вы новичок в этой теме, это может быть очень сложно, но я прошу вас набраться терпения и прочитать еще раз, если потребуется. Я выбрал очень точный язык, чтобы очень тонко изложить идею. К сожалению, впереди много чтения!

  1. Стандартизация: мы стандартизируем данные, чтобы каждая функция имела равный вклад в анализ. Представьте себе линию регрессии, построенную для данных, в которых значения по оси x отличаются небольшими величинами, например, от 0,3 до 0,9 (небольшое изменение данных), а значения по оси Y различаются на суммы, исчисляемые тысячами (большое отклонение данных). Такая линия будет сильно скошена по оси y и не сможет следовать дисперсии данных по оси x. Таким образом, нам нужно привести все признаки к сопоставимому масштабу, чтобы ни один из них не был доминирующим. При стандартизации мы делаем значения признаков равными 0, а стандартное отклонение равным 1.
  2. Ковариация: мы создаем ковариационную матрицу данных, чтобы выяснить, есть ли какая-либо связь между функциями. Потому что функции могут быть сильно коррелированы, что предполагает избыточную информацию. Таким образом, мы пытаемся удалить такие функции. Итак, чтобы получить представление о сильно коррелирующих функциях, мы строим ковариационную матрицу.
  3. Собственные значения и собственные векторы: этот шаг помогает нам определить основные компоненты данных. Основные компоненты — это новые значения признаков, построенные как линейные комбинации или смеси исходных значений признаков. Это делается таким образом, чтобы новые значения признаков были некоррелированными, и большая часть информации была доступна в первом компоненте. Затем оставшиеся в следующем и так далее. Таким образом, мы можем удалить менее важные основные компоненты, чтобы уменьшить размерность данных. Геометрически собственные векторы (векторы характеристик) — это те векторы, которые не меняют своего направления при применении к линейному преобразованию, а собственные значения представляют величину, на которую они масштабируются. Таким образом, как только мы находим собственные векторы ковариационной матрицы, мы получаем направления данных, которые объясняют максимальную величину дисперсии в зависимости от собственных значений (большое собственное значение означает большую дисперсию в направлении этого собственного вектора). А поскольку ковариационная матрица симметрична, собственные векторы будут перпендикулярны. Как удобно! так что теперь мы можем использовать их в качестве осей для объяснения набора данных. Таким образом, мы сортируем собственные векторы на основе их собственных значений в порядке убывания, чтобы сначала получить важные главные компоненты.
  4. Векторы признаков: это те векторы, которые делают разрез. То есть выбранные основные компоненты, которые, по нашему мнению, достаточны для представления данных без потери большого количества информации.
  5. Проекция на главные компоненты: до сих пор мы только что нашли направления, которые лучше всего объясняют данные. Теперь мы переориентируем данные с исходных осей на те, которые описываются главными компонентами. Это даст нам данные с уменьшенными размерами.

Фу! это кажется много работы! Но если вы правильно понимаете математику и используете концепцию векторизации Python, кодирование всего этого будет легкой задачей. Я бы порекомендовал подробно прочитать эти концепции, если вы не знакомы с некоторыми ресурсами Гилберта Стрэнга. Так что давайте не будем больше ждать и пачкать руки!

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

Во-первых, это функция стандартизации. (Шаг 1)

#function for standardizing image
def Standardize(X):
    #calcualte the mean of each column mu   
    mu = np.mean(X, axis = 0) 
    
    #Substract the mean from X
    X = X - mu  
    
    #calcualte the standard deviation of each column std
    std = np.std(X, axis = 0)  
    
    #handleing zero standard deviation
    std_filled = std.copy()
    std_filled[std == 0] = 1.0
    
    #calculate standardized X called Xbar
    Xbar = (X-mu) / std_filled  
    
    return Xbar, mu, std

Затем функция для вычисления собственных значений и собственных векторов и сортировки их на основе собственных значений в порядке убывания. (Шаг 3)

#function for calcualting eigen values and eigen vectors
def eig(S):
    #calculate the eigen values and eigen vectors
    eig_val, eig_vec = np.linalg.eigh(S)  
    
    #sorting them in decrasing order
    sorted_eig  = np.argsort(-eig_val)
    eig_val = eig_val[sorted_eig]
    eig_vec = eig_vec[:, sorted_eig]
    
    return (eig_val, eig_vec)

Затем идет функция для проецирования данных. (Шаг 5). Я знаю, что перескакиваю через ступеньки, но, как я уже сказал, это были коды, которые я ожидал, что студенты напишут, поэтому я собрал их здесь как функции, объединенные вместе. Шаг 2 и шаг 4 выполняются в правильном порядке во фрагменте, следующем за этим фрагментом.

#function for projection matrix
def projection_matrix(B):
    #calculate the projection matrix P
    P = B @ B.T 
    
    return P

Теперь основную функцию выполнять ППШ.

#implementing PCA
def PCA(X, num_components):
    #calculate the data covariance matrix S
    S = np.cov(X)  
    
    #now find eigenvalues and corresponding eigenvectors for S 
    #by implementing eig().
    eig_vals, eig_vecs = eig(S)
    
    #select eigen vectors
    U = eig_vecs[:, range(num_components)]
    
    #reconstruct the images from the lowerdimensional representation
    #to do this, we first need to find the projection_matrix 
    #(which you implemented earlier)
    #which projects our input data onto the vector space spanned 
    #by the eigenvectors
    P = projection_matrix(U) # projection matrix
    
    return P

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

#standardizind
Xbar_train, mu_train, std_train = Standardize(X_train)
Xbar_test, mu_test, std_test = Standardize(X_test)

Теперь я пишу функцию для вычисления среднеквадратичной ошибки (MSE) для обучения.

#function for mean square error
def mse(predict, actual):
    return np.square(predict - actual).sum(axis = 1).mean()

Теперь, используя обучающие данные, я строю проекционную матрицу, которая проецирует заданные данные в пространство признаков собственного вектора, которое было рассчитано с использованием показанной здесь линии (проекция = PCA (Xbar_train.T, num_component). Затем я использую эту проекционную матрицу. сделаны из обученных данных для проецирования тестовых данных в это пространство собственных векторов, чтобы увидеть, что происходит, используя линию (reconst = Xbar_test @ проекция). Можно ли также реконструировать невидимые данные (неизвестная аниме-вайфу) с помощью этой проекционной матрицы? Давайте выясним!

#calculating loss and reconstructing images
loss = []
reconstructions = []
max_components = len(X_train.T)
print("Processing...")
animation = np.arange(1, max_components + 1, 1)
for num_component in range(1, max_components + 1):
    sys.stdout.write("\r" + str(animation[num_component - 1]) + \
                      " / " + str(max_components))
    sys.stdout.flush()
    projection = PCA(Xbar_train.T, num_component)
    reconst = Xbar_test @ projection
    error = mse(reconst, Xbar_test)
    reconstructions.append(reconst)
    loss.append((num_component, error))
print()
print("Done!")

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

#visualizing mse vs number of principal components
trace = go.Scatter(x = loss[:, 0], y = loss[:, 1])
data = [trace]
fig = go.Figure(data)
fig.update_layout(title = "MSE vs number of principal components", 
                  xaxis_title = "Number of principal components", 
                  yaxis_title = "MSE", template = "plotly_dark")
fig.show()

Вы можете ясно видеть уменьшение MSE по мере того, как я увеличиваю количество основных компонентов, используемых для создания проекционной матрицы. Это означает, что качество восстановленных изображений становится все лучше и лучше по мере того, как мы включаем все больше и больше основных компонентов. Почти до нуля! если используются все 250 главных компонент.

Давайте узнаем, как выглядят наши вайфу! Я не могу дождаться!

Во-первых, помните, что данные должны быть нестандартизированы. В противном случае изображения не будут иметь правильных значений пикселей. Мы можем сделать это, используя следующий простой код.

#unstandardizing the reconstructed images
reconstructions = np.asarray(reconstructions)
reconstructions = reconstructions * std_test + mu_test 
loss = np.asarray(loss)
print("Done")

После этого мы можем взглянуть на наш вайфус, используя следующий код.

#plotting the reconstructed images
col = 5
row = 2
size = 10
skip = 20
images = 3


component = 0
titles = []
for loop in range(size):
  titles.append(str(component) + " components")
  component += skip


figs = [0]*images
for num in range(images):
  figs[num] = subplots.make_subplots(rows = row, 
                                cols = col, subplot_titles = titles)
  component = 0
  for r in range(row):
    for c in range(col):
      figs[num].add_trace(go.Image(
             z = color.gray2rgb(
             reconstructions[component, : , num].reshape(64, 64))), 
             row = r + 1, col = c + 1)
      component += skip
  figs[num]["layout"].update(title = "Recontructed Images", 
                             template = "plotly_dark")
for num in range(images):
  figs[num].show()

Вывод должен выглядеть примерно так:

3 вайфуса готовы! Вы можете видеть, что я могу получить несколько очаровательных лиц, просто используя 180 основных компонентов. Это означает, что 70 функций были в основном бесполезными (250–180). Представьте себе, что вы обнаружите, что около 30% функций, которые вы предоставляете своей модели, являются избыточными! Эта информация может помочь повысить производительность модели в целом.

Мы все любим аниме, и я надеюсь, что теперь вы каждый раз будете ассоциировать PCA с ним! Я надеюсь, что я сделал эту тему легкой и интересной для вас, чтобы понять. PCA — очень сложный, но очень важный алгоритм, который чаще помогает, чем нет. Я бы порекомендовал продолжить изучение линейной алгебры, если вы не понимаете, что происходит должным образом. Используйте ресурсы Гилберта Стрэнга, чтобы получить конкретное представление о том, что такое линейные преобразования, собственные значения и собственные векторы, их значение, проекции и т. д.

Это все от меня на этот раз. Пожалуйста, извините меня сейчас. Ответьте на два вопроса в комментариях.

Я надеюсь, что вы узнали что-то из этой статьи, и я желаю вам удачи.

Пока!

Особая благодарность проф. Вишванату Гопалакришнану за помощь в устранении основных ошибок в коде.