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

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

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

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

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

Этот учебник предназначен для людей, которые имеют базовое понимание линейной алгебры и массивов в NumPy и хотят понять, как представлены и могут быть обработаны n-мерные () массивы. В частности, если вы не знаете, как применять общие функции к 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]

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

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

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

>>> from numpy import linalg

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

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

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

где - массив, представляющий изображение в оттенках серого, а и - массивы красного, зеленого и синего каналов, которые мы имели первоначально. Обратите внимание, что мы можем использовать оператор @ (оператор матричного умножения для массивов 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, в данном случае, намного экономичнее в практике, чем построение диагональной матрицы с теми же данными. Для восстановления исходной матрицы мы можем перестроить диагональную матрицу с элементами s на диагонали и с соответствующими размерами для умножения: в нашем случае, должен быть 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)

На графике видно, что хотя у нас 768 сингулярных значений в s, большинство из них (после 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, мы будем использовать функцию 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 значений), мы можем восстановить многие отличительные особенности этого изображения.

Заключение

Конечно, это не лучший метод приближения изображения. Однако на самом деле существует результат в линейной алгебре, который гласит, что приближение, построенное нами выше, является наилучшим приближением к исходной матрице с точки зрения нормы разности. Для получения дополнительной информации см. G. H. Golub and C. F. Van Loan, Matrix Computations, Baltimore, MD, Johns Hopkins University Press, 1985.

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

• Учебник по Python
• Справочник по NumPy
• Учебник по SciPy
• Заметки по лекциям по SciPy
• Словарь MATLAB, R, IDL, NumPy/SciPy