Учебник: Линейная алгебра на 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 работает, первый индекс в третьей оси — это данные пикселей красного цвета для нашего изображения. Мы можем получить доступ к ним, используя синтаксис

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

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

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

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

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