2016-04-25 77 views
2

我最近安装了犰狳并尝试了稀疏矩阵的特征值问题。不幸的是,分解失败是参数'N'(下面的代码)太大e.q.我想知道这里发生了什么。矩阵不是很复杂 - 它具有对角结构。对于太大的矩阵,特征值分解失败,armadillo eigs_sym()

UPDATE

数学也有这个矩阵的问题。它告诉我Arnoldi算法不会收敛。也许我需要在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"); 

}

回答

0

对于大小为2000的小矩阵左右,它往往是时间更容易找到所有的特征向量,因为这些方法都是不敏感接近奇异矩阵。

我换成你的代码

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

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

和编译的程序只需要不到一秒钟就解决了我的已故2015年的MacBook Pro。

+0

前段时间我已经解决了这个问题。问题在于ARPACK库。更具体地说,Armadillo库不允许你指定Krylov基础的大小和最大迭代次数,就像原始的ARPACK库一样。该矩阵的最小特征值之间的距离对于阶数5的基础来说太小(默认Armadillo 2 * nev + 1)。我在Fortran中编写了类似的程序,对于矩阵大小N = 10000,我需要90阶的基本大小以便使用Arnoldi迭代进行收敛。 Mathematica 10.0还允许您指定基向量数和最大迭代次数。有用。 – Bociek