Разложение по собственным значениям не работает с броненосцем eigs_sym() для слишком больших матриц

Я недавно установил броненосец и попробовал проблему собственных значений для разреженных матриц. К сожалению, разложение не выполняется, если параметр «N» (код ниже) слишком велик, например. 1000. Интересно, что здесь происходит. Матрица не очень сложная - она ​​имеет диагональную структуру.

ОБНОВЛЕНИЕ

У Mathematica тоже есть проблемы с этой матрицей. Это говорит мне, что алгоритм Арнольди не сходится. Может быть, мне нужно вручную указать некоторые параметры в подпрограммах Arnoldi arpack, чтобы обеспечить сходимость?

Вот мой код:

#include <armadillo>

int main ()
{
    double N = 1000.0;

    // create matrix
    int kmin = 0;
    int kmax = static_cast<int>( std::floor( N/2.0 ) );
    int dim = (kmax - kmin) + 1;

    // locations and values in sparse matrices
    arma::umat hc_locations (2, 3*dim-2);
    arma::vec hc_values (3*dim-2);

    // diagonal part
    for (int k=0; k<dim; k++)
    {
        hc_locations (0,k) = k;
        hc_locations (1,k) = k;
        hc_values (k) = 2.0/N*static_cast<double>(kmin + k)*( 2.0*( N-2.0*static_cast<double>(k + kmin) ) - 1.0 );

    }
    // upper and lower diagonal
    for (int k=0; k<dim-1; k++)
    {
        hc_locations (0,k+dim) = k;
        hc_locations (1,k+dim) = k+1;
        hc_values (k+dim) = 2.0/N*std::sqrt( ( static_cast<double>(k+1+kmin) ) *
                                             ( static_cast<double>(k+1+kmin) ) *
                                             ( N - static_cast<double>(2*(k+1+kmin)) + 1.0 ) * 
                                             ( N - static_cast<double>(2*(k+1+kmin)) + 2.0 ) );

        hc_locations (0, k+2*dim-1) = k+1;
        hc_locations (1, k+2*dim-1) = k;
        hc_values (k+2*dim-1) = 2.0/N*std::sqrt ( ( static_cast<double>(k+1+kmin) ) * 
                                                   ( static_cast<double>(k+1+kmin) ) *
                                                   ( N - static_cast<double>(2*(k+kmin)) ) *
                                                   ( N - static_cast<double>(2*(k+kmin)) - 1.0 ) );
    }

    arma::sp_mat Ham(hc_locations, hc_values);

    // eigenvalue problem
    arma::vec eigval;
    arma::mat eigvec;

    arma::eigs_sym( eigval, eigvec, Ham, 3, "sa");

}


c++ armadillo
person Bociek    schedule 25.04.2016    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Для небольших матриц размером 2000 или около того часто бывает проще найти все собственные значения и собственные векторы, поскольку эти методы менее восприимчивы к матрицам, близким к сингулярным.

я заменил твой код

arma::eigs_sym( eigval, eigvec, Ham, 3, "sa");

с

arma::mat fullMat = arma::mat(Ham);
arma::eig_sym( eigval, eigvec, fullMat);

и скомпилированная программа на моем Macbook Pro конца 2015 года решается менее чем за секунду.

person anuragm    schedule 02.06.2016
comment
Я решил проблему некоторое время назад. Проблема заключается в библиотеке ARPACK. В частности, библиотека Armadillo не позволяет указать размер базиса Крылова и максимальное количество итераций, как это делает исходная библиотека ARPACK. Расстояние между наименьшими собственными значениями этой матрицы слишком мало для базиса порядка 5 (по умолчанию Armadillo 2*nev+1). Я написал аналогичную программу на фортране, и для размера матрицы N = 10000 мне нужен базовый размер порядка 90, чтобы сходиться с использованием итерации Арнольди. Также Mathematica 10.0 позволяет указать количество базисных векторов и максимальное количество итераций. Оно работает. - person Bociek; 05.06.2016