Использование альтернативного драйвера LAPACK в методе svd numpy?

Я использую numpy.svd для вычисления сингулярных разложений плохо обусловленных матриц. В некоторых особых случаях svd не сходится и вызывает Linalg.Error. Я провел некоторое исследование и обнаружил, что numpy использует процедуру DGESDD из LAPACK. Стандартная реализация имеет жестко заданный лимит итераций в 35 или около того. Если я пытаюсь разложить ту же матрицу в Matlab, все работает нормально, и я думаю, что на это есть две причины: 1. Matlab использует DGESVD вместо DGESDD, который в целом кажется более надежным. 2. Matlab использует в подпрограмме ограничение на число итераций, равное 75. (Они изменили его в исходниках и перекомпилировали.)

Теперь возникает вопрос: есть ли простой способ изменить используемый бэкэнд в numpy с DGESDD на DGESVD без необходимости изменять источник numpy?

Заранее спасибо Миша


numpy exception lapack svd
person Mischa Obrecht    schedule 21.05.2012    source источник
comment
Это маловероятно, но вы можете написать свой собственный модуль f2py, вызывающий DGESVD.   -  person Vladimir F    schedule 21.05.2012
comment
Я надеялся избежать этого, но, похоже, мне все равно придется это сделать...   -  person Mischa Obrecht    schedule 21.05.2012

Ответы (2)


arrow_upward
1
arrow_downward

Мне помогло вычислить только "экономичный размер" SVD эта матрица X:

U,S,V = np.linalg.svd(X, full_matrices=False)
person theV0ID    schedule 03.10.2015

arrow_upward
0
arrow_downward

Я немного опоздал, но, возможно, это поможет кому-то еще...

У меня была похожая проблема в julia.

Я нашел этот подход от список справки R, который должен работать в любой среде, использующей библиотеку lapack:

По сути, если svd(M) терпит неудачу, попробуйте svd(M') и соответствующим образом поменяйте местами полученные U,V.

Вот как я делаю это в Джулии:

try
  U,S,V = svd( E_restricted )
  failed = false
catch
  failed = true
end
if failed
  # try it with matrix transposed
  try
    V,S,U = svd( E_restricted' )
    failed = false
  catch
    failed = true
  end
end
if failed
  error("ERROR: svd(E) and svd(E') failed!")
end
person danny    schedule 05.04.2013