scipy.sparse.linalg.

svds#

scipy.sparse.linalg.svds(A, k=6, ncv=None, tol=0, which='LM', v0=None, maxiter=None, return_singular_vectors=True, solver='arpack', rng=None, опции=None, *, random_state=None)[источник]#

Частичное разложение по сингулярным числам разреженной матрицы.

Вычислить наибольшее или наименьшее k сингулярные значения и соответствующие сингулярные векторы разреженной матрицы A. Порядок, в котором возвращаются сингулярные значения, не гарантирован.

В описаниях ниже пусть M, N = A.shape.

Параметры:

Andarray, разреженная матрица или LinearOperator

Матрица для разложения с плавающей точкой числового типа данных.

kint, по умолчанию: 6

Количество сингулярных значений и сингулярных векторов для вычисления. Должно удовлетворять 1 <= k <= kmax, где kmax=min(M, N) для solver='propack' и kmax=min(M, N) - 1 в противном случае.

ncvint, необязательный

Когда solver='arpack', это количество векторов Ланцоша, сгенерированных. См. ‘arpack’ для подробностей. Когда solver='lobpcg' или solver='propack', этот параметр игнорируется.

tolfloat, опционально

Допуск для сингулярных значений. Ноль (по умолчанию) означает машинную точность.

which{‘LM’, ‘SM’}

Который k сингулярные значения для поиска: либо наибольшие по величине ('LM'), либо наименьшие по величине ('SM') сингулярные значения.

v0ndarray, необязательно

Начальный вектор для итерации; см. документацию по конкретному методу (‘arpack’, 'lobpcg'), или ‘propack’ подробности.

maxiterint, необязательный

Максимальное количество итераций; см. документацию по конкретному методу (‘arpack’, 'lobpcg'), или ‘propack’ подробности.

return_singular_vectors{True, False, "u", "vh"}

Сингулярные значения всегда вычисляются и возвращаются; этот параметр управляет вычислением и возвратом сингулярных векторов.

True: возвращает сингулярные векторы.
False: не возвращать сингулярные векторы.
"u": если M <= N, вычислить только левые сингулярные векторы и вернуть None для правых сингулярных векторов. В противном случае вычислить все сингулярные векторы.
"vh": если M > N, вычислять только правые сингулярные векторы и возвращать None для левых сингулярных векторов. В противном случае вычислять все сингулярные векторы.

Если solver='propack', опция учитывается независимо от формы матрицы.

solver{‘arpack’, ‘propack’, ‘lobpcg’}, опционально

Используемый решатель. ‘arpack’, 'lobpcg', и ‘propack’ поддерживаются. По умолчанию: ‘arpack’.

