Интерполяционное матричное разложение (scipy.linalg.interpolative)#

Добавлено в версии 0.13.

Изменено в версии 1.15.0: Базовые алгоритмы были портированы на Python из оригинального кода на Fortran77. Подробности см. в ссылках ниже.

Интерполяционное разложение (ID) матрицы \(A \in \mathbb{C}^{m \times n}\) ранга \(k \leq \min \{ m, n \}\) является факторизацией

\[A \Pi = \begin{bmatrix} A \Pi_{1} & A \Pi_{2} \end{bmatrix} = A \Pi_{1} \begin{bmatrix} I & T \end{bmatrix},\]

где \(\Pi = [\Pi_{1}, \Pi_{2}]\) является матрицей перестановки с \(\Pi_{1} \in \{ 0, 1 \}^{n \times k}\), т.е., \(A \Pi_{2} = A \Pi_{1} T\). Это может быть эквивалентно записано как \(A = BP\), где \(B = A \Pi_{1}\) и \(P = [I, T] \Pi^{\mathsf{T}}\) являются скелет и интерполяционные матрицы, соответственно.

Если \(A\) не имеет точного ранга \(k\), тогда существует аппроксимация в форме ID такая, что \(A = BP + E\), где \(\| E \| \sim \sigma_{k + 1}\) имеет порядок \((k + 1)\)-е наибольшее сингулярное число \(A\). Обратите внимание, что \(\sigma_{k + 1}\) является наилучшей возможной ошибкой для ранга-\(k\) аппроксимация и, фактически, достигается с помощью сингулярного разложения (SVD) \(A \approx U S V^{*}\), где \(U \in \mathbb{C}^{m \times k}\) и \(V \in \mathbb{C}^{n \times k}\) имеют ортонормированные столбцы и \(S = \mathop{\mathrm{diag}} (\sigma_{i}) \in \mathbb{C}^{k \times k}\) является диагональной с неотрицательными элементами. Основные преимущества использования ID вместо SVD заключаются в следующем:

  • его дешевле построить;

  • он сохраняет структуру \(A\); и

  • эффективнее вычислять с учетом тождественной подматрицы \(P\).

Подпрограммы#

Основная функциональность:

interp_decomp(A, eps_or_k[, rand, rng])

Вычислить ID матрицы.

reconstruct_matrix_from_id(B, idx, proj)

Восстановить матрицу из её ID.

reconstruct_interp_matrix(idx, proj)

Восстановить матрицу интерполяции из ID.

reconstruct_skel_matrix(A, k, idx)

Восстановить скелетную матрицу из ID.

id_to_svd(B, idx, proj)

Преобразовать ID в SVD.

svd(A, eps_or_k[, rand, rng])

Вычислите SVD матрицы через ID.

estimate_spectral_norm(A[, its, rng])

Оценка спектральной нормы матрицы рандомизированным степенным методом.

estimate_spectral_norm_diff(A, B[, its, rng])

Оценка спектральной нормы разности двух матриц рандомизированным степенным методом.

estimate_rank(A, eps[, rng])

Оценить ранг матрицы с заданной относительной точностью с использованием рандомизированных методов.

Следующие вспомогательные функции устарели и будут удалены в SciPy 1.17.0:

seed([seed])

Эта функция исторически использовалась для установки начального значения алгоритмов рандомизации, используемых в scipy.linalg.interpolative функции, написанные на Fortran77.

rand(*shape)

Эта функция исторически использовалась для генерации равномерно распределенных случайных чисел для алгоритмов рандомизации, используемых в scipy.linalg.interpolative функции, написанные на Fortran77.

Ссылки#

Этот модуль использует алгоритмы, найденные в программном пакете ID [R5a82238cdab4-1] Мартинссоном, Рохлиным, Школьницким и Тайгертом, которая представляет собой библиотеку на Fortran для вычисления ID с использованием различных алгоритмов, включая подход rank-revealing QR [R5a82238cdab4-2] и более современные рандомизированные методы, описанные в [R5a82238cdab4-3], [R5a82238cdab4-4], и [R5a82238cdab4-5].

