c++ - Armadillo eigs_sym : decomposition failed

Question

我正在 Visual Studio 中编写 C++ 代码并使用 Armadillo 7.900.1

我没有运气让 eigs_sym 函数在 Armadillo 下工作(我使用的是 Armadillo 附带的 lapack 和 blas 版本)。典型的错误信息如下:

error: lapack::stedc(): failed to compute all eigenvalues
warning: eigs_sym(): decomposition failed

产生这个输出的代码是:

sp_mat T(locations, values);
arma::eigs_sym(eigval, eigvec, T, num_eigs_wanted, "sm", tol);

在哪里

T is a 480,000x480,000 sparse matrix. 

如果 T 很小(即 2000x2000)并且 tol 很高(tol = 5)，我可以让代码正常工作。但是一旦 T 大约为 20000x20000，那么无论使用什么值的 tol 或 num_eigs_wanted 都会出现上述错误。 (我显然改变了“值”和位置“来改变 T 的大小)。

矩阵 T 是对称实数正定矩阵。

当我在 Matlab 中调用 eigs 函数时，完全相同的矩阵没有问题。

有什么想法吗？在我看来， Armadillo 中的 eigs_sym 似乎坏了……人们使用过任何替代品吗？

干杯

Answer 1

不是答案，但评论时间太长了。我可以用下面的代码重现。该矩阵只是一个对角矩阵，对角线上的值从1开始递增。特征值显然是1,2,3,...，矩阵是正定对称的。

#include <iostream>
#include <armadillo>

int main()
{
  constexpr int N = 20000;
  arma::umat locations(2,N);
  arma::vec values(N);

  for (int i = 0; i < N; ++i)
  {
    locations(0,i) = i;
    locations(1,i) = i;
    values(i) = i+1;
  }

  arma::sp_mat T(locations, values);

  int num_eigs_wanted = 1;
  arma::vec eigval(num_eigs_wanted);
  arma::mat eigvec(N,num_eigs_wanted);
  double tol = 1e-6;
  arma::eigs_sym(eigval, eigvec, T, num_eigs_wanted, "sm", tol);

  std::cout << eigval << '\n';
}

编译器标志是

clang++-5.0 -std=c++11 test.cpp -larmadillo

我将 Armadillo 与 MKL 后端一起使用，我收到:

warning: eigs_sym(): decomposition failed
[matrix size: 0x1]

---

这似乎是 Armadillo 的问题，因为当我使用 Eigen 的不受支持的 Arpack 支持时，它工作得很好。

#include <iostream>
#include <Eigen/Sparse>
#include <unsupported/Eigen/ArpackSupport>

int main()
{
  constexpr int N = 20000;

  std::vector< Eigen::Triplet<double> > triplets;
  triplets.reserve(N);
  for (int i = 0; i < N; ++i)
    triplets.push_back( {i,i,i+1} );

  Eigen::SparseMatrix<double> T(N,N);
  T.setFromTriplets(triplets.begin(), triplets.end());

  int num_eigs_wanted = 1;
  double tol = 1e-6;
  Eigen::ArpackGeneralizedSelfAdjointEigenSolver< Eigen::SparseMatrix<double> > eigs_sym;
  eigs_sym.compute(T, num_eigs_wanted, "SM", Eigen::ComputeEigenvectors, tol);
  std::cout << eigs_sym.eigenvalues() << '\n';
}

您可能还想看看 https://github.com/yixuan/spectra/