rng{None, int, numpy.random.Generator, опционально

Если rng передается по ключевому слову, типы, отличные от numpy.random.Generator передаются в numpy.random.default_rng для создания экземпляра Generator. Если rng уже является Generator экземпляр, то предоставленный экземпляр используется. Укажите rng для повторяемого поведения функции.

Если этот аргумент передаётся по позиции или random_state передается по ключевому слову, устаревшее поведение для аргумента random_state применяется:

Если random_state равно None (или numpy.random), numpy.random.RandomState используется синглтон.
Если random_state является int, новый RandomState используется экземпляр, инициализированный с random_state.
Если random_state уже является Generator или RandomState экземпляр, тогда этот экземпляр используется.

Изменено в версии 1.15.0: В рамках SPEC-007 переход от использования numpy.random.RandomState to numpy.random.Generator, этот ключевое слово было изменено с random_state to rng. В переходный период оба ключевых слова будут продолжать работать, хотя только одно может быть указано за раз. После переходного периода вызовы функций, использующие random_state ключевое слово будет выдавать предупреждения. Поведение обоих random_state и rng описаны выше, но только rng ключевое слово должно использоваться в новом коде.

опцииdict, optional

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

Возвращает:

undarray, shape=(M, k): Унитарная матрица, имеющая левые сингулярные векторы в качестве столбцов.
sndarray, shape=(k,): Сингулярные значения.
vhndarray, shape=(k, N): Унитарная матрица, имеющая правые сингулярные векторы в качестве строк.

Примечания

Это наивная реализация, использующая ARPACK или LOBPCG в качестве решателя собственных значений матрицы A.conj().T @ A или A @ A.conj().T, в зависимости от того, какой размер меньше, с последующим методом Рэлея-Ритца как постобработкой; см. Использование нормальной матрицы, в методе Рэлея-Ритца, (2022, 19 ноя.), Википедия, https://w.wiki/4zms.

В качестве альтернативы можно вызвать решатель PROPACK.

Варианты входной матрицы A числовой тип данных может быть ограничен. Только solver="lobpcg" поддерживает все типы данных с плавающей точкой вещественные: ‘np.float32’, ‘np.float64’, ‘np.longdouble’ и комплексные: ‘np.complex64’, ‘np.complex128’, ‘np.clongdouble’. solver="arpack" поддерживает только ‘np.float32’, ‘np.float64’ и ‘np.complex128’.

Примеры

Построить матрицу A из сингулярных значений и векторов.

>>> import numpy as np
>>> from scipy import sparse, linalg, stats
>>> from scipy.sparse.linalg import svds, aslinearoperator, LinearOperator

Построение плотной матрицы A из сингулярных значений и векторов.

>>> rng = np.random.default_rng()
>>> orthogonal = stats.ortho_group.rvs(10, random_state=rng)
>>> s = [1e-3, 1, 2, 3, 4]  # non-zero singular values
>>> u = orthogonal[:, :5]         # left singular vectors
>>> vT = orthogonal[:, 5:].T      # right singular vectors
>>> A = u @ np.diag(s) @ vT

Всего с четырьмя сингулярными значениями/векторами SVD аппроксимирует исходную матрицу.

>>> u4, s4, vT4 = svds(A, k=4)
>>> A4 = u4 @ np.diag(s4) @ vT4
>>> np.allclose(A4, A, atol=1e-3)
True

Со всеми пятью ненулевыми сингулярными значениями/векторами мы можем более точно воспроизвести исходную матрицу.

>>> u5, s5, vT5 = svds(A, k=5)
>>> A5 = u5 @ np.diag(s5) @ vT5
>>> np.allclose(A5, A)
True

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

>>> np.allclose(s5, s)
True

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

>>> (np.allclose(np.abs(u5), np.abs(u)) and
...  np.allclose(np.abs(vT5), np.abs(vT)))
True

Сингулярные векторы также ортогональны.

>>> (np.allclose(u5.T @ u5, np.eye(5)) and
...  np.allclose(vT5 @ vT5.T, np.eye(5)))
True

Если существуют (почти) кратные сингулярные значения, соответствующие индивидуальные сингулярные векторы могут быть неустойчивыми, но всё инвариантное подпространство, содержащее все такие сингулярные векторы, вычисляется точно, что может быть измерено углами между подпространствами через 'subspace_angles'.

>>> rng = np.random.default_rng()
>>> s = [1, 1 + 1e-6]  # non-zero singular values
>>> u, _ = np.linalg.qr(rng.standard_normal((99, 2)))
>>> v, _ = np.linalg.qr(rng.standard_normal((99, 2)))
>>> vT = v.T
>>> A = u @ np.diag(s) @ vT
>>> A = A.astype(np.float32)
>>> u2, s2, vT2 = svds(A, k=2, rng=rng)
>>> np.allclose(s2, s)
True

Углы между отдельными точными и вычисленными сингулярными векторами могут быть не такими малыми. Для проверки используйте:

>>> (linalg.subspace_angles(u2[:, :1], u[:, :1]) +
...  linalg.subspace_angles(u2[:, 1:], u[:, 1:]))
array([0.06562513])  # may vary
>>> (linalg.subspace_angles(vT2[:1, :].T, vT[:1, :].T) +
...  linalg.subspace_angles(vT2[1:, :].T, vT[1:, :].T))
array([0.06562507])  # may vary

В отличие от углов между 2-мерными инвариантными подпространствами, которые охватывают эти векторы, которые малы для правых сингулярных векторов

>>> linalg.subspace_angles(u2, u).sum() < 1e-6
True

а также для левых сингулярных векторов.

>>> linalg.subspace_angles(vT2.T, vT.T).sum() < 1e-6
True

Следующий пример следует примеру ‘sklearn.decomposition.TruncatedSVD’.

>>> rng = np.random.default_rng()
>>> X_dense = rng.random(size=(100, 100))
>>> X_dense[:, 2 * np.arange(50)] = 0
>>> X = sparse.csr_array(X_dense)
>>> _, singular_values, _ = svds(X, k=5, rng=rng)
>>> print(singular_values)
[ 4.3221...  4.4043...  4.4907...  4.5858... 35.4549...]

Функция может быть вызвана без явного построения транспонирования входной матрицы.

>>> rng = np.random.default_rng()
>>> G = sparse.random_array((8, 9), density=0.5, rng=rng)
>>> Glo = aslinearoperator(G)
>>> _, singular_values_svds, _ = svds(Glo, k=5, rng=rng)
>>> _, singular_values_svd, _ = linalg.svd(G.toarray())
>>> np.allclose(singular_values_svds, singular_values_svd[-4::-1])
True

Наиболее эффективный по памяти сценарий - когда ни исходная матрица, ни её транспонированная версия не создаются явно. Наш пример вычисляет наименьшие сингулярные значения и векторы 'LinearOperator', созданного из функции numpy 'np.diff', используемой построчно для соответствия работе 'LinearOperator' со столбцами.

>>> diff0 = lambda a: np.diff(a, axis=0)

Давайте создадим матрицу из 'diff0' для использования только для проверки.

>>> n = 5  # The dimension of the space.
>>> M_from_diff0 = diff0(np.eye(n))
>>> print(M_from_diff0.astype(int))
[[-1  1  0  0  0]
 [ 0 -1  1  0  0]
 [ 0  0 -1  1  0]
 [ 0  0  0 -1  1]]

Матрица ‘M_from_diff0’ является двудиагональной и может быть альтернативно создана напрямую с помощью

>>> M = - np.eye(n - 1, n, dtype=int)
>>> np.fill_diagonal(M[:,1:], 1)
>>> np.allclose(M, M_from_diff0)
True

Его транспонирование

>>> print(M.T)
[[-1  0  0  0]
 [ 1 -1  0  0]
 [ 0  1 -1  0]
 [ 0  0  1 -1]
 [ 0  0  0  1]]

может рассматриваться как матрица инцидентности; см. Матрица инцидентности, (2022, 19 ноя.), Википедия, https://w.wiki/5YXU, линейного графа с 5 вершинами и 4 рёбрами. Нормальная матрица 5x5 M.T @ M следовательно, является

>>> print(M.T @ M)
[[ 1 -1  0  0  0]
 [-1  2 -1  0  0]
 [ 0 -1  2 -1  0]
 [ 0  0 -1  2 -1]
 [ 0  0  0 -1  1]]

графовый лапласиан, тогда как фактически используемая в 'svds' матрица нормального размера 4x4 M @ M.T

>>> print(M @ M.T)
[[ 2 -1  0  0]
 [-1  2 -1  0]
 [ 0 -1  2 -1]
 [ 0  0 -1  2]]

является так называемым рёберным лапласианом; см. Симметричный лапласиан через матрицу инцидентности, в Лапласиан матрицы, (2022, 19 ноя.), Википедия, https://w.wiki/5YXW.

Настройка 'LinearOperator' требует опций 'rmatvec' и 'rmatmat' умножения на транспонированную матрицу M.T, но мы хотим быть матрично-свободными для экономии памяти, поэтому зная, как M.T выглядит как, мы вручную конструируем следующую функцию для использования в rmatmat=diff0t.

>>> def diff0t(a):
...     if a.ndim == 1:
...         a = a[:,np.newaxis]  # Turn 1D into 2D array
...     d = np.zeros((a.shape[0] + 1, a.shape[1]), dtype=a.dtype)
...     d[0, :] = - a[0, :]
...     d[1:-1, :] = a[0:-1, :] - a[1:, :]
...     d[-1, :] = a[-1, :]
...     return d

Мы проверяем, что наша функция 'diff0t' для транспонирования матрицы корректна.

>>> np.allclose(M.T, diff0t(np.eye(n-1)))
True

Теперь мы настраиваем наш матрично-свободный ‘LinearOperator’ под названием ‘diff0_func_aslo’ и для проверки матричный ‘diff0_matrix_aslo’.

>>> def diff0_func_aslo_def(n):
...     return LinearOperator(matvec=diff0,
...                           matmat=diff0,
...                           rmatvec=diff0t,
...                           rmatmat=diff0t,
...                           shape=(n - 1, n))
>>> diff0_func_aslo = diff0_func_aslo_def(n)
>>> diff0_matrix_aslo = aslinearoperator(M_from_diff0)

И проверять как матрицу, так и её транспонирование в 'LinearOperator'.

>>> np.allclose(diff0_func_aslo(np.eye(n)),
...             diff0_matrix_aslo(np.eye(n)))
True
>>> np.allclose(diff0_func_aslo.T(np.eye(n-1)),
...             diff0_matrix_aslo.T(np.eye(n-1)))
True

После проверки настройки 'LinearOperator', мы запускаем решатель.

>>> n = 100
>>> diff0_func_aslo = diff0_func_aslo_def(n)
>>> u, s, vT = svds(diff0_func_aslo, k=3, which='SM')

Квадраты сингулярных чисел и сингулярные векторы известны явно; см. Чистые граничные условия Дирихле, в Собственные значения и собственные векторы второй производной, (2022, 19 ноября), Википедия, https://w.wiki/5YX6, поскольку ‘diff’ соответствует первой производной, и его матрица нормальных уравнений меньшего размера n-1 x n-1 M @ M.T представляют дискретную вторую производную с граничными условиями Дирихле. Мы используем эти аналитические выражения для проверки.

>>> se = 2. * np.sin(np.pi * np.arange(1, 4) / (2. * n))
>>> ue = np.sqrt(2 / n) * np.sin(np.pi * np.outer(np.arange(1, n),
...                              np.arange(1, 4)) / n)
>>> np.allclose(s, se, atol=1e-3)
True
>>> np.allclose(np.abs(u), np.abs(ue), atol=1e-6)
True