Мы рекомендуем пользователю также обратиться к документации для Пакет ID.

[R5a82238cdab4-1]

P.G. Martinsson, V. Rokhlin, Y. Shkolnisky, M. Tygert. “ID: a software package for low-rank approximation of matrices via interpolative decompositions, version 0.2.” http://tygert.com/id_doc.4.pdf.

[R5a82238cdab4-2]

H. Cheng, Z. Gimbutas, P.G. Martinsson, V. Rokhlin. “On the compression of low rank matrices.” SIAM J. Sci. Comput. 26 (4): 1389–1404, 2005. DOI:10.1137/030602678.

[R5a82238cdab4-3]

E. Liberty, F. Woolfe, P.G. Martinsson, V. Rokhlin, M. Tygert. “Randomized algorithms for the low-rank approximation of matrices.” Proc. Natl. Acad. Sci. U.S.A. 104 (51): 20167–20172, 2007. DOI:10.1073/pnas.0709640104.

[R5a82238cdab4-4]

P.G. Martinsson, V. Rokhlin, M. Tygert. «Рандомизированный алгоритм для разложения матриц.» Appl. Comput. Harmon. Anal. 30 (1): 47–68, 2011. DOI:10.1016/j.acha.2010.02.003.

[R5a82238cdab4-5]

F. Woolfe, E. Liberty, V. Rokhlin, M. Tygert. «Быстрый рандомизированный алгоритм для аппроксимации матриц». Appl. Comput. Harmon. Anal. 25 (3): 335–366, 2008. DOI:10.1016/j.acha.2007.12.002.

Учебник#

Инициализация#

Первый шаг — импортировать scipy.linalg.interpolative выполнив команду:

>>> import scipy.linalg.interpolative as sli

Теперь давайте построим матрицу. Для этого рассмотрим матрицу Гильберта, которая, как известно, имеет низкий ранг:

>>> from scipy.linalg import hilbert
>>> n = 1000
>>> A = hilbert(n)

Мы также можем сделать это явно через:

>>> import numpy as np
>>> n = 1000
>>> A = np.empty((n, n), order='F')
>>> for j in range(n):
...     for i in range(n):
...         A[i,j] = 1. / (i + j + 1)

Обратите внимание на использование флага order='F' в numpy.empty. Это создает матрицу в порядке, совместимом с Fortran, что важно для избежания копирования данных при передаче в бэкенд.

Затем мы определяем процедуры умножения для матрицы, рассматривая её как scipy.sparse.linalg.LinearOperator:

>>> from scipy.sparse.linalg import aslinearoperator
>>> L = aslinearoperator(A)

Это автоматически настраивает методы, описывающие действие матрицы и её сопряжённой на вектор.

Вычисление ID#

У нас есть несколько вариантов алгоритмов для вычисления ID. Они в основном делятся по двум дихотомиям:

  1. как матрица представлена, т.е. через её элементы или через её действие на вектор; и

  2. приближать ли её к фиксированной относительной точности или к фиксированному рангу.

Мы последовательно рассмотрим каждый вариант ниже.

Во всех случаях ID представлен тремя параметрами:

  1. ранг k;

  2. массив индексов idx; и

  3. коэффициенты интерполяции proj.

ID задаётся соотношением np.dot(A[:,idx[:k]], proj) == A[:,idx[k:]].

Из элементов матрицы#

Сначала мы рассматриваем матрицу, заданную через её элементы.

Чтобы вычислить ID с фиксированной точностью, введите:

>>> eps = 1e-3
>>> k, idx, proj = sli.interp_decomp(A, eps)

где eps < 1 является желаемой точностью.

Для вычисления ID до фиксированного ранга используйте:

>>> idx, proj = sli.interp_decomp(A, k)

где k >= 1 является желаемым рангом.

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

>>> k, idx, proj = sli.interp_decomp(A, eps, rand=False)

и:

>>> idx, proj = sli.interp_decomp(A, k, rand=False)

соответственно.

Из действия матрицы#

