Учебник: Линейная алгебра на n-мерных массивах

Предварительные знания

Прежде чем читать этот учебник, вы должны немного знать Python. Если вы хотите освежить свои знания, взгляните на учебник по Python.

Если вы хотите выполнить примеры в этом учебнике, на вашем компьютере также должны быть установлены matplotlib и SciPy.

Профиль обучающегося

Этот учебник предназначен для людей, которые имеют базовое понимание линейной алгебры и массивов в NumPy и хотят понять, как представлены и могут быть обработаны n-мерные (\(n>=2\)) массивы. В частности, если вы не знаете, как применять общие функции к n-мерным массивам (без использования циклов for), или если вы хотите понять свойства осей и формы для n-мерных массивов, этот учебник может оказаться полезным.

Цели обучения

После этого учебника вы сможете:

• Понять разницу между одномерными, двумерными и n-мерными массивами в NumPy;
• Понять, как применять некоторые операции линейной алгебры к n-мерным массивам без использования циклов for;
• Понять свойства осей и формы для n-мерных массивов.

Содержание

В этом учебнике мы будем использовать разложение матрицы из линейной алгебры, разложение по сингулярным значениям, для генерации сжатого приближения изображения. Будем использовать изображение face из модуля scipy.misc:

>>> from scipy import misc
>>> img = misc.face()

Теперь img — это массив NumPy, как мы можем видеть, используя функцию type:

>>> type(img)
<class 'numpy.ndarray'>

Мы можем отобразить изображение с помощью функции matplotlib.pyplot.imshow:

>>> import matplotlib.pyplot as plt
>>> plt.imshow(img)

Свойства формы, оси и массива

Обратите внимание, что в линейной алгебре размерность вектора относится к количеству элементов в массиве. В NumPy он вместо этого определяет количество осей. Например, одномерный массив — это вектор, такой как [1, 2, 3], двумерный массив — матрица и так далее.

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

>>> img.shape
(768, 1024, 3)

Вывод — это кортеж из трех элементов, что означает, что это трехмерный массив. На самом деле, поскольку это цветное изображение, и мы использовали функцию imread для его чтения, данные организованы в три двумерных массива, представляющих цветовые каналы (в данном случае красный, зеленый и синий — RGB). Вы можете увидеть это, посмотрев на форму выше: она указывает, что у нас есть массив из 3 матриц, каждая с формой 768x1024.

Кроме того, используя свойство ndim этого массива, мы можем видеть, что

>>> img.ndim
3

NumPy называет каждую размерность axis. Из-за того, как работает imread, первый индекс в 3-й оси — это данные красного пикселя для нашего изображения. Мы можем получить к ним доступ, используя синтаксис

>>> img[:, :, 0]
array([[121, 138, 153, ..., 119, 131, 139],
       [ 89, 110, 130, ..., 118, 134, 146],
       [ 73,  94, 115, ..., 117, 133, 144],
       ...,
       [ 87,  94, 107, ..., 120, 119, 119],
       [ 85,  95, 112, ..., 121, 120, 120],
       [ 85,  97, 111, ..., 120, 119, 118]], dtype=uint8)

Из вывода выше мы видим, что каждое значение в img[:,:,0] является целым числом от 0 до 255, представляющим уровень красного цвета в каждом соответствующем пикселе изображения (имейте в виду, что это может быть другим, если вы используете собственное изображение вместо scipy.misc.face).

Как ожидалось, это матрица 768x1024:

>>> img[:, :, 0].shape
(768, 1024)

Поскольку мы собираемся выполнять операции линейной алгебры с этими данными, может быть интереснее иметь вещественные числа от 0 до 1 в каждом элементе матриц для представления значений RGB. Мы можем сделать это, задав

>>> img_array = img / 255

Эта операция, деление массива на скаляр, работает благодаря правилам распространения NumPy. (Обратите внимание, что в реальных приложениях было бы лучше использовать, например, утилиту img_as_float из scikit-image).

Вы можете проверить, что вышеуказанное работает, выполнив некоторые тесты; например, запросив максимальные и минимальные значения для этого массива:

>>> img_array.max(), img_array.min()
(1.0, 0.0)

или проверив тип данных в массиве:

>>> img_array.dtype
dtype('float64')

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

>>> red_array = img_array[:, :, 0]
>>> green_array = img_array[:, :, 1]
>>> blue_array = img_array[:, :, 2]

Операции по оси

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

Чтобы продолжить, импортируйте подмодуль линейной алгебры из NumPy:

>>> from numpy import linalg

