Странное поведение при использовании Arpack EigenLibrary Wrapper для маленьких матриц

Я хочу решить следующий обобщенный EVP, используя ArpackWrapper библиотеки Eigen:

введите здесь описание изображения

К_е - СПД. Обычно K_g неопределенное и единственное число, но для этого MVP оно просто неопределенное. Кроме того, меня интересуют наименьшие собственные значения. Для больших систем 300kx300k я получил приемлемые результаты, но для этого небольшого примера результаты кажутся странными. Первые 4 результирующих собственных значения читаются

 0.8987
-0.720851 
 0.607632 
 0.729297

Для менее запрошенных собственных значений solverEig.info() возвращает Eigen::NoConvergence.

Если я использую Maple для вычисления собственных значений, которые я получаю

 522.991427951073 
-175.66558721639944
  66.23707710214939
  -7.756864956770603
 355.6461914072188

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

Поэтому мой вопрос: почему это не дает правильного результата? Это проблема arpack или собственной оболочки? Или, скорее всего, мое неправильное использование/понимание arpack или собственной оболочки.

Версии:

Эйген 3.3.7

Arpack https://github.com/opencollab/arpack-ng приводит к одному и тому же результату поведение

Код:

#include <iostream>
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include "eigen3/unsupported/Eigen/ArpackSupport"
#include "eigen3/unsupported/Eigen/SparseExtra"
using namespace std;
int main() {

    Eigen::SparseMatrix<double> ke; 
    Eigen::SparseMatrix<double> kg; 

    Eigen::loadMarket(ke,"ke.txt"); //5x5 Matrix
    Eigen::loadMarket(kg,"kg.txt"); //5x5 Matrix


    Eigen::ArpackGeneralizedSelfAdjointEigenSolver<Eigen::SparseMatrix<double>> solverEig;

    solverEig.compute(ke,-kg,4,"SM",Eigen::ComputeEigenvectors);


    cout<<solverEig.eigenvalues()<<endl;
    return 0;
}

Файлы:

Ke.txt

%%MatrixMarket matrix coordinate  real general
5 5 19
1 1 3.2621670111997303820317029021680355072021484375000000000000000000e+01
2 1 1.4310835055998653686515353911090642213821411132812500000000000000e+01
3 1 -2.0000000000000000000000000000000000000000000000000000000000000000e+00
5 1 -1.4310835055998653686515353911090642213821411132812500000000000000e+01
1 2 1.4310835055998653686515353911090642213821411132812500000000000000e+01
2 2 8.7155417527999333060506614856421947479248046875000000000000000000e+01
3 2 0.0000000000000000000000000000000000000000000000000000000000000000e+00
5 2 -7.1554175279993268432576769555453211069107055664062500000000000000e+00
1 3 -2.0000000000000000000000000000000000000000000000000000000000000000e+00
2 3 0.0000000000000000000000000000000000000000000000000000000000000000e+00
3 3 8.2000000000000000000000000000000000000000000000000000000000000000e+01
4 3 -8.0000000000000000000000000000000000000000000000000000000000000000e+01
3 4 -8.0000000000000000000000000000000000000000000000000000000000000000e+01
4 4 8.0000000000000000000000000000000000000000000000000000000000000000e+01
5 4 0.0000000000000000000000000000000000000000000000000000000000000000e+00
1 5 -1.4310835055998653686515353911090642213821411132812500000000000000e+01
2 5 -7.1554175279993268432576769555453211069107055664062500000000000000e+00
4 5 0.0000000000000000000000000000000000000000000000000000000000000000e+00
5 5 8.7155417527999333060506614856421947479248046875000000000000000000e+01

кг.txt

