Я работал над реализацией решения для разреженных неопределенных систем в Python (обсуждалось здесь), и я пытался перестроить функцию nullspace, которая использует стандартную функцию numpy svd (numpy.linalg.svd
) в кулинарной книге SciPy, используя версию scipy.sparse svd (scipy.sparse.linalg.svds
), но он выводит разные левые и правые сингулярные векторы для примеров, которые я запускал, включая матрицы:
[[1,1,0,0,0],[0,0,1,1,0],[1,1,1,1,1]]
[[1,0,1],[1,1,0],[0,1,1]]
Почему эти два решателя производят два разных выхода svd для матриц выше? Что я могу сделать для обеспечения того же выхода?
Вот пример: table - это csc_matrix
такая, что
table.todense() = matrix([[1,1,0,0,0],[0,0,1,1,0],[1,1,1,1,1]],dtype=int64)
Итак, следующие выходы кода
numpy.linalg.svd(table.todense()) =
[[ -3.64512933e-01 7.07106781e-01 -6.05912800e-01]
[ -3.64512933e-01 -7.07106781e-01 -6.05912800e-01]
[ -8.56890100e-01 2.32635116e-16 5.15499134e-01]]
-----------------------------------------------------
[ 2.58873755 1.41421356 0.54629468]
-----------------------------------------------------
[[ -4.7181e-01 -4.7181e-01 -4.7181e-01 -4.7181e-01 -3.3101e-01]
[5e-01 5e-01 -5e-01 -5e-01 6.16450329e-17]
[-1.655e-01 -1.655e-01 -1.655e-01 -1.655e-01 9.436e-01]
[5e-01 -5e-01 -5e-01 5e-01 -1.77302319e-16]
[-5e-01 5e-01 -5e-01 5e-01 2.22044605e-16]]
И следующее
scipy.sparse.linalg.svds(table,k=2)=
[[ 7.07106781e-01, -3.64512933e-01],
[ -7.07106781e-01, -3.64512933e-01],
[ 2.73756255e-18, -8.56890100e-01]]
-------------------------------------
[ 1.41421356, 2.58873755]
-------------------------------------
[[ 5e-01, 5e-01, -5e-01, -5e-01, 1.93574904e-18],
[ -4.71814e-01, -4.71814e-01, -4.71814e-01, -4.71814e-01, -3.31006e-01]]
Обратите внимание, что существует немало значений, которые перекрываются между двумя решениями. Кроме того, функция scipy.sparse.linalg.svds
не позволяет k быть больше или равна min(table.shape)
, поэтому я выбрал k = 2.
Результат в вопросе, который вы опубликовали, выглядит хорошо для меня. В вызове numpy вы вычисляете каждое сингулярное значение, а в scipy-коде вы вычисляете только верхние k сингулярных значений, и они соответствуют вершине k из вывода numpy.
Редкий верхний k svd не позволит вам рассчитать каждое исключительное значение, потому что если вы хотите это сделать, тогда вы можете просто использовать полную функцию svd.
Ниже я включил код для вас, чтобы проверить это самостоятельно. Предостережение заключается в том, что в то время как numpy и scipy full svds могут воссоздать исходную матрицу достаточно хорошо, верхний k svd не может. Это потому, что вы выбрасываете данные. Обычно это прекрасно, учитывая, что вы в порядке, будучи достаточно близко. Проблема заключается в том, что SVD, если используется с верхним k, может использоваться как приближение низкого ранга исходной матрицы, а не как замена.
Для ясности мой опыт в этом заключается в реализации python, параллельной версии этой статьи для оригинального автора, A Sparse Plus с низким показателем экспоненциальной модели языка для ограниченных сценариев ресурсов.
import numpy as np
from scipy import linalg
from scipy.sparse import linalg as slinalg
x = np.array([[1,1,0,0,0],[0,0,1,1,0],[1,1,1,1,1]],dtype=np.float64)
npsvd = np.linalg.svd(x)
spsvd = linalg.svd(x)
sptop = slinalg.svds(x,k=2)
print "np"
print "u: ", npsvd[0]
print "s: ", npsvd[1]
print "v: ", npsvd[2]
print "\n=========================\n"
print "sp"
print "u: ", spsvd[0]
print "s: ", spsvd[1]
print "v: ", spsvd[2]
print "\n=========================\n"
print "sp top k"
print "u: ", sptop[0]
print "s: ", sptop[1]
print "v: ", sptop[2]
nptmp = np.zeros((npsvd[0].shape[1],npsvd[2].shape[0]))
nptmp[np.diag_indices(np.min(nptmp.shape))] = npsvd[1]
npreconstruct = np.dot(npsvd[0], np.dot(nptmp,npsvd[2]))
print npreconstruct
print "np close? : ", np.allclose(npreconstruct, x)
sptmp = np.zeros((spsvd[0].shape[1],spsvd[2].shape[0]))
sptmp[np.diag_indices(np.min(sptmp.shape))] = spsvd[1]
spreconstruct = np.dot(spsvd[0], np.dot(sptmp,spsvd[2]))
print spreconstruct
print "sp close? : ", np.allclose(spreconstruct, x)
sptoptmp = np.zeros((sptop[0].shape[1],sptop[2].shape[0]))
sptoptmp[np.diag_indices(np.min(sptoptmp.shape))] = sptop[1]
sptopreconstruct = np.dot(sptop[0], np.dot(sptoptmp,sptop[2]))
print sptopreconstruct
print "sp top close? : ", np.allclose(sptopreconstruct, x)