Для извлечения информации из заданной матрицы мы можем использовать SVD для получения 3 массивов, которые можно перемножить, чтобы получить исходную матрицу. Согласно теории линейной алгебры, для данной матрицы \(A\), можно вычислить следующее произведение:

где \(U\) и \(V^T\) квадратные, а \(\Sigma\) имеет тот же размер, что и \(A\). \(\Sigma\) — диагональная матрица и содержит сингулярные значения \(A\), упорядоченные от наибольшего к наименьшему. Эти значения всегда неотрицательны и могут использоваться как показатель «важности» некоторых признаков, представленных матрицей \(A\).

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

где \(Y\) — массив, представляющий изображение в оттенках серого, а \(R, G\) и \(B\) — массивы каналов красного, зеленого и синего, которые у нас были изначально. Заметьте, что мы можем использовать оператор @ (оператор умножения матриц для массивов NumPy, см. numpy.matmul) для этого:

>>> img_gray = img_array @ [0.2126, 0.7152, 0.0722]

Теперь img_gray имеет форму

>>> img_gray.shape
(768, 1024)

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

В нашем случае мы приближаем серо-шкальную часть изображения, поэтому будем использовать цветовую карту gray:

>>> plt.imshow(img_gray, cmap="gray")

Теперь, применяя функцию linalg.svd к этой матрице, мы получаем следующее разложение:

>>> U, s, Vt = linalg.svd(img_gray)

Давайте проверим, что это ожидалось:

>>> U.shape, s.shape, Vt.shape
((768, 768), (768,), (1024, 1024))

Обратите внимание, что s имеет особую форму: у него только одна размерность. Это означает, что некоторые функции линейной алгебры, которые ожидают двумерные массивы, могут не работать. Например, из теории можно было бы ожидать, что s и Vt совместимы для умножения. Однако это не так, так как у s нет второй оси. Выполнение

>>> s @ Vt
Traceback (most recent call last):
  ...
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0,
with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 1024 is different from
768)

приводит к ValueError. Это происходит потому, что использование одномерного массива для s, в этом случае, гораздо экономичнее на практике, чем построение диагональной матрицы с теми же данными. Чтобы восстановить исходную матрицу, мы можем перестроить диагональную матрицу \(\Sigma\) с элементами s по диагонали и с соответствующими размерами для умножения: в нашем случае, \(\Sigma\) должна быть 768x1024, так как U является 768x768, а Vt является 1024x1024.

>>> import numpy as np
>>> Sigma = np.zeros((768, 1024))
>>> for i in range(768):
...     Sigma[i, i] = s[i]

Теперь мы хотим проверить, является ли восстановленная U @ Sigma @ Vt близкой к исходной img_gray матрице.

Приближение

Модуль linalg включает функцию norm, которая вычисляет норму вектора или матрицы, представленных в массиве NumPy. Например, из объяснения SVD выше, мы ожидали бы, что норма разницы между img_gray и восстановленным произведением SVD будет небольшой. Как ожидалось, вы должны увидеть что-то вроде

>>> linalg.norm(img_gray - U @ Sigma @ Vt)
1.3926466851808837e-12

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

Мы также могли использовать функцию numpy.allclose, чтобы убедиться, что восстановленное произведение действительно близко к нашей исходной матрице (разность между двумя массивами мала):

>>> np.allclose(img_gray, U @ Sigma @ Vt)
True

Чтобы увидеть, является ли приближение приемлемым, мы можем проверить значения в s:

>>> plt.plot(s)

На графике мы видим, что хотя в s у нас 768 сингулярных значений, большинство из них (после 150-го элемента и далее) довольно маленькие. Поэтому может быть целесообразно использовать только информацию, относящуюся к первым (скажем, 50) сингулярным значениям, чтобы построить более экономичное приближение к нашему изображению.

Идея состоит в том, чтобы считать все, кроме первых k сингулярных значений в Sigma (которые такие же, как в s ) нулями, сохранив U и Vt нетронутыми, и вычислить произведение этих матриц как приближение.

Например, если мы выберем

>>> k = 10

мы можем построить приближение, выполнив

>>> approx = U @ Sigma[:, :k] @ Vt[:k, :]

Обратите внимание, что мы должны были использовать только первые k строки Vt, так как все остальные строки будут умножены на нули, соответствующие сингулярным значениям, которые мы исключили из этого приближения.

>>> plt.imshow(approx, cmap="gray")

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

Применение ко всем цветам

Теперь мы хотим выполнить ту же операцию, но для всех трех цветов. Наш первый импульс — повторить ту же операцию, что и выше, для каждой матрицы цветов по отдельности. Однако NumPy’s broadcasting позаботится об этом за нас.

