Броненосец: eigs_gen для наименьшего собственного значения

Я использую eigs_gen броненосца, чтобы найти наименьшее алгебраическое собственное значение разреженной матрицы.

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

#include <iostream>
#include <armadillo>

using namespace std;
using namespace arma;


int
main(int argc, char** argv)
  {
  cout << "Armadillo version: " << arma_version::as_string() << endl;

  sp_mat A(5,5);

  A(1,2) = -1;
  A(2,1) = -1;
  A(3,4) = -1;
  A(4,3) = -1;

  cx_vec eigval;
  cx_mat eigvec;

  eigs_gen(eigval, eigvec, A, 1, "sr");  // find smallest eigenvalue ---> INCORRECT RESULTS
  eigval.print("Smallest real eigval:");

  eigs_gen(eigval, eigvec, A, 2, "sr");  // find 2 smallest eigenvalues ---> ALMOST CORRECT RESULTS 
  eigval.print("Two smallest real eigvals:");

  return 0;
  }

Моя команда компиляции:

 g++ file.cpp -o file.exe -O2 -I/path-to-armadillo/armadillo-4.600.3/include -DARMA_DONT_USE_WRAPPER -lblas -llapack -larpack

Результат:

Armadillo version: 4.600.3 (Off The Reservation)
Smallest real eigval:
    (+1.000e+00,+0.000e+00)
Two smallest real eigvals:
    (-1.000e+00,+0.000e+00)
    (-1.164e-17,+0.000e+00)

Любая идея о том, почему это происходит и как это преодолеть, приветствуется.

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

Обновление: добавлена ​​конструкция тестовой матрицы, которая после внесения Райаном изменений, включивших в библиотеку параметр "sa", похоже, не сходится:

    #define ARMA_64BIT_WORD
    #include <armadillo>
    #include <iostream>
    #include <vector>
    #include <stdio.h>

    using namespace arma;
    using namespace std;

    int main(){

    size_t l(3), ls(l*l*l);
    sp_mat A = sprandn<sp_mat>(ls, ls, 0.01);
    sp_mat B = A.t()*A;

    vec eigval;
    mat eigvec;
    eigs_sym(eigval, eigvec, B, 1, "sa");

    return 0;

    }

Интересующие размеры матрицы намного больше, например. ls = 8000 - 27000, и это не совсем матрица, построенная здесь, но я предполагаю, что проблема должна быть той же.


sparse-matrix armadillo eigenvalue
person MaviPranav    schedule 20.01.2015    source источник

Ответы (2)


arrow_upward
3
arrow_downward

Я считаю, что проблема здесь в том, что вы используете eigs_gen() (который вызывает DNAUPD) на симметричной матрице. ARPACK отмечает, что DNAUPD не предназначен для симметричных матриц, но не уточняет, что произойдет, если вы все равно будете использовать симметричные матрицы:

ПРИМЕЧАНИЕ. Если линейный оператор «OP» является вещественным и симметричным относительно вещественной положительной полуопределенной симметричной матрицы B, т. е. B*OP = (OP')*B, то вместо него следует использовать подпрограмму ssaupd.

(из http://www.mathkeisan.com/usersguide/man/dnaupd.html )

Я изменил внутренний код Armadillo, чтобы передать «sa» (наименьшее алгебраическое значение) в вызовы ARPACK в eigs_sym() (sp_auxlib_meat.hpp), и мне удалось получить правильные собственные значения. Я отправил исправление, чтобы сделать поддержку "sa" и "la" доступной для eigs_sym(), что, я думаю, должно решить вашу проблему после выпуска новой версии (или в какой-то момент в будущем).

person ryan    schedule 27.01.2015
comment
Спасибо, на самом деле это очень желанное и существенное улучшение библиотеки. - person MaviPranav; 18.02.2015
comment
на самом деле, похоже, есть проблемы с конвергенцией для опции sa, отсутствующие для опции lm. Я включил обновление с кодом, показывающим конструкцию тестовой матрицы, которая, похоже, не сходится для меня. - person MaviPranav; 20.02.2015
comment
Я бы предположил, что проблема здесь не в том, что Armadillo обертывает ARPACK, а в численных методах, которые использует ARPACK. Лучшее исследование свойств вашей входной матрицы и того, подходит ли неявно перезапущенный метод Арнольди для вашей конкретной матрицы, вероятно, следующее, что нужно сделать. - person ryan; 24.02.2015

arrow_upward
1
arrow_downward

Проблема с повторяющимися собственными значениями; если я изменю первые два матричных элемента на

A(1,2) = -1.00000001;
A(2,1) = -1.00000001;

ожидаемые результаты получены.

person MaviPranav    schedule 21.01.2015