Теперь рассмотрим матрицу, заданную через её действие на вектор как scipy.sparse.linalg.LinearOperator.

Чтобы вычислить ID с фиксированной точностью, введите:

>>> k, idx, proj = sli.interp_decomp(L, eps)

Для вычисления ID до фиксированного ранга используйте:

>>> idx, proj = sli.interp_decomp(L, k)

Эти алгоритмы рандомизированы.

Восстановление ID#

Приведённые выше процедуры ID не выводят матрицы скелета и интерполяции явно, а вместо этого возвращают соответствующую информацию в более компактной (и иногда более полезной) форме. Чтобы построить эти матрицы, напишите:

>>> B = sli.reconstruct_skel_matrix(A, k, idx)

для скелетной матрицы и:

>>> P = sli.reconstruct_interp_matrix(idx, proj)

для матрицы интерполяции. Тогда ID-аппроксимация может быть вычислена как:

>>> C = np.dot(B, P)

Это также может быть построено напрямую с использованием:

>>> C = sli.reconstruct_matrix_from_id(B, idx, proj)

без необходимости сначала вычислять P.

Альтернативно, это также можно сделать явно, используя:

>>> B = A[:,idx[:k]]
>>> P = np.hstack([np.eye(k), proj])[:,np.argsort(idx)]
>>> C = np.dot(B, P)

Вычисление SVD#

ID может быть преобразован в SVD с помощью команды:

>>> U, S, V = sli.id_to_svd(B, idx, proj)

Аппроксимация SVD тогда:

>>> approx = U @ np.diag(S) @ V.conj().T

SVD также можно вычислить «с нуля», объединив оба шага ID и преобразования в одну команду. В соответствии с различными алгоритмами ID выше, существуют соответствующие различные алгоритмы SVD, которые можно использовать.

Из элементов матрицы#

Сначала мы рассматриваем алгоритмы сингулярного разложения для матрицы, заданной своими элементами.

Для вычисления SVD с фиксированной точностью введите:

>>> U, S, V = sli.svd(A, eps)

Для вычисления SVD с фиксированным рангом используйте:

>>> U, S, V = sli.svd(A, k)

Оба алгоритма используют случайную выборку; для детерминированных версий укажите ключевое слово rand=False как выше.

Из действия матрицы#

Теперь рассмотрим матрицу, заданную в терминах ее действия на вектор.

Для вычисления SVD с фиксированной точностью введите:

>>> U, S, V = sli.svd(L, eps)

Для вычисления SVD с фиксированным рангом используйте:

>>> U, S, V = sli.svd(L, k)

Вспомогательные процедуры#

Также доступно несколько служебных подпрограмм.

Для оценки спектральной нормы матрицы используйте:

>>> snorm = sli.estimate_spectral_norm(A)

Этот алгоритм основан на рандомизированном степенном методе и поэтому требует только произведений матрица-вектор. Количество итераций можно задать с помощью ключевого слова its (по умолчанию: its=20). Матрица интерпретируется как scipy.sparse.linalg.LinearOperator, но также допустимо предоставлять его как numpy.ndarray, в этом случае он тривиально преобразуется с помощью scipy.sparse.linalg.aslinearoperator.

Тот же алгоритм также может оценить спектральную норму разности двух матриц A1 и A2 следующим образом:

>>> A1, A2 = A**2, A
>>> diff = sli.estimate_spectral_norm_diff(A1, A2)

Это часто полезно для проверки точности приближения матрицы.

Некоторые процедуры в scipy.linalg.interpolative требует оценки ранга матрицы также. Это можно сделать с помощью:

>>> k = sli.estimate_rank(A, eps)

или:

>>> k = sli.estimate_rank(L, eps)

в зависимости от представления. Параметр eps управляет определением численного ранга.

Наконец, генерация случайных чисел, необходимая для всех рандомизированных процедур, может быть контролируема путём предоставления псевдослучайных генераторов NumPy с фиксированным сидом. См. numpy.random.Generator и numpy.random.default_rng для получения дополнительной информации.

Замечания#

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