Если наш массив имеет более двух измерений, то SVD можно применить ко всем осям сразу. Однако функции линейной алгебры в NumPy ожидают увидеть массив вида (N, M, M), где первая ось представляет количество матриц.

В нашем случае,

>>> img_array.shape
(768, 1024, 3)

поэтому нам нужно переставить ось в этом массиве, чтобы получить форму, подобную (3, 768, 1024). К счастью, функция numpy.transpose может это сделать за нас:

np.transpose(x, axes=(i, j, k))

указывает, что ось будет переупорядочена таким образом, что окончательная форма транспонированного массива будет переупорядочена в соответствии с индексами (i, j, k).

Давайте посмотрим, как это работает для нашего массива:

>>> img_array_transposed = np.transpose(img_array, (2, 0, 1))
>>> img_array_transposed.shape
(3, 768, 1024)

Теперь мы готовы применить SVD:

>>> U, s, Vt = linalg.svd(img_array_transposed)

Наконец, чтобы получить полное приближённое изображение, нам нужно собрать эти матрицы в приближение. Теперь обратите внимание, что

>>> U.shape, s.shape, Vt.shape
((3, 768, 768), (3, 768), (3, 1024, 1024))

Чтобы построить окончательную матрицу приближения, мы должны понять, как работает умножение по различным осям.

Произведения с n-мерными массивами

Если вы работали ранее только с одномерными или двумерными массивами в NumPy, вы можете использовать numpy.dot и numpy.matmul (или оператор @ ) взаимозаменяемо. Однако для n-мерных массивов они работают по-разному. Для получения более подробной информации см. документацию numpy.matmul.

Теперь, чтобы построить наше приближение, нам сначала нужно убедиться, что наши сингулярные значения готовы для умножения, поэтому мы создаём нашу Sigma матрицу аналогично тому, что мы делали раньше. Массив Sigma должен иметь размер (3, 768, 1024). Чтобы добавить сингулярные значения на диагональ Sigma, мы будем использовать функцию numpy.fill_diagonal из NumPy, используя каждую из 3 строк в s в качестве диагонали для каждой из 3 матриц в Sigma:

>>> Sigma = np.zeros((3, 768, 1024))
>>> for j in range(3):
...     np.fill_diagonal(Sigma[j, :, :], s[j, :])

Теперь, если мы хотим перестроить полное SVD (без приближения), мы можем сделать

>>> reconstructed = U @ Sigma @ Vt

Обратите внимание, что

>>> reconstructed.shape
(3, 768, 1024)

и

>>> plt.imshow(np.transpose(reconstructed, (1, 2, 0)))

должно дать вам изображение, неотличимое от исходного (хотя мы можем ввести ошибки с плавающей запятой для этого восстановления). В действительности, вы можете увидеть сообщение об ошибке, которое говорит “Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).” Это ожидается от манипуляций, которые мы только что выполнили над исходным изображением.

Теперь, чтобы сделать приближение, мы должны выбрать только первые k сингулярные значения для каждого цветового канала. Это можно сделать с помощью следующего синтаксиса:

>>> approx_img = U @ Sigma[..., :k] @ Vt[..., :k, :]

Вы можете видеть, что мы выбрали только первые k компоненты последней оси для Sigma (это означает, что мы использовали только первые k столбцы каждой из трёх матриц в стеке), и что мы выбрали только первые k компоненты во второй-последней оси Vt (это означает, что мы выбрали только первые k строк каждой матрицы в стеке Vt и все столбцы). Если вы не знакомы с синтаксисом эллипса, это заполнитель для других осей. Для получения более подробной информации см. документацию по Индексированию.

Теперь,

>>> approx_img.shape
(3, 768, 1024)

что не является правильной формой для отображения изображения. Наконец, переупорядочив оси обратно в нашу исходную форму (768, 1024, 3), мы можем увидеть наше приближение:

>>> plt.imshow(np.transpose(approx_img, (1, 2, 0)))

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

Заключение

Конечно, это не лучший метод приближения изображения. Однако, на самом деле, в линейной алгебре есть результат, который говорит о том, что приближение, которое мы построили выше, является наилучшим, которое мы можем получить для исходной матрицы с точки зрения нормы разности. Для получения дополнительной информации см. Г. Х. Голуб и К. Ф. Ван Лоан, Матричные вычисления, Балтимор, МД, Издательство Университета Джона Хопкинса, 1985.

Дополнительные материалы

• Руководство по Python
• Справочник по командам
• Руководство по SciPy
• SciPy Lecture Notes
• Словарь Matlab, R, IDL, NumPy/SciPy