Реализация SVD на Python с рандомизированной линейной алгеброй
Матричная декомпозиция - мощный инструмент для решения многих задач машинного обучения, который широко используется для сжатия данных, уменьшения размерности и обучения разреженности, и это лишь некоторые из них. Во многих случаях для аппроксимации матрицы данных структурой низкого ранга оптимальным выбором часто считается разложение по сингулярным значениям (SVD). Однако точный и эффективный SVD для больших матриц данных (например, матрицы размером 8k на 10k) является сложной вычислительной задачей. Чтобы разрешить SVD в этой ситуации, было разработано множество алгоритмов с применением методов рандомизированной линейной алгебры. Одним из наиболее важных алгоритмов является рандомизированный SVD, который с точки зрения конкуренции эффективен для разложения любой большой матрицы с относительно низким рангом.
Этот пост представит предварительную и важную идею рандомизированной SVD. Чтобы помочь читателям лучше понять рандомизированный SVD, мы также предоставляем соответствующую реализацию Python в этом посте. Кроме того, здесь можно найти блокнот Jupyter этого поста.
Предварительный
Формула SVD
Начнем с напоминания о концепции СВД. Как вы, возможно, уже знаете, SVD - одна из самых важных формул разложения в линейной алгебре. Для любой данной матрицы A SVD имеет вид
A = UΣV^T
где матрицы U и V состоят из левых и правых сингулярных векторов соответственно. Диагональные элементы Σ являются сингулярными значениями.
Пример небольшой матрицы
Возьмем, к примеру, матрицу 3 на 3, мы можем вычислить SVD, используя numpy.linalg.svd() в Python. Давайте посмотрим:
import numpy as np
A = np.array([[1, 3, 2],
[5, 3, 1],
[3, 4, 5]])
u, s, v = np.linalg.svd(A, full_matrices = 0)
print('Left singular vectors:')
print(u)
print()
print('Singular values:')
print(s)
print()
print('Right singular vectors:')
print(v)
print()
В этом случае сингулярные значения равны 9,3427, 3,2450 и 1,0885.
Left singular vectors: [[-0.37421754 0.28475648 -0.88253894] [-0.56470638 -0.82485997 -0.02669705] [-0.7355732 0.48838486 0.46948087]] Singular values: [9.34265841 3.24497827 1.08850813] Right singular vectors: [[-0.57847229 -0.61642675 -0.53421706] [-0.73171177 0.10269066 0.67383419] [ 0.36051032 -0.78068732 0.51045041]]
Рандомизированный СВД
Основная идея
Рандомизированный SVD можно разбить на три основных этапа. Для любой заданной матрицы m -by- n A, если мы наложим целевой ранг k с k ‹min (m, n), тогда первый шаг, как показано на рисунке 2, -
- 1) сгенерировать гауссову матрицу случайных чисел Ω размером n- на -k,
- 2) вычислить новую матрицу m на k Y,
- и 3) применить QR-разложение к матрице Y.
Обратите внимание, что первый шаг должен вернуть матрицу m -by- k Q.
Затем второй шаг, показанный на рисунке 3, - это
- 4) выведите матрицу k -by- n B путем умножения транспонированной матрицы Q и матрица A вместе,
- и 5) вычислить SVD матрицы B. Здесь вместо вычисления SVD исходной матрицы A, B это меньшая матрица для работы.
Поскольку сингулярные значения (т. Е. Σ) и правые сингулярные векторы (т. Е. V) матрицы B также являются сингулярными значениями и правыми сингулярными векторами исходной матрицы A, мы должны сохранить вычисленные сингулярные значения и правые сингулярные векторы. матрицей B на этом этапе.
Как показано на рисунке 3, если мы объединим матрицу Q, полученную на первом этапе, с левыми сингулярными векторами B , мы можем получить левые сингулярные векторы (т. е. U) матрицы A в третьем шаг.
Пример небольшой матрицы
Несмотря на то, что мы узнали основную идею рандомизированного SVD выше, было бы не совсем понятно, если бы не было интуитивно понятного примера. С этой целью мы следуем вышеупомянутой небольшой матрице SVD.
Во-первых, давайте попробуем написать Python-функцию рандомизированного SVD. Здесь мы будем использовать две функции Numpy, то есть np.linalg.qr() и np.linalg.svd().
import numpy as np
def rsvd(A, Omega):
Y = A @ Omega
Q, _ = np.linalg.qr(Y)
B = Q.T @ A
u_tilde, s, v = np.linalg.svd(B, full_matrices = 0)
u = Q @ u_tilde
return u, s, v
Теперь давайте протестируем его с помощью матрицы 3 на 3 (rank = 2 для обозначения k с помощью k ‹min (m, n )):
np.random.seed(1000)
A = np.array([[1, 3, 2],
[5, 3, 1],
[3, 4, 5]])
rank = 2
Omega = np.random.randn(A.shape[1], rank)
u, s, v = rsvd(A, Omega)
print('Left singular vectors:')
print(u)
print()
print('Singular values:')
print(s)
print()
print('Right singular vectors:')
print(v)
print()
Результатом этого рандомизированного примера SVD является:
Left singular vectors: [[ 0.38070859 0.60505354] [ 0.56830191 -0.74963644] [ 0.72944767 0.26824507]] Singular values: [9.34224023 3.02039888] Right singular vectors: [[ 0.57915029 0.61707064 0.53273704] [-0.77420021 0.21163814 0.59650929]]
Напомним, что сингулярные значения этой матрицы - 9,3427, 3,2450 и 1,0885. В этом случае рандомизированный SVD имеет первые два сингулярных значения: 9.3422 и 3.0204.
Мы видим, что первые сингулярные значения, вычисленные этими двумя алгоритмами SVD, очень близки. Однако второе сингулярное значение рандомизированного SVD имеет небольшую погрешность. Есть ли другой способ улучшить этот результат? И как?
Ответ положительный!
Рандомизированный SVD с степенной итерацией
Для повышения качества рандомизированного SVD можно напрямую использовать метод степенных итераций. Для получения более подробной информации об итерации мощности см. Стр. 39 в [1], а также реализацию Matlab на стр. 40.
В следующих кодах Python power_iteration() - это функция для итеративного вычисления матрицы m на k Y ( значение по умолчанию power_iter равно 3), а затем получить матрицу m на k Q с помощью QR-разложения.
import numpy as np
def power_iteration(A, Omega, power_iter = 3):
Y = A @ Omega
for q in range(power_iter):
Y = A @ (A.T @ Y)
Q, _ = np.linalg.qr(Y)
return Q
def rsvd(A, Omega):
Q = power_iteration(A, Omega)
B = Q.T @ A
u_tilde, s, v = np.linalg.svd(B, full_matrices = 0)
u = Q @ u_tilde
return u, s, v
Давайте протестируем нашу новую функцию rsvd():
np.random.seed(1000)
A = np.array([[1, 3, 2],
[5, 3, 1],
[3, 4, 5]])
rank = 2
Omega = np.random.randn(A.shape[1], rank)
u, s, v = rsvd(A, Omega)
print('Left singular vectors:')
print(u)
print()
print('Singular values:')
print(s)
print()
print('Right singular vectors:')
print(v)
print()
Результат:
Left singular vectors: [[ 0.37421757 0.28528579] [ 0.56470638 -0.82484381] [ 0.73557319 0.48810317]] Singular values: [9.34265841 3.24497775] Right singular vectors: [[ 0.57847229 0.61642675 0.53421706] [-0.73178429 0.10284774 0.67373147]]
Напомним, что:
- Особые значения SVD: 9,3427, 3,2450 и 1,0885.
- Особые значения рандомизированного SVD без степенной итерации: 9.3422 и 3.0204.
- Особые значения рандомизированного SVD с степенной итерацией: 9.3427 и 3,2450.
Как видите, рандомизированный SVD с степенной итерацией обеспечивает чрезвычайно точные сингулярные значения.
Сжатие изображения
Как упомянуто выше, можно сжать (низкоранговую) матрицу сигналов с помощью SVD или рандомизированного SVD. Фактически, способ сжатия изображения с помощью SVD довольно прост: взять SVD изображения напрямую и сохранить только доминирующие сингулярные значения и левый / правый сингулярные векторы. В терминах рандомизированного SVD мы можем сначала предопределить количество доминирующих сингулярных значений, а затем получить сингулярные значения и левые / правые сингулярные векторы с помощью рандомизированного SVD.
Для нашей оценки мы выбрали цветное изображение Лены в качестве наших данных. Размер этого изображения 256 × 256 × 3. Здесь мы строим матрицу A размером 256 × 256, выбирая только зеленый канал.
- Использование SVD напрямую
- Вместо этого используется рандомизированный SVD
Из этого эксперимента по сжатию изображений мы видим, что:
- 1) По сравнению с SVD, рандомизированный SVD также может производить точное сжатие с заданным низким рангом (здесь мы устанавливаем
rank = 50). - 2) Рандомизированный SVD удобен для вычислений. При указании
rank = 50общее время ЦП рандомизированного SVD составляет около 11,6 мс, а общее время ЦП SVD составляет 31,5 мс.
Резюме
В этом посте вы открыли для себя метод рандомизированной линейной алгебры для SVD. В частности, вы узнали:
- Основная идея рандомизированной СВД.
- Как реализовать рандомизированный SVD в Python.
У вас есть вопросы?
использованная литература
[1] Стивен Л. Брантон, Дж. Натан Куц (2019). Наука и инженерия, управляемые данными: машинное обучение, динамические системы и управление. Стр. 37–41.
[2] Н. Бенджамин Эрихсон, Сергей Воронин, Стивен Л. Брантон, Дж. Натан Куц (2016). Рандомизированное разложение матриц с использованием R. arXiv: 1608.02148. [PDF]