%%MatrixMarket matrix coordinate  real general
5 5 19
1 1 -2.1389117730512491322159007722802925854921340942382812500000000000e-01
2 1 0.0000000000000000000000000000000000000000000000000000000000000000e+00
3 1 1.8999822924102233168142106478626374155282974243164062500000000000e-01
5 1 0.0000000000000000000000000000000000000000000000000000000000000000e+00
1 2 0.0000000000000000000000000000000000000000000000000000000000000000e+00
2 2 -2.1389117730512491322159007722802925854921340942382812500000000000e-01
3 2 0.0000000000000000000000000000000000000000000000000000000000000000e+00
5 2 -9.3693312577685552988704387189500266686081886291503906250000000000e-02
1 3 1.8999822924102233168142106478626374155282974243164062500000000000e-01
2 3 0.0000000000000000000000000000000000000000000000000000000000000000e+00
3 3 -3.8974822924103957877406401166808791458606719970703125000000000000e-01
4 3 1.9975000000001727484821856251073768362402915954589843750000000000e-01
3 4 1.9975000000001727484821856251073768362402915954589843750000000000e-01
4 4 1.8161606213113259955527212241577217355370521545410156250000000000e-01
5 4 0.0000000000000000000000000000000000000000000000000000000000000000e+00
1 5 0.0000000000000000000000000000000000000000000000000000000000000000e+00
2 5 -9.3693312577685552988704387189500266686081886291503906250000000000e-02
4 5 0.0000000000000000000000000000000000000000000000000000000000000000e+00
5 5 4.7505937470883541351440726430155336856842041015625000000000000000e-01

c++ c++17 arpack eigenvalue eigen3
person Alex M    schedule 31.03.2020    source источник
comment
не ваш вопрос, но вы проверили собственные векторы, Ke v + Kg v λ ~ 0?   -  person denis    schedule 09.04.2020
comment
собственные векторы также не совпадают   -  person Alex M    schedule 11.08.2020

Ответы (2)


arrow_upward
1
arrow_downward

Возможно, в вашем случае параметр NCV ARPACK по умолчанию равен NEV (nb собственного вектора для вычисления) или размеру матрицы N, которые оба малы, хотя рекомендуется, чтобы этот параметр был не меньше 20 (Я бы поставил по крайней мере 70 или 100 на самом деле ..). Смотрите эти ссылки:

https://www.caam.rice.edu/software/ARPACK/UG/node136.html (NCV, NEV и примечание 4)

https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html#scipy.sparse.linalg.eigs (описание ncv)

person Thibault Cimic    schedule 12.08.2020
comment
ncv Количество сгенерированных векторов Ланцоша ncv должно быть больше k; рекомендуется, чтобы ncv › 2*k. По умолчанию: мин(n, макс(2*k + 1, 20)). Я тоже читал об этом. Следовательно, в случае матриц размера 5 он выберет ncv равным 5. Таким образом, ожидается, что результат будет необоснованным, или мне все же следует ожидать разумных результатов? Если нет, то решением будет проверка размера матрицы и переключения, например. к DenseGeneralEigenSolver Eigen? - person Alex M; 13.08.2020
comment
В том-то и дело, что да, ncv = 5 слишком мало, чтобы такие итерационные методы давали хорошие результаты, лучше разложить на множители явно, не так дорого. Вы можете попробовать заставить ncv быть 50 или 100, но не уверены, что вы получите. Вы можете проверить книгу Юсефа Саада о проблеме собственных значений: www-users.cs. umn.edu/~saad/books.html - person Thibault Cimic; 13.08.2020

arrow_upward
0
arrow_downward

Странно, что вы получаете настоящие собственные значения, а не сложную пару — возможно, вы запускали eigsh для симметричных задач, которые рассматривают только нижние половины?

from scipy.sparse.linalg import eigs  # arpack

v0 = np.ones( K.shape[0] )  # else random

evals, V = eigs( A=K, M=M, k=3, tol=0, v0=v0 )  # k < n-1
print( "eigs evals:", evals )
eigcheck( K, evals, V, M=M )

K: 
 [[32.6217 14.3108 -2 0 -14.3108]
 [14.3108 87.1554 0 0 -7.15542]
 [-2 0 82 -80 0]
 [0 0 -80 80 0]
 [-14.3108 -7.15542 0 0 87.1554]]
M: 
 [[-0.213891 0 0.189998 0 0]
 [0 -0.213891 0 0 -0.0936933]
 [0.189998 0 -0.389748 0.19975 0]
 [0 0 0.19975 0.181616 0]
 [0 -0.0936933 0 0 0.475059]]
eigs evals: [ -1301.07 -2136.22j   -1301.07 +2136.22j   36.0965 -2477.99j ]
eigcheck : |Av - λv|_max [2e-13 2e-13 2e-13]

Не могли бы вы указать мне на подобные большие проблемы в Интернете? Спасибо

person denis    schedule 12.08.2020